gauss_elim.cpp

斯坦福Energy211/CME211课《c++编程——地球科学科学家和工程师》的课件

// If function is to be called by FORTRAN 
// code, add _ to end of function name here, 
// but leave it out in FORTRAN code.  Also,
// precede declaration by extern "C" so that
// name-mangling won't occur and FORTRAN
// linker can find this function
extern "C" void gauss_elim_( double A[][4], double L[][4], double U[][4], int *pn )
{
    int n = *pn;

    for ( int i = 0; i < n; i++ ) {
        for ( int j = 0; j < n; j++ ) {
            L[i][j] = 0.0;
            U[j][i] = A[j][i];
        }
        L[i][i] = 1.0;
    }

	// the matrices that are passed to this
	// function are stored in column order,
	// so we work with their transposes
    for ( int j = 0; j < n - 1; j++ ) {
        for ( int i = j + 1; i < n; i++ ) {
            double m = U[j][i] / U[j][j];
            for ( int k = j; k < n; k++ )
                U[k][i] -= m * U[k][j];
            L[j][i] = m;
        }
    }
}