mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-26 21:55:55 +01:00
Save current state
This commit is contained in:
Daniel Richtmann
Daniel Richtmann
parent
53cfa44d7a
commit
789e892865
@ -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 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
|
||||
|
||||
// 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<Field> &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<Field> V(MaxIterations + 1, src); */
|
||||
/* std::vector<std::complex<double>> y(MaxIterations + 1, 0.); */
|
||||
/* std::vector<std::complex<double>> gamma(MaxIterations + 1, 0.); */
|
||||
/* std::vector<std::complex<double>> c(MaxIterations + 1, 0.); */
|
||||
/* std::vector<std::complex<double>> 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<Field> &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<Field> v(m + 1, src);
|
||||
Eigen::MatrixXcd H = Eigen::MatrixXcd::Zero(m + 1, m);
|
||||
|
||||
std::vector<std::complex<double>> y(m + 1, 0.);
|
||||
std::vector<std::complex<double>> gamma(m + 1, 0.);
|
||||
std::vector<std::complex<double>> c(m + 1, 0.);
|
||||
std::vector<std::complex<double>> 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("<<j<<","<<j<<")" << H(j,j) << std::endl;
|
||||
std::cout << "H("<<j+1<<","<<j<<")" << H(j+1,j) << std::endl;
|
||||
|
||||
// apply new Givens rotation
|
||||
H(j, j) = nu;
|
||||
H(j + 1, j) = 0.;
|
||||
|
||||
/* ORDERING??? */
|
||||
gamma[j + 1] = -s[j] * gamma[j];
|
||||
gamma[j] = std::conj(c[j]) * gamma[j];
|
||||
|
||||
/* for(auto k = 0; k <= j+1 ; ++k) */
|
||||
/* std::cout << "k " << k << "nu " << nu << " c["<<k<<"]" << c[k]<< " s["<<k<<"]" << s[k] << " gamma["<<k<<"]" << gamma[k] << std::endl; */
|
||||
|
||||
std::cout << GridLogIterative << "GeneralisedMinimalResidual: Iteration "
|
||||
<< j << " residual " << std::abs(gamma[j + 1]) << std::endl; //" target "
|
||||
/* << TargetResSq << std::endl; */
|
||||
if(std::abs(gamma[j + 1]) / sqrt(beta) < Tolerance) {
|
||||
IterationsToComplete = j;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// backward substitution
|
||||
computeSolution(y, gamma, H, v, psi, IterationsToComplete);
|
||||
|
||||
std::cout << "GMRES: End of operator()" << std::endl;
|
||||
}
|
||||
|
||||
private:
|
||||
/* void qrUpdate(std::vector<std::complex<double>> &gamma, */
|
||||
/* std::vector<std::complex<double>> &c, */
|
||||
/* std::vector<std::complex<double>> &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<std::complex<double>> & y,
|
||||
std::vector<std::complex<double>> const &gamma,
|
||||
Eigen::MatrixXcd const & H,
|
||||
std::vector<Field> 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 ];
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user