采用部分主元法求方阵A逆矩阵的C代码

//方法优点：直接、易懂、稳定

//方法缺点：存储量大，速度较慢

///////////////////////////////////////////////////////////////////////////////////////////

#include
#include

//打印矩阵
void PrintOut(double M[4][4])
{
 int i,j;
 for(i = 0; i<4;  i++) //行
 {
  printf("\n");
  for(j=0; j<4; j++)
  {
   printf("%4.3f\t",M[i][j]);
  }
 }
 printf("\n");
}

//采用部分主元法的高斯消去法求方阵A的逆矩阵B
//
//A 方阵  (IN)
//n 方阵的阶数 (IN)
//B 方阵  (OUT)
//返回true 说明正确计算出逆矩阵
//返回false说明矩阵A不存在逆矩阵

bool gaussj(double A[4][4], double B[4][4])
{
 int i,j,k;
 double lMax,temp;

 //临时矩阵存放A
 double tt[4][4];
 for(i=0;i<4;i++)
 {
  for(j=0;j<4;j++)
  {
   tt[i][j] = A[i][j];
  }
 }

 //初始化B为单位阵
 for(i=0;i<4;i++)
 {
  for(j=0;j<4;j++)
  {
   if(i!=j)B[i][j] = 0;
   else B[i][j] = 1;
  }
 }
 for(i=0;i<4;i++)
 {
  //寻找主元
  lMax = tt[i][i];
  k = i;
  for(j=i+1;j<4;j++) //扫描从i+1到n的各行
  {
   if( fabs(tt[j][i]) > fabs(lMax))
   {
    lMax = tt[j][i];
    k = j;
   }
  }
  //如果主元所在行不是第i行，进行行交换
  if(k!=i)
  {
   for(j=0;j<4;j++)
   {
    temp = tt[i][j] ;
    tt[i][j] = tt[k][j];
    tt[k][j] = temp;
    //B伴随计算
    temp = B[i][j];
    B[i][j] = B[k][j];
    B[k][j] = temp;
   }
  }
  //判断主元是否是0，如果是，则矩阵A不是满秩矩阵，不存在逆矩阵
  if(tt[i][i] == 0) return false;
  //消去A的第i列除去i行以外的各行元素
  temp = tt[i][i];
  for(j=0;j<4;j++)
  {
   tt[i][j] = tt[i][j] / temp; //主对角线上元素变成1
   B[i][j] = B[i][j] / temp; //伴随计算
  }

  for(j=0;j<4;j++) //行 0 -> n
  {
   if(j!=i)  //不是第i行
   {
    temp = tt[j][i];
    for(k=0;k<4;k++) // j行元素 - i行元素* j行i列元素
    {
     tt[j][k] = tt[j][k] - tt[i][k] * temp;
     B[j][k] = B[j][k] - B[i][k] * temp;
    }
   }
  }
 }
 return true;
}

//计算两个矩阵的乘积
// AB -> C
void Multi(double A[4][4], double B[4][4], double C[4][4])
{
 int i,j,k;
 for(i=0;i<4;i++)
 {
  for(j=0;j<4;j++)
  {
   C[i][j] = 0;
   for(k=0;k<4;k++)
   {
    C[i][j] += A[i][k] * B[k][j];
   }
  }
 }
}

//测试
void  main()
{
 double A[4][4] = { {1 ,3, -5,7},
      {0 ,1, 2,-3},
      {0 ,0, 1, 2},
      {0, 0, 0, 1}
     };
 double B[4][4];
 double C[4][4];

 printf("矩阵:");
 PrintOut(A);

 if(!gaussj(A,B))
 {
  printf("没有逆矩阵！\n");
 }
 else
 {
  printf("\n的逆矩阵是:");
  PrintOut(B);

  printf("\nA*B = ");
  Multi(A,B,C);
  PrintOut(C);
 }
 printf("\nTest Completed!\n");
}

//End

//iwaswzq 2004/11/22

--------------------next---------------------