矩阵运算库

#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