Matematyka.pl

Kod: Zaznacz cały

#include<cmath>

void swap(double& a,double& b){

double tmp;

tmp=a;

a=b;

b=tmp;

}

class matrix{

private:

int row,col;

double** data;

public:

matrix();

matrix(int,int);

~matrix();

int getrow();

int getcol();

void setrow(int);

void setcol(int);

void setdata(double**);

void getdata(double**);

matrix swapcols(int,int);

matrix swaprows(int,int);

matrix concat(matrix);

matrix copy();

matrix transpose();

matrix add(matrix);

matrix multiply(double);

matrix multiply(matrix);

matrix inverse(double,int&);

double det(double);

int rank(double);

};





matrix::matrix(){

int i;

data=new  double*[col+1];

for(i=0;i<=col;i++);

data[i]=new double[row+1];

}



matrix::matrix(int r,int c){

int i;

data=new double*[c+1];

for(i=0;i<=c;i++)

data[i]=new double[r+1];

row=r;

col=c;

}



matrix::~matrix(){

int i;

for(i=1;i<=col;i++)

delete data[i];

}



int matrix::getrow(){

return row;

}



int matrix::getcol(){

return col;

}



void matrix::setrow(int r){

row=r;

}



void matrix::setcol(int c){

col=c;

}

void matrix::setdata(double** d){

int i,j;

for(i=1;i<=row;i++)

for(j=1;j<=col;j++)

data[i][j]=d[i][j];

}

void matrix::getdata(double** d){

int i,j;

for(i=1;i<=row;i++)

for(j=1;j<=col;j++)

d[i][j]=data[i][j];

}

matrix matrix::swapcols(int p,int q){

int i,j;

double tmp;

matrix b=matrix(row,col);

for(i=0;i<=row;i++)

for(j=0;j<=col;j++)

b.data[i][j]=data[i][j];

for(i=0;i<=row;i++){

tmp=b.data[i][p];

b.data[i][p]=b.data[i][q];

b.data[i][q]=tmp;

}

return b;

}



matrix matrix::swaprows(int p,int q){

int i,j;

double tmp;

matrix b=matrix(row,col);

for(i=0;i<=row;i++)

for(j=0;j<=col;j++)

b.data[i][j]=data[i][j];

for(i=0;i<=row;i++){

tmp=b.data[p][i];

b.data[p][i]=b.data[q][i];

b.data[q][i]=tmp;

}

return b;

}



matrix matrix::concat(matrix b){

int i,j;

matrix c=matrix(row,col+b.col);

for(i=0;i<=c.row;i++)

for(j=0;j<=col;j++)

c.data[i][j]=data[i][j];

for(i=0;i<=c.row;i++)

for(j=col+1;j<=c.col;j++)

c.data[i][j]=b.data[i][j];

return c;

}



matrix matrix::copy(){

int i,j;

matrix b=matrix(row,col);

for(i=0;i<=row;i++)

for(j=0;j<=col;j++)

b.data[i][j]=data[i][j];

return b;

}



matrix matrix::transpose(){

int i,j;

matrix b=matrix(col,row);

b.row=col;

b.col=row;

for(j=0;j<=col;j++)

for(i=0;i<=row;i++)

b.data[i][j]=data[j][i];

return b;

}



matrix matrix::add(matrix b){

int i,j;

matrix c=matrix(row,col);

c.row=row;

c.col=col;

for(i=0;i<=row;i++)

for(j=0;j<=col;j++)

c.data[i][j]=data[i][j]+b.data[i][j];

return c;

}



matrix matrix::multiply(double x){

int i,j;

matrix c=matrix(row,col);

c.row=row;

c.col=col;

for(i=0;i<=row;i++)

for(j=0;j<=col;j++)

c.data[i][j]=data[i][j]*x;

return c;

}



matrix matrix::multiply(matrix b){

int i,j,k;

double s;

matrix c=matrix(row,b.col);

c.row=row;

c.col=b.col;

for(i=0;i<=c.row;i++){

for(k=0;k<=c.col;k++){

s=0.0;

for(j=0;j<=col;j++)

s+=data[i][j]*b.data[j][k];

c.data[i][k]=s;

}

}

return c;

}





matrix matrix::inverse(double eps,int& err){

int i,j,k,l,blad;

double maxA,d,e;

matrix b=matrix(row,row);

b=copy();

int* M;

M=new int[row+1];

blad=0;

for(i=1;i<=row;i++){

maxA=0.0;

for(j=i;j<=row;j++){

d=b.data[j][i];

if(fabs(maxA)<fabs(d)){

maxA=d;

k=j;

}

}

if(fabs(maxA)<eps) blad=1;

if(!blad){

M[i]=k;

b.data[k][i]=1.0;

for(j=1;j<=row;j++){

d=(double)(b.data[k][j]/maxA);

b.data[k][j]=b.data[i][j];

b.data[i][j]=d;

}

for(j=1;j<=row;j++)

if(j!=i){

d=b.data[j][i];

b.data[j][i]=0.0;

for(l=1;l<=row;l++){

e=d*b.data[i][l];

b.data[j][l]-=e;

}

}

}

}

if(!blad){

for(i=row;i>=1;i--){

k=M[i];

if(k!=i)

for(j=1;j<=row;j++){

d=b.data[j][i];

b.data[j][i]=b.data[j][k];

b.data[j][k]=d;

}

}

}

delete M;

err=blad;

return b;

}

double matrix::det(double eps){

int i,j,k,m;

double s,t;

m=0;

for(i=1;i<=row;i++){

t=fabs(data[i][i]);

k=i;

for(j=i+1;j<=row;j++)

if(fabs(data[j][i])>t) {

t=fabs(data[j][i]);

k=j;

}

if(t<eps){

return 0.0;

}

if(i==k)

t=data[i][i];

else{

m=(m+1)%2;

for(j=row;j>=i;j--){

t=data[k][j];

data[k][j]=data[i][j];

data[i][j]=t;

}

}

t=(double)(1/t);

for(j=i+1;j<=row;j++){

s=-data[j][i]*t;

for(k=i+1;k<=row;k++)

data[j][k]+=data[i][k]*s;

}

}

s=1.0;

for(i=1;i<=row;i++)

s*=data[i][i];

if(m==0) return s;

else return -s;

}

int matrix::rank(double e){

int i,j,k,r;

int imax,jmax;

bool pivot;



for(j=1;j<=col+1;j++) data[0][j]=(double)j;



i=1;



do{



imax=i;

for(k=i+1;k<=row;k++) if (fabs(data[k][i])>fabs(data[imax][i])) imax=k;

pivot=true;

if(fabs(data[imax][i])>e)

  for(j=1;j<=col+1;j++) swap(data[i][j],data[imax][j]);

  else{

  jmax=i;

  for(k=i+1;k<=col;k++) if(fabs(data[i][k])>fabs(data[i][jmax])) jmax=k;



  if(fabs(data[i][jmax])>e)

      for(k=0;k<=row;k++) swap(data[k][i],data[k][jmax]);

  else pivot=false;

  }

  if(pivot){

  for(j=i+1;j<=col+1;j++) data[i][j]/=(double)data[i][i];

   data[i][i]=1.0;

  for(k=i+1;k<=row;k++){

   for(j=i+1;j<=col+1;j++) data[k][j]-=data[i][j]*data[k][i];

    data[k][i]=0.0;

   }

   i++;

   }

   }

   while(!(!pivot||i>row||i>col));

   r=i-1;

   return r;

   }