//数据类型定义 typedef struct{ double *tab; int row,col; }Matrix; Matrix newMatrix(int m,int n){ Matrix a; assert(m>=0&&n>=0); a.row=m;a.col=n; assert(a.tab=(double *)calloc(m*n,sizeof(double))); return a; } Matrix newMatrixByArray(double *a,int m,int n){ Matrix b=newMatrix(m,n); memcpy(b.tab,a,n*m*sizeof(double)); return b; } void clrMatrix(Matrix a){ free(a,tab); } //运算实现 Matrix transform(Matrix a){ //矩阵转置 int i,j,m=a.col,n=a.row; Matrix at=newMatrix(m,n); for(i=0;i<n;i++){ for(j=0;j<n;j++){ at.tab[i*n+j]=a.tab[j*m+i]; } } return at; } Matrix matrixAdd(Matrix a,Matrix b){ //矩阵相加 Matrix c; int i,j,m=a.row,n=a.col; assert(b.row==m&&b.col==n); c=newMatrix(m,n); for(i=0;i<m;i++){ for(j=0;j<n;j++){ c.tab[i*n+j]=a.tab } } } Matrix matrixMultiply(Matrix a,Matrix b){ //矩阵相乘 Matrix c; int i,j,k; assert(a.col==b.row); c=newMatrix(a.row,b.col); for(i=0;i<a.row;i++){ for(j=0;j<b.col;j++){ for(k=0;k<a.col;k++){ c.tab[i*c.col+j]+=a.tab[i*a.col+k]*b.tab[k*b.col+j]; } } } return c; } //LUP分解实现矩阵求逆 double *lupSolve(Matrix lu,int *pi,double *b){ //矩阵求逆 int i,j,n=lu.row; double *x=(double *)calloc(n,sizeof(double)); for(i=0;i<n;i++){ x[i]=b[pi[i]]; for(j=0;j<i;j++){ x[i]=x[i]-lu.tab[i*n+j]*x[j]; } } for(i=n-1;i>=0;i--){ for(j=i+1;j<n;j++){ x[i]=x[i]-lu.tab[i*n+j]*x[j]; } x[i]=x[i]/lu.tab[i*n+i]; } return x; } Matrix matrixInverse(Matrix a){ double *c=NULL,*xi=NULL; int n=a.row,*pi=(int *)calloc(n,sizeof(int)),i; Matrix b=newMatrix(n,n),lu,x; lu=lup(a,pi); for(i=0;i<n;i++){ if(e) free(e); e=(double *)calloc(n,sizeof(double)); e[i]=1.0; if(xi) free(xi); xi=luSolve(lu,pi,e); memcpy(b.tab+i*n,xi,n*sizeof(double)); } x=transform(b); clrMatrix(lu);clrMatrix(b); free(pi);free(e);free(xi); return c; }