C语言实现矩阵的四则运算

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