看模板,寻找的最好理解,最好用的矩阵基本操作的模板
#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;
}
版权声明:本文为博主原创文章,未经博主允许不得转载。
原文:http://blog.csdn.net/li1004133206/article/details/48060485