Function GaussElimination

plugin/seq/plotPDF.cpp:748–796

// End of SimplePDFModule
//----------------------------------------------------------------------


void GaussElimination(double *const x, double *const *const a, const int N)
{
    const double EPS = 1e-10;

    // forward elimination
    for(int k = 0; k < N-1; k++){

	{ // begin pivoting
	    int row = k;
	    
	    for(int r = k+1; r < N; r++)
		if(fabs(a[row][k]) < fabs(a[r][k]))
		    row = r;

	    if(fabs(a[row][k]) < EPS){
		std::cerr << "singular matrix : " << row << std::endl;
		exit(1);
	    }

	    // k == row: no pivoting
	    if(k != row){
		for(int j = 0; j <= N; j++){
		    const double tmp = a[k][j];
		    a[k][j] = a[row][j];
		    a[row][j] = tmp;
		}
	    }
	} // end pivoting

	const double d = 1 / a[k][k];
        for(int i = k+1; i < N; i++){
            for(int j = k+1; j <= N; j++)
                a[i][j] -= a[i][k] * a[k][j] *d;
            a[i][k] = 0;
        }
    }

    // backward substitution
    for(int i = N-1; i >= 0; i--){
	for(int j = i+1; j < N; j++)
	    a[i][N] -= a[i][j]*a[j][N];
	a[i][N] /= a[i][i];
    }

    for(int i = 0; i < N; i++)
	x[i] = (fabs(a[i][N]) < EPS)? 0: a[i][N];

    return;
}