高斯消去法求反矩陣 in C++

在過去解二元(兩個未知數)一次方程式時，
一般來說，假設兩式有解的情況下，
會很直觀地先將某一個未知數(假設為X)消除，
為了要消除X，進而使用兩個式子的X的最小公倍數係數，
將兩個式子的X係數變成相同，就可以將兩條式子相減，
後續可以得到Y係數的解，最後再將Y係數的解帶回某一式子，
即可得到X,Y的解。

以上是兩個未知數的情況下求解常用的方法，
可是在多個未知數(N-Dimension)呢？

德國的數學家－高斯(Johann Karl Friedrich Gauß)，
發明了一套方法，一次可以解多個未知數的方法，
也就是本文要介紹的高斯消去法(Gauss Elimination)。

其演算過程如下：

將S的Eq(1)，同除以a1後得新Eq(1’)，將Eq(1’)乘以-a2並加入Eq(2)，再將Eq(1’)乘以-a3加入Eq(3)，會將原本的a2以及a3的位置消為0，再重複同樣的過程，就可以將S的形式轉變成T，而e1&e2&e3就是解。

將方程式轉成矩陣形式如下：

B矩陣就是解。

C++的程式碼如下：

#include 
#include 
using namespace std ;


#define MATRIX_MAXSIZE 10
double MATRIX[MATRIX_MAXSIZE][MATRIX_MAXSIZE] ;
double INVMATRIX[MATRIX_MAXSIZE][2 * MATRIX_MAXSIZE] ;

int MATRIX_SIZE = 3 ;

void PRINT_MATRIX(double *_data , int M , int N , int matrixY) {

 for ( int i = 0 ; i < M ; ++i ) {
  double *ptr = _data + i * matrixY;
  for ( int j = 0 ; j < N ; ++j ) {
   cout << setprecision(4) << ptr[j] << "\t";
  }
  if ( matrixY == 2 * MATRIX_MAXSIZE ) {
   for ( int j = 0 ; j < N ; ++j ) {
    cout << setprecision(4) << ptr[j + MATRIX_MAXSIZE] << "\t";
   }
  }
  cout << endl ;
 }
} ;

void PRINT_ORIMATRIX() {
 PRINT_MATRIX(MATRIX[0] , MATRIX_SIZE , MATRIX_SIZE , MATRIX_MAXSIZE) ;
} ;

void PRINT_INVMATRIX() {
 PRINT_MATRIX(INVMATRIX[0] , MATRIX_SIZE , MATRIX_SIZE , MATRIX_MAXSIZE * 2) ;
} ;

void INIT_MATRIX() {
 for ( int i = 0 ; i < MATRIX_MAXSIZE ; ++i ) {
  for ( int j = 0 ; j < MATRIX_MAXSIZE ; ++j ) {
   MATRIX[i][j] = 0.0 ;
   INVMATRIX[i][j] = 0.0 ;
   if ( i == j ) {
    INVMATRIX[i][j + MATRIX_MAXSIZE] = 1.0 ;
   }else {
    INVMATRIX[i][j + MATRIX_MAXSIZE] = 0.0 ;
   }
  }
 }
} ;

void RAND_SETMATRIX() {
 std::cin >> MATRIX_SIZE ;
 INIT_MATRIX() ;
 srand(time(NULL)) ;

 for ( int i = 0 ; i < MATRIX_SIZE ; ++i ) {
  for ( int j = 0 ; j < MATRIX_SIZE ; ++j ) {
   int num = rand() % 250 ;
   MATRIX[i][j] = num;
   INVMATRIX[i][j] = num ;
  }
 }

 //test swap function...
 //INVMATRIX[1][0] = MATRIX[1][0] = MATRIX[0][0] ;
 //INVMATRIX[1][1] = MATRIX[1][1] = MATRIX[0][1] ;

} ;

void VERIFY_MATRIX() {
 for ( int i = 0 ; i < MATRIX_SIZE ; ++i ) {
  for ( int j = 0 ; j < MATRIX_SIZE ; ++j ) {
   double dotResult = 0 ;
   for ( int k = 0 ; k < MATRIX_SIZE ; ++k ) {
    dotResult += (MATRIX[i][k] * INVMATRIX[k][j + MATRIX_MAXSIZE]);
   }
   cout << setprecision(4) << dotResult << "\t" ;
  }
  cout << endl ;
 } 
} ;

bool GAUSSIAN_ELIMINATION() {
 bool noAns = false ;
 for ( int i = 0 ; i < MATRIX_SIZE && !noAns; ++i ) {
  int j ;
  if ( INVMATRIX[i][i] == 0 ) {
   //find another pivot
   for ( j = i + 1 ; j < MATRIX_SIZE ; ++j ) {
    if ( INVMATRIX[j][i] != 0 ) {
     //find pivot
     //swap two row
     for ( int k = 0 ; k < 2 * MATRIX_MAXSIZE ; ++k ) {
      double tmp = INVMATRIX[i][k] ;
      INVMATRIX[i][k] = INVMATRIX[j][k] ;
      INVMATRIX[j][k] = tmp ;
     }
     break ;
    }
   }
   if ( j == MATRIX_SIZE ) {
    noAns = true ;
   }
  }
  if ( !noAns ) {
   double factor = 1.0 / INVMATRIX[i][i] ;
   for ( j = 0 ; j < MATRIX_SIZE ; ++j ) {
    INVMATRIX[i][j] *= factor;
    INVMATRIX[i][j + 10] *= factor;
   }
   for ( j = 0 ; j < MATRIX_SIZE ; ++j ) {
    factor = -INVMATRIX[j][i];
    //avoid self-elimination 
    if ( i != j ) {
     for ( int k = 0 ; k < MATRIX_SIZE ; ++k ) {
      //
      INVMATRIX[j][k] += factor * INVMATRIX[i][k];
      INVMATRIX[j][k+MATRIX_MAXSIZE] += factor * INVMATRIX[i][k+MATRIX_MAXSIZE];
     }
    }
   }
  }
 }
 return !noAns;
} ;



void main() {
 RAND_SETMATRIX() ;
 cout << "ORI: \r\n";
 PRINT_ORIMATRIX() ;
 cout << "INV_MATRIX: \r\n";
 PRINT_INVMATRIX() ;

 bool ans = GAUSSIAN_ELIMINATION() ;
 if ( ans ) {
  cout << "AFTER ELIMINATION: \r\n" ;
  PRINT_INVMATRIX() ;
  cout << "VERIFY: \r\n" ;
  VERIFY_MATRIX() ;
 }else {
  cout << "NO ANS. \r\n" ;
 }
 
 system("pause") ;
} ;