#ifndef _MATRIX_H_
/********************************
* matrix.h *
* 矩阵运算库 *
* edit 2004-8-9 *
* 版本1.0 *
*********************************
矩阵计算的几个函数
atV计算转置
abV计算a×b
abaddV计算a+b
absubV计算a-b
inputaV输入矩阵
printaV输出矩阵
求逆
int a_1V(double *in,double *out,int n)
转置
int atV(double *in,double *out,int m,int n)
矩阵相乘
int abV(double *a,double *b,double *out,int am,int an,int bm,int bn)
矩阵相加
int abaddV(double *a,double *b,double *out,int m,int n)
矩阵相减
int absubV(double *a,double *b,double *out,int m,int n)
输出
printaV(double *a,int m,int n)
输入
inputaV(double *a,int m,int n)
*/
/* 转置一个矩阵,矩阵存在一个一维数组in中,输出也存在一个一维数组out
矩阵列数和维数由n指定,
调用可以使用数组,例
二维数组调用方式
atV(double in[n][n],double out[n][m],int m,int n)
一维调用方式
atV(double in[k],double out[k],int m,int n)
指针调用方式
atV(double *in,double *out,int m,int n)
*/
#include <malloc.h>
#include <math.h>
int atV(double *in,double *out,int m,int n)
{
int i,j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
*(out+ j*m +i)=*(in+ i*n +j);
}
return 1;
}
/* 矩阵相乘
out
am an 为矩阵a的行数和列数
bm nn 为矩阵b的行数和列数
矩阵列数和维数由n指定
调用可以使用数组,例
二维数组调用方式
avV(double a[am][an],double b[bm][bn],double out[am][bn],int am,int an,int bm,
int bn)
一维调用方式
avV(double a[k],double b[l],double out[j],int am,int an,int bm,int bn)
指针调用方式
avV(double *a,double *b,double *out,int am,int an,int bm,int bn)
*/
int abV(double *matrixa,double *matrixb,double *out,int am,int an,int bm,int bn)
{
int i,j,k;
double temp;
if(an!=bm) {return 0;} /*矩阵不可以计算*/
for(i=0;i<am;i++)
{
for(j=0;j<bn;j++)
{
for(k=0,temp=0;k<an;k++)
temp=temp + *(matrixa+i*an+k) * *(matrixb+k*bn+j);
*(out+i*bn+j) = temp;
}
}
return 1;
}
/*输入一个m*n的矩阵,也可以输出m*n数组以及一维数组按m*n输出*/
int inputaV(double *matrixa,int m,int n)
{
int i,j;
for(i=0;i<m;i++)
{
printf("请输入第%d行/n",i+1);
for(j=0;j<n;j++)
{
scanf("%lf",(matrixa+i*n+j));
}
}
return 1;
}
/*输入一个m*n矩阵,也可以输出m*n数组以及一维数组按m*n输出*/
printaV(double *matrixa,int m,int n)
{
int i,j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
printf("%lf ",*(matrixa+i*n+j));
printf("/n");
}
printf("/n");
}
void initiones(double *q,int n) /*形成单位矩阵*/
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{ if(i==j)
q[i*n+j]=1;
else
q[i*n+j]=0;
}
}
void swap(double *r,int i,int line,int n) /*换行*/
{
int j;
double m;
for(j=0;j<n;j++)
{ m=r[i*n+j];
r[i*n+j]=r[line*n+j];
r[line*n+j]=m;
}
}
int solve(double *p,double *q,int n) /*形成下三角矩阵*/
{
int i ,j,k,m,line;
double max,temp,savep;
for (i=0;i<n;i++)
{ max=fabs(p[i*n+i]);
temp=p[i*n+i];
line=i;
for (j=i+1;j<n;j++)
{ if( fabs(p[j*n+i]) > max ) /*找出绝对值最大元*/
{line=j;
max=fabs(p[j*n+i]);
temp=p[j*n+i]; }
}
if(max<=1e-5) /*三角矩阵的对角元等于0则无解*/
{
return 0;}
if(line!=i)
{ swap(p,i,line,n);
swap(q,i,line,n);}
for (k=0;k<n;k++)
{p[i*n+k]/=temp;
q[i*n+k]/=temp;}
for(k=i+1;k<n;k++) /*消元,成下三角阵;*/
{ savep=p[k*n+i];
for( m=0;m<n;m++)
{ q[k*n+m]=q[k*n+m]-q[i*n+m]*savep;
p[k*n+m]=p[k*n+m]-p[i*n+m]*savep;}
}
}
return 1;
}
void backloop(double *p,double *q,int n) /*回带,形成单位矩阵*/
{
int i,j,k;
double savep;
for(i=n-1;i>0;i--)
{
for(j=i-1;j>=0;j--)
{
savep=p[j*n+i];
p[j*n+i]=p[j*n+i]-p[i*n+i]*savep;
for(k=0;k<n;k++)
q[j*n+k]=q[j*n+k]-q[i*n+k]*savep;
}
}
}
/*求a的逆阵,in为输入阵,out为输出的逆阵,n为维数*/
int a_1V(double *in,double *out,int n)
{
double *p;
int i,j;
p=(double *)malloc(sizeof(double)*(n*n));
for(i=0;i<n;i++)
for(j=0;j<n;j++) p[j*n+i]=in[j*n+i];
initiones(out,n);
if(!solve(p,out,n)) {return 0;}
backloop(p,out,n);
return 1;
}
/*求矩阵和,a,b为输入阵,out为和的矩阵,m,n为行和列*/
int abaddV(double *matrixa,double *matrixb,double *out,int m,int n)
{
int i,j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
*(out+i*n+j)=*(matrixa+i*n+j)+*(matrixb+i*n+j);
}
return 1;
}
/*求矩阵差,a,b为输入阵,out为差的矩阵,m,n为行和列*/
int absubV(double *matrixa,double *matrixb,double *out,int m,int n)
{
int i,j;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
*(out+i*n+j)=*(matrixa+i*n+j)-*(matrixb+i*n+j);
}
return 1;
}
#endif