mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
QMR implemented, preserve even if not used much
This commit is contained in:
Peter Boyle
parent
3d2fe80780
commit
58a31f0763
371
Grid/algorithms/iterative/QuasiMinimalResidual.h
Normal file
371
Grid/algorithms/iterative/QuasiMinimalResidual.h
Normal file
@ -0,0 +1,371 @@
|
|||||||
|
/*************************************************************************************
|
||||||
|
|
||||||
|
Grid physics library, www.github.com/paboyle/Grid
|
||||||
|
|
||||||
|
Source file: ./lib/algorithmsf/iterative/QuasiMinimalResidual.h
|
||||||
|
|
||||||
|
Copyright (C) 2019
|
||||||
|
|
||||||
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
|
|
||||||
|
This program is free software; you can redistribute it and/or modify
|
||||||
|
it under the terms of the GNU General Public License as published by
|
||||||
|
the Free Software Foundation; either version 2 of the License, or
|
||||||
|
(at your option) any later version.
|
||||||
|
|
||||||
|
This program is distributed in the hope that it will be useful,
|
||||||
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
|
GNU General Public License for more details.
|
||||||
|
|
||||||
|
You should have received a copy of the GNU General Public License along
|
||||||
|
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||||
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||||
|
|
||||||
|
See the full license in the file "LICENSE" in the top level distribution
|
||||||
|
directory
|
||||||
|
*************************************************************************************/
|
||||||
|
/* END LEGAL */
|
||||||
|
#pragma once
|
||||||
|
|
||||||
|
NAMESPACE_BEGIN(Grid);
|
||||||
|
|
||||||
|
template<class Field>
|
||||||
|
RealD innerG5ProductReal(Field &l, Field &r)
|
||||||
|
{
|
||||||
|
Gamma G5(Gamma::Algebra::Gamma5);
|
||||||
|
Field tmp(l.Grid());
|
||||||
|
// tmp = G5*r;
|
||||||
|
G5R5(tmp,r);
|
||||||
|
ComplexD ip =innerProduct(l,tmp);
|
||||||
|
std::cout << "innerProductRealG5R5 "<<ip<<std::endl;
|
||||||
|
return ip.real();
|
||||||
|
}
|
||||||
|
|
||||||
|
template<class Field>
|
||||||
|
class QuasiMinimalResidual : public OperatorFunction<Field> {
|
||||||
|
public:
|
||||||
|
using OperatorFunction<Field>::operator();
|
||||||
|
|
||||||
|
bool ErrorOnNoConverge;
|
||||||
|
RealD Tolerance;
|
||||||
|
Integer MaxIterations;
|
||||||
|
Integer IterationCount;
|
||||||
|
|
||||||
|
QuasiMinimalResidual(RealD tol,
|
||||||
|
Integer maxit,
|
||||||
|
bool err_on_no_conv = true)
|
||||||
|
: Tolerance(tol)
|
||||||
|
, MaxIterations(maxit)
|
||||||
|
, ErrorOnNoConverge(err_on_no_conv)
|
||||||
|
{};
|
||||||
|
|
||||||
|
#if 1
|
||||||
|
void operator()(LinearOperatorBase<Field> &LinOp, const Field &b, Field &x)
|
||||||
|
{
|
||||||
|
RealD resid;
|
||||||
|
IterationCount=0;
|
||||||
|
|
||||||
|
RealD rho, rho_1, xi, gamma, gamma_1, theta, theta_1;
|
||||||
|
RealD eta, delta, ep, beta;
|
||||||
|
|
||||||
|
GridBase *Grid = b.Grid();
|
||||||
|
Field r(Grid), d(Grid), s(Grid);
|
||||||
|
Field v(Grid), w(Grid), y(Grid), z(Grid);
|
||||||
|
Field v_tld(Grid), w_tld(Grid), y_tld(Grid), z_tld(Grid);
|
||||||
|
Field p(Grid), q(Grid), p_tld(Grid);
|
||||||
|
|
||||||
|
Real normb = norm2(b);
|
||||||
|
|
||||||
|
LinOp.Op(x,r); r = b - r;
|
||||||
|
|
||||||
|
assert(normb> 0.0);
|
||||||
|
|
||||||
|
resid = norm2(r)/normb;
|
||||||
|
if (resid <= Tolerance) {
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
v_tld = r;
|
||||||
|
y = v_tld;
|
||||||
|
rho = norm2(y);
|
||||||
|
|
||||||
|
// Take Gamma5 conjugate
|
||||||
|
// Gamma G5(Gamma::Algebra::Gamma5);
|
||||||
|
// G5R5(w_tld,r);
|
||||||
|
// w_tld = G5* v_tld;
|
||||||
|
w_tld=v_tld;
|
||||||
|
z = w_tld;
|
||||||
|
xi = norm2(z);
|
||||||
|
|
||||||
|
gamma = 1.0;
|
||||||
|
eta = -1.0;
|
||||||
|
theta = 0.0;
|
||||||
|
|
||||||
|
for (int i = 1; i <= MaxIterations; i++) {
|
||||||
|
|
||||||
|
// Breakdown tests
|
||||||
|
assert( rho != 0.0);
|
||||||
|
assert( xi != 0.0);
|
||||||
|
|
||||||
|
v = (1. / rho) * v_tld;
|
||||||
|
y = (1. / rho) * y;
|
||||||
|
|
||||||
|
w = (1. / xi) * w_tld;
|
||||||
|
z = (1. / xi) * z;
|
||||||
|
|
||||||
|
ComplexD Zdelta = innerProduct(z, y); // Complex?
|
||||||
|
std::cout << "Zdelta "<<Zdelta<<std::endl;
|
||||||
|
delta = Zdelta.real();
|
||||||
|
|
||||||
|
y_tld = y;
|
||||||
|
z_tld = z;
|
||||||
|
|
||||||
|
if (i > 1) {
|
||||||
|
p = y_tld - (xi * delta / ep) * p;
|
||||||
|
q = z_tld - (rho * delta / ep) * q;
|
||||||
|
} else {
|
||||||
|
p = y_tld;
|
||||||
|
q = z_tld;
|
||||||
|
}
|
||||||
|
|
||||||
|
LinOp.Op(p,p_tld); // p_tld = A * p;
|
||||||
|
ComplexD Zep = innerProduct(q, p_tld);
|
||||||
|
ep=Zep.real();
|
||||||
|
std::cout << "Zep "<<Zep <<std::endl;
|
||||||
|
// Complex Audit
|
||||||
|
assert(abs(ep)>0);
|
||||||
|
|
||||||
|
beta = ep / delta;
|
||||||
|
assert(abs(beta)>0);
|
||||||
|
|
||||||
|
v_tld = p_tld - beta * v;
|
||||||
|
y = v_tld;
|
||||||
|
|
||||||
|
rho_1 = rho;
|
||||||
|
rho = norm2(y);
|
||||||
|
LinOp.AdjOp(q,w_tld);
|
||||||
|
w_tld = w_tld - beta * w;
|
||||||
|
z = w_tld;
|
||||||
|
|
||||||
|
xi = norm2(z);
|
||||||
|
|
||||||
|
gamma_1 = gamma;
|
||||||
|
theta_1 = theta;
|
||||||
|
|
||||||
|
theta = rho / (gamma_1 * beta);
|
||||||
|
gamma = 1.0 / sqrt(1.0 + theta * theta);
|
||||||
|
std::cout << "theta "<<theta<<std::endl;
|
||||||
|
std::cout << "gamma "<<gamma<<std::endl;
|
||||||
|
|
||||||
|
assert(abs(gamma)> 0.0);
|
||||||
|
|
||||||
|
eta = -eta * rho_1 * gamma* gamma / (beta * gamma_1 * gamma_1);
|
||||||
|
|
||||||
|
if (i > 1) {
|
||||||
|
d = eta * p + (theta_1 * theta_1 * gamma * gamma) * d;
|
||||||
|
s = eta * p_tld + (theta_1 * theta_1 * gamma * gamma) * s;
|
||||||
|
} else {
|
||||||
|
d = eta * p;
|
||||||
|
s = eta * p_tld;
|
||||||
|
}
|
||||||
|
|
||||||
|
x =x+d; // update approximation vector
|
||||||
|
r =r-s; // compute residual
|
||||||
|
|
||||||
|
if ((resid = norm2(r) / normb) <= Tolerance) {
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
std::cout << "Iteration "<<i<<" resid " << resid<<std::endl;
|
||||||
|
}
|
||||||
|
assert(0);
|
||||||
|
return; // no convergence
|
||||||
|
}
|
||||||
|
#else
|
||||||
|
// QMRg5 SMP thesis
|
||||||
|
void operator()(LinearOperatorBase<Field> &LinOp, const Field &b, Field &x)
|
||||||
|
{
|
||||||
|
// Real scalars
|
||||||
|
GridBase *grid = b.Grid();
|
||||||
|
|
||||||
|
Field r(grid);
|
||||||
|
Field p_m(grid), p_m_minus_1(grid), p_m_minus_2(grid);
|
||||||
|
Field v_m(grid), v_m_minus_1(grid), v_m_plus_1(grid);
|
||||||
|
Field tmp(grid);
|
||||||
|
|
||||||
|
RealD w;
|
||||||
|
RealD z1, z2;
|
||||||
|
RealD delta_m, delta_m_minus_1;
|
||||||
|
RealD c_m_plus_1, c_m, c_m_minus_1;
|
||||||
|
RealD s_m_plus_1, s_m, s_m_minus_1;
|
||||||
|
RealD alpha, beta, gamma, epsilon;
|
||||||
|
RealD mu, nu, rho, theta, xi, chi;
|
||||||
|
RealD mod2r, mod2b;
|
||||||
|
RealD tau2, target2;
|
||||||
|
|
||||||
|
mod2b=norm2(b);
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// Initial residual
|
||||||
|
/////////////////////////
|
||||||
|
LinOp.Op(x,tmp);
|
||||||
|
r = b - tmp;
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// \mu = \rho = |r_0|
|
||||||
|
/////////////////////////
|
||||||
|
mod2r = norm2(r);
|
||||||
|
rho = sqrt( mod2r);
|
||||||
|
mu=rho;
|
||||||
|
|
||||||
|
std::cout << "QuasiMinimalResidual rho "<< rho<<std::endl;
|
||||||
|
/////////////////////////
|
||||||
|
// Zero negative history
|
||||||
|
/////////////////////////
|
||||||
|
v_m_plus_1 = Zero();
|
||||||
|
v_m_minus_1 = Zero();
|
||||||
|
p_m_minus_1 = Zero();
|
||||||
|
p_m_minus_2 = Zero();
|
||||||
|
|
||||||
|
// v0
|
||||||
|
v_m = (1.0/rho)*r;
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// Initial coeffs
|
||||||
|
/////////////////////////
|
||||||
|
delta_m_minus_1 = 1.0;
|
||||||
|
c_m_minus_1 = 1.0;
|
||||||
|
c_m = 1.0;
|
||||||
|
s_m_minus_1 = 0.0;
|
||||||
|
s_m = 0.0;
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// Set up convergence check
|
||||||
|
/////////////////////////
|
||||||
|
tau2 = mod2r;
|
||||||
|
target2 = mod2b * Tolerance*Tolerance;
|
||||||
|
|
||||||
|
for(int iter = 0 ; iter < MaxIterations; iter++){
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// \delta_m = (v_m, \gamma_5 v_m)
|
||||||
|
/////////////////////////
|
||||||
|
delta_m = innerG5ProductReal(v_m,v_m);
|
||||||
|
std::cout << "QuasiMinimalResidual delta_m "<< delta_m<<std::endl;
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// tmp = A v_m
|
||||||
|
/////////////////////////
|
||||||
|
LinOp.Op(v_m,tmp);
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// \alpha = (v_m, \gamma_5 temp) / \delta_m
|
||||||
|
/////////////////////////
|
||||||
|
alpha = innerG5ProductReal(v_m,tmp);
|
||||||
|
alpha = alpha/delta_m ;
|
||||||
|
std::cout << "QuasiMinimalResidual alpha "<< alpha<<std::endl;
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// \beta = \rho \delta_m / \delta_{m-1}
|
||||||
|
/////////////////////////
|
||||||
|
beta = rho * delta_m / delta_m_minus_1;
|
||||||
|
std::cout << "QuasiMinimalResidual beta "<< beta<<std::endl;
|
||||||
|
|
||||||
|
/////////////////////////
|
||||||
|
// \tilde{v}_{m+1} = temp - \alpha v_m - \beta v_{m-1}
|
||||||
|
/////////////////////////
|
||||||
|
v_m_plus_1 = tmp - alpha*v_m - beta*v_m_minus_1;
|
||||||
|
|
||||||
|
///////////////////////////////
|
||||||
|
// \rho = || \tilde{v}_{m+1} ||
|
||||||
|
///////////////////////////////
|
||||||
|
rho = sqrt( norm2(v_m_plus_1) );
|
||||||
|
std::cout << "QuasiMinimalResidual rho "<< rho<<std::endl;
|
||||||
|
|
||||||
|
///////////////////////////////
|
||||||
|
// v_{m+1} = \tilde{v}_{m+1}
|
||||||
|
///////////////////////////////
|
||||||
|
v_m_plus_1 = (1.0 / rho) * v_m_plus_1;
|
||||||
|
|
||||||
|
////////////////////////////////
|
||||||
|
// QMR recurrence coefficients.
|
||||||
|
////////////////////////////////
|
||||||
|
theta = s_m_minus_1 * beta;
|
||||||
|
gamma = c_m_minus_1 * beta;
|
||||||
|
epsilon = c_m * gamma + s_m * alpha;
|
||||||
|
xi = -s_m * gamma + c_m * alpha;
|
||||||
|
nu = sqrt( xi*xi + rho*rho );
|
||||||
|
c_m_plus_1 = fabs(xi) / nu;
|
||||||
|
if ( xi == 0.0 ) {
|
||||||
|
s_m_plus_1 = 1.0;
|
||||||
|
} else {
|
||||||
|
s_m_plus_1 = c_m_plus_1 * rho / xi;
|
||||||
|
}
|
||||||
|
chi = c_m_plus_1 * xi + s_m_plus_1 * rho;
|
||||||
|
|
||||||
|
std::cout << "QuasiMinimalResidual coeffs "<< theta <<" "<<gamma<<" "<< epsilon<<" "<< xi<<" "<< nu<<std::endl;
|
||||||
|
std::cout << "QuasiMinimalResidual coeffs "<< chi <<std::endl;
|
||||||
|
|
||||||
|
////////////////////////////////
|
||||||
|
//p_m=(v_m - \epsilon p_{m-1} - \theta p_{m-2}) / \chi
|
||||||
|
////////////////////////////////
|
||||||
|
p_m = (1.0/chi) * v_m - (epsilon/chi) * p_m_minus_1 - (theta/chi) * p_m_minus_2;
|
||||||
|
|
||||||
|
////////////////////////////////////////////////////////////////
|
||||||
|
// \psi = \psi + c_{m+1} \mu p_m
|
||||||
|
////////////////////////////////////////////////////////////////
|
||||||
|
x = x + ( c_m_plus_1 * mu ) * p_m;
|
||||||
|
|
||||||
|
////////////////////////////////////////
|
||||||
|
//
|
||||||
|
////////////////////////////////////////
|
||||||
|
mu = -s_m_plus_1 * mu;
|
||||||
|
delta_m_minus_1 = delta_m;
|
||||||
|
c_m_minus_1 = c_m;
|
||||||
|
c_m = c_m_plus_1;
|
||||||
|
s_m_minus_1 = s_m;
|
||||||
|
s_m = s_m_plus_1;
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Could use pointer swizzle games.
|
||||||
|
////////////////////////////////////
|
||||||
|
v_m_minus_1 = v_m;
|
||||||
|
v_m = v_m_plus_1;
|
||||||
|
p_m_minus_2 = p_m_minus_1;
|
||||||
|
p_m_minus_1 = p_m;
|
||||||
|
|
||||||
|
|
||||||
|
/////////////////////////////////////
|
||||||
|
// Convergence checks
|
||||||
|
/////////////////////////////////////
|
||||||
|
z1 = RealD(iter+1.0);
|
||||||
|
z2 = z1 + 1.0;
|
||||||
|
tau2 = tau2 *( z2 / z1 ) * s_m * s_m;
|
||||||
|
std::cout << " QuasiMinimumResidual iteration "<< iter<<std::endl;
|
||||||
|
std::cout << " QuasiMinimumResidual tau bound "<< tau2<<std::endl;
|
||||||
|
|
||||||
|
// Compute true residual
|
||||||
|
mod2r = tau2;
|
||||||
|
if ( 1 || (tau2 < (100.0 * target2)) ) {
|
||||||
|
LinOp.Op(x,tmp);
|
||||||
|
r = b - tmp;
|
||||||
|
mod2r = norm2(r);
|
||||||
|
std::cout << " QuasiMinimumResidual true residual is "<< mod2r<<std::endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
if ( mod2r < target2 ) {
|
||||||
|
|
||||||
|
std::cout << " QuasiMinimumResidual has converged"<<std::endl;
|
||||||
|
return;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
};
|
||||||
|
|
||||||
|
NAMESPACE_END(Grid);
|
||||||
Loading…
Reference in New Issue
Block a user