From 789e892865817151329c3d5804f538b91320583d Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Tue, 18 Jul 2017 17:57:13 +0200 Subject: [PATCH] Save current state --- .../iterative/GeneralisedMinimalResidual.h | 389 ++++++++++++------ 1 file changed, 253 insertions(+), 136 deletions(-) diff --git a/lib/algorithms/iterative/GeneralisedMinimalResidual.h b/lib/algorithms/iterative/GeneralisedMinimalResidual.h index 0bdd43ad..453071c7 100644 --- a/lib/algorithms/iterative/GeneralisedMinimalResidual.h +++ b/lib/algorithms/iterative/GeneralisedMinimalResidual.h @@ -29,7 +29,6 @@ directory #ifndef GRID_GENERALISED_MINIMAL_RESIDUAL_H #define GRID_GENERALISED_MINIMAL_RESIDUAL_H -/////////////////////////////////////////////////////////////////////////////////////////////////////// // from Y. Saad - Iterative Methods for Sparse Linear Systems, PP 172 // Compute r0 = b − Ax0 , β := ||r0||2 , and v1 := r0 /β // For j = 1, 2, ..., m Do: @@ -47,156 +46,274 @@ directory // want to solve Ax = b -> A = LinOp, psi = x, b = src -namespace Grid -{ -template< class Field > -class GeneralisedMinimalResidual : public OperatorFunction< Field > -{ -public: - bool ErrorOnNoConverge; // Throw an assert when GMRES fails to converge, - // defaults to True. - RealD Tolerance; - Integer MaxIterations; - Integer IterationsToComplete; // Number of iterations the GMRES took to - // finish. Filled in upon completion +namespace Grid { - GeneralisedMinimalResidual( RealD tol, - Integer maxit, - bool err_on_no_conv = true ) - : Tolerance( tol ) - , MaxIterations( maxit ) - , ErrorOnNoConverge( err_on_no_conv ){}; +template +class GeneralisedMinimalResidual : public OperatorFunction { + public: + bool ErrorOnNoConverge; // Throw an assert when GMRES fails to converge, + // defaults to True. + RealD Tolerance; + Integer MaxIterations; + Integer IterationsToComplete; // Number of iterations the GMRES took to + // finish. Filled in upon completion - // want to solve Ax = b -> A = LinOp, psi = x, b = src + GeneralisedMinimalResidual(RealD tol, + Integer maxit, + bool err_on_no_conv = true) + : Tolerance(tol), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv){}; - void operator()( LinearOperatorBase< Field > &LinOp, - const Field & src, - Field & psi ) - { - std::cout << GridLogMessage - << "GeneralisedMinimalResidual: Start of operator()" - << std::endl; - psi.checkerboard = src.checkerboard; - conformable( psi, src ); + // want to solve Ax = b -> A = LinOp, psi = x, b = src - Field r( src ); - Field mmv( src ); + /* void */ + /* operator()(LinearOperatorBase &LinOp, const Field &src, Field &psi) + * { */ + /* typedef typename Eigen::MatrixXcd MyMatrix; */ + /* typedef typename Eigen::VectorXcd MyVector; */ - std::vector< Field > v( MaxIterations + 1, src ); + /* Field r(src); */ + /* Field w(src); */ + /* Field mmv(src); */ - RealD beta{}; - RealD b{}; - RealD d{}; + /* std::vector V(MaxIterations + 1, src); */ + /* std::vector> y(MaxIterations + 1, 0.); */ + /* std::vector> gamma(MaxIterations + 1, 0.); */ + /* std::vector> c(MaxIterations + 1, 0.); */ + /* std::vector> s(MaxIterations + 1, 0.); */ - Eigen::MatrixXcd H - = Eigen::MatrixXcd::Zero( MaxIterations + 1, MaxIterations ); + /* int m = MaxIterations; */ - // Compute r0 = b − Ax0 , β := ||r0||2 , and v1 := r0 /β - LinOp.Op( psi, mmv ); + /* RealD gamma0{}; */ - r = src - mmv; - beta = norm2( r ); - V[ 0 ] = ( 1 / beta ) * r; + /* MyMatrix H = Eigen::MatrixXcd::Zero(MaxIterations + 1, MaxIterations); */ - for( auto j = 0; j < MaxIterations; ++j ) - { - LinOp.Op( V[ j ], mmv ); + /* RealD normPsiSq = norm2(psi); */ + /* RealD normSrcSq = norm2(src); */ + /* RealD TargetResSq = Tolerance * Tolerance * normSrcSq; */ - for( auto i = 0; i < j; ++i ) - { - std::cout - << GridLogMessage - << "GeneralisedMinimalResidual: End of inner iteration " - << i << std::endl; - H( i, j ) = innerProduct( mmv, v[ i ] ); - mmv = mmv - H( i, j ) * V[ i ]; - } + /* LinOp.Op(psi, mmv); */ - H( j + 1, j ) = norm2( mmv ); + /* r = src - mmv; */ + /* gamma[0] = norm2(r); */ + /* std::cout << gamma[0] << std::endl; */ + /* gamma0 = std::real(gamma[0]); */ + /* V[0] = (1. / gamma[0]) * r; */ - std::cout << GridLogMessage << "GeneralisedMinimalResidual: H" - << j + 1 << "," << j << "= " << H( j + 1, j ) - << std::endl; - if( H( j + 1, j ) == 0. ) - { - IterationsToComplete = j; - break; - } + /* std::cout << GridLogMessage << std::setprecision(4) */ + /* << "GeneralisedMinimalResidual: psi " << normPsiSq */ + /* << std::endl; */ + /* std::cout << GridLogMessage << std::setprecision(4) */ + /* << "GeneralisedMinimalResidual: src " << normSrcSq */ + /* << std::endl; */ + /* std::cout << GridLogMessage << std::setprecision(4) */ + /* << "GeneralisedMinimalResidual: target " << TargetResSq */ + /* << std::endl; */ + /* std::cout << GridLogMessage << std::setprecision(4) */ + /* << "GeneralisedMinimalResidual: r " << gamma0 << + * std::endl; */ - V[ j + 1 ] = ( 1. / H( j + 1, j ) ) * mmv; - std::cout << GridLogMessage - << "GeneralisedMinimalResidual: End of outer iteration " - << j << std::endl; - } - std::cout << GridLogMessage - << "GeneralisedMinimalResidual: End of operator()" - << std::endl; + /* std::cout */ + /* << GridLogIterative << std::setprecision(4) */ + /* << "GeneralisedMinimalResidual: before starting to iterate residual " + */ + /* << gamma0 << " target " << TargetResSq << std::endl; */ + + /* for(auto j = 0; j < m; ++j) { */ + /* LinOp.Op(V[j], w); */ + + /* for(auto i = 0; i <= j; ++i) { */ + /* H(i, j) = innerProduct(V[i], w); */ + /* w = w - H(i, j) * V[i]; */ + /* } */ + + /* H(j + 1, j) = norm2(w); */ + /* V[j + 1] = (1. / H(j + 1, j)) * w; */ + + /* if(std::abs(H(j + 1, j)) > 1e-15) { */ + /* qrUpdate(gamma, c, s, H, j); */ + /* } */ + + /* /\* std::cout << GridLogMessage << "GeneralisedMinimalResidual: H( " + * *\/ */ + /* /\* << j + 1 << "," << j << " ) = " << H( j + 1, j ) *\/ */ + /* /\* << std::endl; *\/ */ + + /* std::cout << GridLogIterative << "GeneralisedMinimalResidual: Iteration + * " */ + /* << j << " residual " << std::abs(gamma[j + 1]) << " target " + */ + /* << TargetResSq << std::endl; */ + /* if(std::abs(gamma[j + 1]) / gamma0 < Tolerance) { */ + /* IterationsToComplete = j; */ + /* break; */ + /* } */ + /* } */ + /* computeSolution(y, gamma, H, V, psi, IterationsToComplete); */ + /* std::cout << GridLogMessage */ + /* << "GeneralisedMinimalResidual: End of operator() after " */ + /* << IterationsToComplete << " iterations" << std::endl; */ + + /* RealD normSrc = sqrt(normSrcSq); */ + /* RealD resnorm = sqrt(norm2(mmv)); */ + /* RealD true_residual = resnorm / srcnorm; */ + /* Field result = mmv; */ + /* Field Dx(src); */ + /* Field tmp(src); */ + + /* // Test the correctness */ + /* LinOp.Op(result, Dx); */ + + /* tmp = Dx - src; */ + + /* std::cout << norm2(tmp) << " " << norm2(tmp) / gamma0 << std::endl; */ + /* } */ + + void + operator()(LinearOperatorBase &LinOp, const Field &src, Field &psi) { + std::cout << "GMRES: Start of operator()" << std::endl; + + int m = MaxIterations; + + Field r(src); + Field w(src); + Field Dpsi(src); + Field Dv(src); + + std::vector v(m + 1, src); + Eigen::MatrixXcd H = Eigen::MatrixXcd::Zero(m + 1, m); + + std::vector> y(m + 1, 0.); + std::vector> gamma(m + 1, 0.); + std::vector> c(m + 1, 0.); + std::vector> s(m + 1, 0.); + + LinOp.Op(psi, Dpsi); + r = src - Dpsi; + + RealD beta = norm2(r); + gamma[0] = beta; + + std::cout << "beta " << beta << std::endl; + + v[0] = (1. / beta) * r; + + // Begin iterating + for(auto j = 0; j < m; ++j) { + LinOp.Op(v[j], Dv); + w = Dv; + + for(auto i = 0; i <= j; ++i) { + H(i, j) = innerProduct(v[i], w); + w = w - H(i, j) * v[i]; + } + + H(j + 1, j) = norm2(w); + v[j + 1] = (1. / H(j + 1, j)) * w; + + // end of arnoldi process, begin of givens rotations + // apply old Givens rotation + for(auto i = 0; i < j ; ++i) { + auto tmp = -s[i] * H(i, j) + c[i] * H(i + 1, j); + H(i, j) = std::conj(c[i]) * H(i, j) + std::conj(s[i]) * H(i + 1, j); + H(i + 1, j) = tmp; + } + + // compute new Givens Rotation + ComplexD nu = sqrt(std::norm(H(j, j)) + std::norm(H(j + 1, j))); + c[j] = H(j, j) / nu; + s[j] = H(j + 1, j) / nu; + std::cout << "nu" << nu << std::endl; + std::cout << "H("<> &gamma, */ + /* std::vector> &c, */ + /* std::vector> &s, */ + /* Eigen::MatrixXcd & H, */ + /* int j) { */ + /* ComplexD beta{}; */ + /* // update QR factorization */ + /* // apply previous Givens rotation */ + /* for(auto i = 0; i < j; i++) { */ + /* beta = -s[i] * H(i, j) + c[i] * H(i + 1, j); */ + /* H(i, j) = std::conj(c[i]) * H(i, j) + std::conj(s[i]) * H(i + 1, + * j); */ + /* H(i + 1, j) = beta; */ + /* } */ + + /* // compute current Givens rotation */ + /* beta = sqrt(std::norm(H(j, j)) + std::norm(H(j + 1, j))); */ + /* s[j] = H(j + 1, j) / beta; */ + /* c[j] = H(j, j) / beta; */ + /* /\* std::cout << "beta= " << beta << std::endl; *\/ */ + /* /\* std::cout << "s[j]= " << s[ j ] << std::endl; *\/ */ + /* /\* std::cout << "c[j]= " << c[ j ] << std::endl; *\/ */ + + /* /\* std::cout << "gamma[j+1]= " << gamma[ j + 1 ] << std::endl; *\/ */ + /* /\* std::cout << "gamma[j]= " << gamma[ j ] << std::endl; *\/ */ + /* // update right column */ + /* gamma[j + 1] = -s[j] * gamma[j]; */ + /* gamma[j] = std::conj(c[j]) * gamma[j]; */ + /* /\* std::cout << "gamma[j+1]= " << gamma[ j + 1 ] << std::endl; *\/ */ + /* /\* std::cout << "gamma[j]= " << gamma[ j ] << std::endl; *\/ */ + + /* // apply current Givens rotation */ + /* H(j, j) = beta; */ + /* H(j + 1, j) = 0.; */ + /* /\* std::cout << "H(j,j)= " << H( j, j ) << std::endl; *\/ */ + /* /\* std::cout << "H(j+1,j)= " << H( j + 1, j ) << std::endl; *\/ */ + /* } */ + + void computeSolution(std::vector> & y, + std::vector> const &gamma, + Eigen::MatrixXcd const & H, + std::vector const & v, + Field & x, + int j) { + for(auto i = j; i >= 0; i--) { + y[i] = gamma[i]; + for(auto k = i + 1; k <= j; k++) + y[i] -= H(i, k) * y[k]; + y[i] /= H(i, i); + } + + /* if(true) // TODO ??? */ + /* { */ + /* for(auto i = 0; i <= j; i++) */ + /* x = x + v[i] * y[i]; */ + /* } else { */ + x = y[0] * v[0]; + for(auto i = 1; i <= j; i++) + x = x + v[i] * y[i]; + /* } */ + } }; } #endif - -// Note: The DD-αAMG codebase turns around the Hessenberg matrix - -void arnoldiStep() -{ - w = D * V[ j ]; - - for( auto i = 0; i <= j; ++i ) - H( i, j ) = innerProduct( V[ j + 1 ], w ); - w = w - H( i, j ) * V[ i ]; - - H( j + 1, j ) = norm2( w ); - - V[ j + 1 ] = w / H( j + 1, j ); -} - -void qr_update_PRECISION() -{ - // update QR factorization - // apply previous Givens rotation - for( auto i = 0; i < j; i++ ) - { - beta = -s[ i ] * H( i, j ) + c[ i ] * H( i + 1, j ); - H( i, j ) = std::conj( c[ i ] ) * H( i, j ) - + std::conj( s[ i ] ) * H( i + 1, j ); - H( i + 1, j ) = beta; - } - - // compute current Givens rotation - beta = sqrt( std::norm( H( j, j ) ) + std::norm( H( j, j + 1 ) ) ); - s[ j ] = H( j + 1, j ) / beta; - c[ j ] = H( j, j ) / beta; - - // update right column - gamma[ j + 1 ] = -s[ j ] * gamma[ j ]; - gamma[ j ] = std::conj( c[ j ] ) * gamma[ j ]; - - // apply current Givens rotation - H( j, j ) = beta; - H( j + 1, j ) = 0; -} - -// check -void compute_solution_PRECISION() -{ - for( auto i = j; i >= 0; i-- ) - { - y[ i ] = gamma[ i ]; - for( auto k = i + 1; k <= j; k++ ) - y[ i ] -= H( i, k ) * y[ k ]; - y[ i ] /= H( i, i ); - } - - if( true ) // TODO ??? - { - for( i = 0; i <= j; i++ ) - x = x + V[ i ] * y[ i ]; - } - else - { - x = y[ 0 ] * V[ 0 ]; - for( i = 1; i <= j; i++ ) - x = x + V[ i ] * y[ i ]; - } -}