Metoda eliminacji Gaussa

/* This version preserves the matrix 'a' and the vector 'b'. */

int cged2( complex < double > ** a, complex < double > * b,
complex < double > * x, int n )
{
    complex < double > tmp, pvt, * t, ** aa, * bb;
    int i, j, k, itmp, retval;
   
    // Initialize return value for successful execution.
    retval = 0;
   
    // Create and initialize working storage
    aa = new complex < double > *[ n ];
    bb = new complex < double >[ n ];
    for( i = 0; i < n; i++ ) {
        aa[ i ] = new complex < double >[ n ];
        for( j = 0; j < n; j++ ) {
            aa[ i ][ j ] = a[ i ][ j ];
        }
        bb[ i ] = b[ i ];
    }
   
    // Main loop
    for( i = 0; i < n; i++ ) { // outer loop on rows
        pvt = aa[ i ][ i ]; // get pivot value
        if( !abs( pvt ) ) {
            for( j = i + 1; j < n; j++ ) {
                if(( pvt = abs( aa[ j ][ i ] ) ) != 0.0 ) break;
               
            }
            if( !abs( pvt ) ) {
                retval = 1;
                goto _100; // pull the plug!
            }
            t = aa[ j ]; // swap matrix rows...
            aa[ j ] = aa[ i ];
            aa[ i ] = t;
            tmp = bb[ j ]; // ...and result vector
            bb[ j ] = bb[ i ];
            bb[ i ] = tmp;
        }
        for( k = i + 1; k < n; k++ ) { // (virtual) eliminate column
            tmp = aa[ k ][ i ] / pvt;
            for( j = i + 1; j < n; j++ ) {
                aa[ k ][ j ] -= tmp * aa[ i ][ j ];
            }
            bb[ k ] -= tmp * bb[ i ];
        }
    }
    // Do back substitution
    for( i = n - 1; i >= 0; i-- ) {
        x[ i ] = bb[ i ];
        for( j = n - 1; j > i; j-- ) {
            x[ i ] -= aa[ i ][ j ] * x[ j ];
        }
        x[ i ] /= aa[ i ][ i ];
    }
    _100:
    for( i = 0; i < n; i++ ) {
        delete[] aa[ i ];
    }
    delete[] aa;
    delete[] bb;
    return retval;
}