solve.cpp

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

#include <iostream>
#include <iomanip>

using namespace std;

// How to declare FORTRAN functions:
// use extern "C" to stop compiler from
// "name mangling", and add _ to end of
// function name as it appears in the
// FORTRAN code
extern "C" void gauss_elim_( double A[][4], double L[][4], double U[][4], int *n );

void PrintMatrix( const char *sName, double A[][4] ) {
    cout << sName << " = " << endl;
    int n = 4;
    for ( int i = 0; i < n; i++ ) {
        for ( int j = 0; j < n; j++ ) 
            cout << setw(10) << A[i][j];
        cout << endl;
    }
}

int main() {

    int n = 4;
    int i,j;

    double A[4][4];
    double L[4][4];
    double U[4][4];

    double b[4];
    double x[4];
    double y[4];

    double h = 1.256637061435917;

    for ( i = 0; i < n; i++ ) {
        for ( j = 0; j < n; j++ )
            A[i][j] = 0.0;
        A[i][i] = -2.0 / (h*h);
        if ( i > 0 ) {
            A[i][i-1] = 1.0 / (h*h);
            A[i-1][i] = 1.0 / (h*h);
        }
    }

    b[0] = 0.951056516295154;
    b[1] = 0.587785252292473;
    b[2] = -b[1];
    b[3] = -b[0];

	// How to call FORTRAN functions: make 
	// sure all arguments are pointers, 
	// because FORTRAN uses pass-by-reference
    gauss_elim_(A,L,U,&n);

    PrintMatrix( "A", A );
    PrintMatrix( "L", L );
    PrintMatrix( "U", U );
    
    for ( i = 0; i < n; i++ ) {
        y[i] = b[i];
        for ( j = 0; j < i; j++ ) 
            y[i] -= L[i][j] * y[j];
    }

    for ( i = n - 1; i >= 0; i-- ) {
        x[i] = y[i];
        for ( j = i + 1; j < n; j++ )
            x[i] -= U[i][j] * x[j];
        x[i] /= U[i][i];
    }

    cout << "x = " << endl;
    for ( i = 0; i < n; i++ )
        cout << setw(10) << x[i] << endl;

    return 0;
}