//数据类型定义
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;
}