藍影

關於 Gauss Jordan (不知道是不是這麼拼) 消去法
基本上大家都知道怎麼做，
這裡的 code 分二段，一段是用高斯消去法求反矩陣(MINVERSE)，
另一段是用高期消去法求連立方程式之解
以下的 code 便是實現其作法
若對 code 有問題，歡迎回覆進行討論

// =======================================
// author  : edison.shih.
// filename: matrix.cpp
// date    : 2010/2/7
// mail    : forget72@hotmail.com
// ** all copyright reverve **
// =======================================
// ==============================
// 計算反矩陣
/*
use gauss jodarn method.
return value:
-1: have a row all 0
0 : exist inverse
1 : have a col all 0
*/
int MINVERSE( double **a,
             unsigned dim,
             double **_result)
{
    unsigned i, j;
    double x=0.0;
    SetUnitMatrix(_result, dim);
   
    double **tmp = (double**)malloc(sizeof(double*)*dim);
    for(i=0; i<dim; i++) tmp[i] = (double*)malloc(sizeof(double)*dim);
   
    CopyMatrix(tmp, a, dim, dim);
   
    for(i=0; i<dim; i++){
       
        // 該列全為0, return -1
        if(IsRi0(tmp, i, dim)) return -1;
       
        if(tmp[i][i]==0.0){
            // 交換至該元素為非0
            bool flag = false;
            for(j=i+1; j<dim; j++){
                Rij(tmp, i, j, dim);
                Rij(_result, i, j, dim);
                if(tmp[i][i]!=0.0) {
                    break;
                    flag = true;
                }
            }
            // 交換不到非0 元素，return 1
            if(flag==false) return 1;
        }
       
        // 以下精華
        x = 1.0 / tmp[i][i];
        Rix(tmp, i, x, dim);
        Rix(_result, i, x, dim);
       
        for(j=0; j<dim; j++){
            if(i==j) continue;
            x = -tmp[j][i] / tmp[i][i];
            Rijx(tmp, i, j, x, dim);
            Rijx(_result, i, j, x, dim);
        }
    }
   
    for(i=0; i<dim; i++) free(tmp[i]);
    free(tmp);
    return 0;
}

// ==============================
// 使用 gauss jordan 求連立方程式
/*
return  value:
-1: 無解
0 : 有一解
1 : 無限多解
*/
int GAUSSJORDAN(double **a,
                double *ans,
                unsigned dim)
{
   
    unsigned i, j;
    double x=0.0;
   
    double **tmp = (double**)malloc(sizeof(double*)*dim);
    for(i=0; i<dim; i++) tmp[i] = (double*)malloc(sizeof(double)*(dim+1));
   
    CopyMatrix(tmp, a, dim, dim+1);
   
    for(i=0; i<dim; i++){
       
        // 該列全為0, 無限多解, return 1
        if(IsRi0(tmp, i, dim+1)) return -1;
       
        if(tmp[i][i]==0.0){
            // 交換至該元素為非0
            bool flag = false;
            for(j=i+1; j<dim; j++){
                Rij(tmp, i, j, dim+1);
                if(tmp[i][i]!=0.0) {
                    break;
                    flag = true;
                }
            }
            // 交換不到非0 元素，無解, return -1
            if(flag==false) return -1;
        }
       
        // 以下精華
        x = 1.0 / tmp[i][i];
        Rix(tmp, i, x, dim+1);
       
        for(j=0; j<dim; j++){
            if(i==j) continue;
            x = -tmp[j][i] / tmp[i][i];
            Rijx(tmp, i, j, x, dim+1);
        }
    }
    for(i=0; i<dim; i++) {
        ans[i] = tmp[i][dim];
    }
   
    for(i=0; i<dim; i++) free(tmp[i]);
    free(tmp);
    return 0;
}