C语言矩阵求逆(源代码)

#include

double * MatrixOpp(double *A,int m,int n); /*矩阵求逆*/

double * MatrixInver(double *A,int m,int n); /*矩阵转置*/

double Surplus(double A[],int m,int n); /*求矩阵行列式*/

double * MatrixOpp(double A[],int m,int n) /*矩阵求逆*/

{

int i,j,x,y,k;

double *SP=NULL,*AB=NULL,*B=NULL,X,*C;

SP=(double *)malloc(m*n*sizeof(double));

AB=(double *)malloc(m*n*sizeof(double));

B=(double *)malloc(m*n*sizeof(double));

X=Surplus(A,m,n);

X=1/X;

for(i=0;i for(j=0;j {

for(k=0;k B[k]=A[k];

{

for(x=0;x B[i*n+x]=0;

for(y=0;y B[m*y+j]=0;

B[i*n+j]=1;

SP[i*n+j]=Surplus(B,m,n);

AB[i*n+j]=X*SP[i*n+j];

}

}

C=MatrixInver(AB,m,n);

return C;

}

double * MatrixInver(double A[],int m,int n) /*矩阵转置*/

{

int i,j;

double *B=NULL;

B=(double *)malloc(m*n*sizeof(double));

for(i=0;i for(j=0;j B[i*m+j]=A[j*n+i];

return B;

}

double Surplus(double A[],int m,int n) /*求矩阵行列式*/

{

int i,j,k,p,r;

double X,temp=1,temp1=1,s=0,s1=0;

if(n==2)

{

for(i=0;i for(j=0;j if((i+j)%2) temp1*=A[i*n+j];

else temp*=A[i*n+j];

X=temp-temp1;

}

else

{

for(k=0;k {

for(i=0,j=k;i temp*=A[i*n+j];

if(m-i)

{

for(p=m-i,r=m-1;p> 0;p--,r--)

temp*=A[r*n+p-1];

}

s+=temp;

temp=1;

}

for(k=n-1;k> =0;k--)

{

for(i=0,j=k;i =0;i++,j--)

temp1*=A[i*n+j];

if(m-i)

{for(p=m-1,r=i;r temp1*=A[r*n+p];}

s1+=temp1;

temp1=1;

}

X=s-s1;

}

return X;

}

/* Test */

int main()

{

int i,j;

double arr[5][5], *result, *t=arr;

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

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

scanf( "%lf

IZE]; //值

};

CMatrix InvMatrix(CMatrix ); //矩阵求逆

CMatrix operator*(CMatrix& a,CMatrix &b); //矩阵乘法

CVector operator*(CMatrix& a,CVector &b);//矩阵乘以向量

CMatrix operator-(CMatrix &a,CMatrix &b); //矩阵减法

void Eliminate(CMatrix m,CVector &b,CVector &solution); //选主元消元法

void JacobiIter(CMatrix m,CVector b,CVector &solution); //解方程组的Jacibi迭代法

#endif

CMatrix InvMatrix(CMatrix a){ //逆

if(a.Row() <=0 || a.Col() <=0)

throw "矩阵输入错误!无法执行 ";

if(a.Row()!=a.Col()) throw "非方阵,无法求逆! ";

CMatrix rhs(a.Row(),a.Col());

int i,j,k,t;

double ax,temp,s;

int flow;

for(i=0;i for(j=0;j rhs.value[i][j]=0;

rhs.value[i][i]=1;

}

for(k=0;k s=0;

for(j=k;j if(abs(a.value[j][k])> =s){

s=abs(a.value[j][k]);

flow=j;

}

}

for(t=0;t temp=a.value[k][t];

a.value[k][t]=a.value[flow][t];

a.value[flow][t]=temp;

temp=rhs.value[k][t];

rhs.value[k][t]=rhs.value[flow][t];

rhs.value[flow][t]=temp;

}

for(i=k+1;i ax=a.value[i][k]/a.value[k][k];

for(t=0;t a.value[i][t]=a.value[i][t]-ax*a.value[k][t];

rhs.value[i][t]=rhs.value[i][t]-ax*rhs.value[k][t];

}

}

}

//********得到一个上三角阵，下面将之化为对角阵***********

for(i=a.Row()-2;i> =0;i--){

for(j=i+1;j ax=a.value[i][j]/a.value[j][j];

for(t=0;t a.value[i][t]=a.value[i][t]-ax*a.value[j][t];

rhs.value[i][t]=rhs.value[i][t]-ax*rhs.value[j][t];

}

}

}

for(i=0;i for(j=0;j rhs.value[i][j]/=a.value[i][i];

}

for(j=0;j a.value[i][j]/=a.value[i][i];

}

}

return rhs;

}