//============================================================================ // // Program/Module // from // Selfverifying Solvers for Dense Systems of // Linear Equations Realized in C-XSC // // Carlos Holbig and Walter Kraemer, 2003. // // Program developed in C-XSC by Bernardo Frederes Kramer Alcalde, // Paulo Sergio Morandi Junior and Carlos Amaral Holbig // //============================================================================ //--------------------------------------------------------------------------- // File: lss_aprx (implementation) // Purpose: Compute an approximate inverse of the square matrix A // Global functions: // MINV(): Computes the approximate inverse of a square matrix A using the // Gauss-Jordan algorithm. //--------------------------------------------------------------------------- #include using namespace cxsc; using namespace std; const real eps = 1e-15; //--------------------------------------------------------------------------- // MINV computes approximate inverse of the matrix A by use of the // Gauss-Jordan algorithm. // A : input matrix and outpout of approximate inverse // err : Error indicator = 0 everything ok, = 1 matrix singular //--------------------------------------------------------------------------- void MINV( rmatrix& A, int& error ) { int i, j, k, au1, au2, ao1, ao2; bool ok; real h; rvector aux( Lb(A,ROW),Ub(A,ROW) ); intvector v( Lb(A,ROW),Ub(A,ROW) ); int l; au1 = Lb(A,ROW); ao1 = Ub(A,ROW); au2 = Lb(A,COL); ao2 = Ub(A,COL); ok = ( (ao1 - au1) == (ao2 - au2) )? true : false; // square matrix ? // Vector that will control what's happening with the column for(i = au1; i <= ao1; i++) v[i] = i; i = au1-1; // Row index j = au2-1; // Column index while ( (j < ao2) && ok ) { i++; j++; // i and j run simultaniously // pivot search if ( j < ao2) { l = i; for( k = i+1; k <= ao1; k++) if( abs(A[k][j] ) > abs( A[l][j] ) ) l = k; // row interchange if(i != l) { aux = A[l]; A[l] = A[i]; A[i] = aux; k = v[l]; v[l] = v[i]; v[i] = k; } } // transformation if( A[i][j] == 0 ) ok = false; else { h = 1.0 / A[i][j]; Col(A,j) = h * rvector( Col(A,j) ); A[i][j] = h; for( k = au2; k <= ao2; k++) if( k != j) { for( l = au1; l <= ao1; l++) if( l != i ) { h = A[l][k] - A[l][j] * A[i][k]; if( h == 0.0) A[l][k] *= eps; else A[l][k] = h; } A[i][k] *= -A[i][j]; // A[i][k] = -A[i][k] * A[i][j] } } } // column interchange if(ok) for( k = au1; k <= ao1; k++) { l = au2 + k - au1; while( v[k] != k ) { i = v[k]; v[k] = v[i]; v[i] = i; j = au2 + i - au1; aux = rvector( Col(A,j) ); Col(A,j) = Col(A,l); Col(A,l) = aux; } } error = !ok; } // end MINV