矩阵基本操作（加减乘、求逆、转置）

看模板，寻找的最好理解，最好用的矩阵基本操作的模板

#define MAXN 100

#define zero(x) (fabs(x)<1e-10)

struct mat

{

int n,m;

double data[MAXN][MAXN];

};

///矩阵加减乘

int add(mat& c,const mat& a,const mat& b)

{

int i,j,k;

if (a.m!=b.m||a.n!=b.n)

return 0;

c.n=a.n;

c.m=a.m;

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

for (j=0; j<c.m; j++)

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

return 1;

}

int jian(mat& c,const mat& a,const mat& b)

{

int i,j,k;

if (a.m!=b.m||a.n!=b.n)

return 0;

c.n=a.n,c.m=a.m;

memset(c.data,0,sizeof(c.data));

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

for (j=0; j<c.m; j++)

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

return 1;

}

int mul(mat& c,const mat& a,const mat& b)

{

int i,j,k;

if (a.m!=b.n)

return 0;

c.n=a.n,c.m=b.m;

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

for (j=0; j<c.m; j++)

for (c.data[i][j]=k=0; k<a.m; k++)

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

return 1;

}

///矩阵的逆

int inv(mat& a)

{

int i,j,k,is[MAXN],js[MAXN];

double t;

if (a.n!=a.m)

return 0;

for (k=0; k<a.n; k++)

{

for (t=0,i=k; i<a.n; i++)

for (j=k; j<a.n; j++)

if (fabs(a.data[i][j])>t)

t=fabs(a.data[is[k]=i][js[k]=j]);

if (zero(t))

return 0;

if (is[k]!=k)

for (j=0; j<a.n; j++)

swap(a.data[k][j],a.data[is[k]][j]);

if (js[k]!=k)

for (i=0; i<a.n; i++)

swap(a.data[i][k],a.data[i][js[k]]);

a.data[k][k]=1/a.data[k][k];

for (j=0; j<a.n; j++)

if (j!=k)

a.data[k][j]*=a.data[k][k];

for (i=0; i<a.n; i++)

if (i!=k)

for (j=0; j<a.n; j++)

if (j!=k)

a.data[i][j]-=a.data[i][k]*a.data[k][j];

for (i=0; i<a.n; i++)

if (i!=k)

a.data[i][k]*=-a.data[k][k];

}

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

{

for (j=0; j<a.n; j++)

if (js[k]!=k)

swap(a.data[k][j],a.data[js[k]][j]);

for (i=0; i<a.n; i++)

if (is[k]!=k)

swap(a.data[i][k],a.data[i][is[k]]);

}

return 1;

}

///转置矩阵，eg2*2的矩阵：[1 2 3 4]->[-2,1,1.5,-0.5]

double det(const mat& a)

{

int i,j,k,sign=0;

double b[MAXN][MAXN],ret=1,t;

if (a.n!=a.m)

return 0;

for (i=0; i<a.n; i++)

for (j=0; j<a.m; j++)

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

for (i=0; i<a.n; i++)

{

if (zero(b[i][i]))

{

for (j=i+1; j<a.n; j++)

if (!zero(b[j][i]))

break;

if (j==a.n)

return 0;

for (k=i; k<a.n; k++)

swap(b[i][k],b[j][k]);

sign++;

}

ret*=b[i][i];

for (k=i+1; k<a.n; k++)

b[i][k]/=b[i][i];

for (j=i+1; j<a.n; j++)

for (k=i+1; k<a.n; k++)

b[j][k]-=b[j][i]*b[i][k];

}

if (sign&1)

ret=-ret;

return ret;

}

int main()

{

mat a,b,c;

while(scanf("%d%d",&a.n,&a.m)!=EOF)

{

for(int i=0; i<a.n; i++)

for(int j=0; j<a.m; j++)

scanf("%lf",&a.data[i][j]);

x=inv(a);

b.n=2,b.m=1;

b.data[0][0]=2;

b.data[1][0]=3;

y=mul(c,a,b);

for(int i=0; i<c.n; i++)

{

for(int j=0; j<c.m; j++)

printf("%lf ",c.data[i][j]);

printf("\n");

}

double m=det(a);

// cout<<m<<endl;

for(int i=0; i<a.n; i++)

{

for(int j=0; j<a.m; j++)

printf("%lf ",a.data[i][j]);

printf("\n");

}

}

return 0;

}