(((silence)))

bool gaussianBackSubstitution(NA::MatDP& m, NA::MatDP& b)
    // 有回代过程的高斯消元法
    // 参数：
    //   m: N x N 的系数矩阵
    //   b: A x N 的右端向量矩阵
    // 返回值：
    //   是否有解
    // 返回时 m 和 b 都已被破坏
    // m 是消去之后的矩阵
    // b 是解向量
{
    if (m.rows()!=m.cols()) {
        return false;
    }
    if (m.rows()!=b.rows()) {
        return false;
    }
    const int n_row=m.rows();
    const double eps=1.0e-10;
    double row_bigest=0, tmp;
    int   i,j,k,nth;

    // 向下消去，过程与高斯-若当方法相同
    for( i=0; i<n_row-1; ++i) {
        row_bigest=0;
        // 选取主元 {{{
        for( j=i; j<n_row; ++j) {
            if(row_bigest<(tmp=fabs(m[j][i]))) {
                row_bigest=tmp;
                nth=j;
            }
        }
        if(fabs(row_bigest)<eps)
            return false;
        if(nth>i){
            m.swapRow(i,nth);
            b.swapRow(i,nth);
        }
        // }}}
        // 消去
        for( j=i+1; j<n_row; ++j) {
            tmp=m[j][i]/m[i][i];
            m[j][i]=0.0;
            for( k=i+1; k<n_row; ++k )
                m[j][k]-=m[i][k]*tmp;
            for( k=0; k<b.cols(); ++k )
                b[j][k]-=b[i][k]*tmp;
        }
    }

    // 回代过程
    for( i=0; i<b.cols(); ++i ) { // 一次可以解出一组方程
        for( j=b.rows()-1; j>=0; --j ) {
            for( k=j+1; k<b.rows(); ++k )
                b[j][i]-=b[k][i]*m[j][k]; // 将已算出的结果(b[k][i])带入
            b[j][i]/=m[j][j];
        }
    }
    return true;
}

// 顺便，还有高斯-若当消去法：

bool gaussianJordan(NA::MatDP& m, NA::MatDP& b)
    // 高斯-若当 消去法
    // 参数：
    //   m: N x N 的系数矩阵
    //   b: A x N 的右端向量矩阵
    // 返回值：
    //   是否有解
    // 返回时 m 和 b 都已被破坏
    // m 是消去之后的矩阵
    // b 是解向量
{
    if (m.rows()!=m.cols()) {
        return false;
    }
    if (m.rows()!=b.rows()) {
        return false;
    }
    const int n_row=m.rows();
    const double eps=1.0e-10;
    double row_bigest=0, tmp;
    int   i,j,k,nth;

    // 向下消去
    for( i=0; i<n_row-1; ++i) {
        row_bigest=0;
        // 选取主元 {{{
        for( j=i; j<n_row; ++j) {
            if(row_bigest<(tmp=fabs(m[j][i]))) {
                row_bigest=tmp;
                nth=j;
            }
        }
        if(fabs(row_bigest)<eps)
            return false;
        if(nth>i){
            m.swapRow(i,nth);
            b.swapRow(i,nth);
        }
        // }}}
        // 消去
        for( j=i+1; j<n_row; ++j) {
            tmp=m[j][i]/m[i][i];
            m[j][i]=0.0;
            for( k=i+1; k<n_row; ++k )
                m[j][k]-=m[i][k]*tmp;
            for( k=0; k<b.cols(); ++k )
                b[j][k]-=b[i][k]*tmp;
        }
    }

    // 向上消去
    for( i=n_row-1; i>=0; --i ) {
        tmp=m[i][i];
        if(fabs(tmp)<eps)
            return false;
        m[i][i]=1.0;
        for( j=i+1; j<n_row; ++j ) {
            m[i][j] /= tmp;
        }
        for( j=0; j<b.cols(); ++j ) {
            b[i][j] /= tmp;
        }
        for( j=i-1; j>=0; --j ) {
            tmp=m[j][i];
            m[j][i]=0.0;
            for( k=i+1; k<n_row; ++k ) {
                m[j][k]-=m[i][k]*tmp;
                if(fabs(m[j][k])<eps)
                    m[j][k]=0.0;
            }
            for( k=0; k<b.cols(); ++k ) {
                b[j][k]-=b[i][k]*tmp;
                if(fabs(b[j][k])<eps)
                    b[j][k]=0.0;
            }
        }
    }

    return true;
}