xiaoyuan.txt

列选主元高斯消去法求解线性方程组,列选主元高斯消去法求解线性方程组

Matrix Matrix::GaussSolve(const Matrix b)
{
 if (r != b.GetRow())  //判断是否方阵
 {
  AfxMessageBox("Gauss Solve Error: 非方阵!");
  exit(1);
 }

 Matrix mtxA(*this); //备份，不改变原矩阵的值。如果原来系数矩阵不需要保留，则可不用设置
 Matrix mtxB(b);

 long row;
 double dMax=0.0;

//消元

 for (long k=0; k<r-1; k++)
 {
  dMax=0.0;
  long row1;
  for (row = k; row < r; row++)    //选主元
   if(fabs(mtxA(row,k)) > dMax)
   {
    dMax = fabs(mtxA(row,k));
    row1=row;
   }
  if (dMax < 1e-8)   //若主元为0，则其右下角子阵全为零，奇异
  {
   AfxMessageBox("Gauss Solve Error:矩阵奇异!");
   exit(1);
  }
  if (row1 !=k)    //交换行，交换一次，行列式变号
  {
   mtxA.SwapRow(k, row1);
   mtxB.SwapRow(k, row1);
  }

  for (row=k+1; row<r; row++)  //计算第k次变换的计算系数
  {
   double temp;
   temp = mtxA(row,k)/mtxA(k,k);

   for (long col=k; col<c; col++)
   {
    double t1,t2,t3;
    t1=mtxA(row,col);
    t2=mtxA(k,col);
    t3=t1-temp*t2;
    mtxA.SetElement(row,col, mtxA(row,col)-temp* mtxA(k,col));
   }
   mtxB.SetElement(row,0,mtxB(row,0)-temp*mtxB(k,0));
  }
 }
 

///回代
 Matrix result(mtxB);
 
 for(row=r-1;row>=0;row--)
 {
  for(k=r-1;k>row;k--)
  {
   double temp;
   temp=result(row,0)-mtxA(row,k)*result(k,0);
   result.SetElement(row,0,temp);
  }
  result.SetElement(row,0,result(row,0)/mtxA(row,row));
 }
 return result;
}