Hello everyone.
I wrote a simple diagonal preconditioned conjugate gradient solver and I followed the steps provided below.

 ``123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140`` ``````#include #include #include #include #include #include #include #include using namespace std; const double SMALL = 1.0E-30; // used to stop divide-by-zero const double NEARZERO = 1.0E-10; // interpretation of "zero" for printing purposes using vec = vector; // vector using matrix = vector; // matrix (=collection of (row) vectors) //====================================================================== void print( const string &title, const matrix &A ) { if ( title != "" ) cout << title << '\n'; for ( auto &row : A ) { for ( auto x : row ) cout << setw( 15 ) << ( abs( x ) < NEARZERO ? 0.0 : x ); cout << '\n'; } } //====================================================================== void print( const string &title, const vec &A ) { if ( title != "" ) cout << title << '\n'; for ( auto x : A ) cout << setw( 15 ) << ( abs( x ) < NEARZERO ? 0.0 : x ); cout << '\n'; } //====================================================================== matrix matmul( const matrix &A, const matrix &B ) // Matrix times matrix { int rowsA = A.size(), colsA = A[0].size(); int rowsB = B.size(), colsB = B[0].size(); assert( colsA == rowsB ); matrix C( rowsA, vec( colsB, 0.0 ) ); for ( int i = 0; i < rowsA; i++ ) { for ( int j = 0; j < colsB; j++ ) { // scalar product of ith row of A and jth column of B for ( int k = 0; k < colsA; k++ ) C[i][j] += A[i][k] * B[k][j]; } } return C; } //====================================================================== vec matmul( const matrix &A, const vec &V ) // Matrix times vector { int rowsA = A.size(), colsA = A[0].size(); int rowsV = V.size(); assert( colsA == rowsV ); vec C( rowsA, 0.0 ); for ( int i = 0; i < rowsA; i++ ) { for ( int k = 0; k < colsA; k++ ) C[i] += A[i][k] * V[k]; } return C; } //====================================================================== double innerProduct( const vec &U, const vec &V ) // Inner product of U and V { return inner_product( U.begin(), U.end(), V.begin(), 0.0 ); } //====================================================================== double norm( const vec &V ) // L2 norm of V { return sqrt( innerProduct( V, V ) ); } //====================================================================== bool conjugateGradient( const matrix &A, const vec &B, vec &X ) { double TOLERANCE = 1.0e-10; int n = A.size(); X = vec( n, 0.0 ); vec R = B, P = B; for ( int k = 0; k < n; k++ ) { vec Rold = R; // Store previous residual vec AP = matmul( A, P ); double alpha = innerProduct( R, R ) / max( innerProduct( P, AP ), SMALL ); for ( int i = 0; i < n; i++ ) { X[i] += alpha * P[i]; // Next estimate of solution R[i] -= alpha * AP[i]; // Residual } if ( norm( R ) < TOLERANCE ) break; // Convergence test double beta = innerProduct( R, R ) / max( innerProduct( Rold, Rold ), SMALL ); for ( int i = 0; i < n; i++ ) P[i] = R[i] + beta * P[i]; // Next gradient } return true; } //====================================================================== int main() { matrix A = { { 7, 3, 1 }, { 3, 7, 0 }, { 1, 0, 10 } }; vec B = { 1, -5, 2 }; vec X; if ( conjugateGradient( A, B, X ) ) { print( "X:" , X ); cout << "\n\nCheck:\n"; print( "B:" , B ); print( "AX:", matmul( A, X ) ); } else { cout << "Unable to solve\n"; } }``````
Last edited on
Change line 31
`double err=0.0;`
to
`double err=10.0;`

Otherwise, you will never loop, since the test condition isn't met with your initial value of err.
Last edited on