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;
}