mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 15:55:37 +00:00
e9cc21900f
gives 8x (!) speed up on Haswell laptop vs. standard CG for 8 RHS solves. 166 iterations vs. 537 iterations so algorithmic gain + 2x in flop rate gain. Better than a slap in the face with a wet kipper.
597 lines
19 KiB
C++
597 lines
19 KiB
C++
/*************************************************************************************
|
|
|
|
Grid physics library, www.github.com/paboyle/Grid
|
|
|
|
Source file: ./lib/algorithms/iterative/BlockConjugateGradient.h
|
|
|
|
Copyright (C) 2017
|
|
|
|
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
|
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 */
|
|
#ifndef GRID_BLOCK_CONJUGATE_GRADIENT_H
|
|
#define GRID_BLOCK_CONJUGATE_GRADIENT_H
|
|
|
|
|
|
namespace Grid {
|
|
|
|
enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS };
|
|
|
|
//////////////////////////////////////////////////////////////////////////
|
|
// Block conjugate gradient. Dimension zero should be the block direction
|
|
//////////////////////////////////////////////////////////////////////////
|
|
template <class Field>
|
|
class BlockConjugateGradient : public OperatorFunction<Field> {
|
|
public:
|
|
|
|
|
|
typedef typename Field::scalar_type scomplex;
|
|
|
|
int blockDim ;
|
|
int Nblock;
|
|
|
|
BlockCGtype CGtype;
|
|
bool ErrorOnNoConverge; // throw an assert when the CG fails to converge.
|
|
// Defaults true.
|
|
RealD Tolerance;
|
|
Integer MaxIterations;
|
|
Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
|
|
|
|
BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true)
|
|
: Tolerance(tol),
|
|
CGtype(cgtype),
|
|
blockDim(_Orthog),
|
|
MaxIterations(maxit),
|
|
ErrorOnNoConverge(err_on_no_conv){};
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Thin QR factorisation (google it)
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
void ThinQRfact (Eigen::MatrixXcd &m_rr,
|
|
Eigen::MatrixXcd &C,
|
|
Eigen::MatrixXcd &Cinv,
|
|
Field & Q,
|
|
const Field & R)
|
|
{
|
|
int Orthog = blockDim; // First dimension is block dim; this is an assumption
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
//Dimensions
|
|
// R_{ferm x Nblock} = Q_{ferm x Nblock} x C_{Nblock x Nblock} -> ferm x Nblock
|
|
//
|
|
// Rdag R = m_rr = Herm = L L^dag <-- Cholesky decomposition (LLT routine in Eigen)
|
|
//
|
|
// Q C = R => Q = R C^{-1}
|
|
//
|
|
// Want Ident = Q^dag Q = C^{-dag} R^dag R C^{-1} = C^{-dag} L L^dag C^{-1} = 1_{Nblock x Nblock}
|
|
//
|
|
// Set C = L^{dag}, and then Q^dag Q = ident
|
|
//
|
|
// Checks:
|
|
// Cdag C = Rdag R ; passes.
|
|
// QdagQ = 1 ; passes
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
sliceInnerProductMatrix(m_rr,R,R,Orthog);
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Cholesky from Eigen
|
|
// There exists a ldlt that is documented as more stable
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
Eigen::MatrixXcd L = m_rr.llt().matrixL();
|
|
|
|
C = L.adjoint();
|
|
Cinv = C.inverse();
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Q = R C^{-1}
|
|
//
|
|
// Q_j = R_i Cinv(i,j)
|
|
//
|
|
// NB maddMatrix conventions are Right multiplication X[j] a[j,i] already
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
// FIXME:: make a sliceMulMatrix to avoid zero vector
|
|
sliceMulMatrix(Q,Cinv,R,Orthog);
|
|
}
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Call one of several implementations
|
|
////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
void operator()(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
|
|
{
|
|
if ( CGtype == BlockCGrQ ) {
|
|
BlockCGrQsolve(Linop,Src,Psi);
|
|
} else if (CGtype == BlockCG ) {
|
|
BlockCGsolve(Linop,Src,Psi);
|
|
} else if (CGtype == CGmultiRHS ) {
|
|
CGmultiRHSsolve(Linop,Src,Psi);
|
|
} else {
|
|
assert(0);
|
|
}
|
|
}
|
|
|
|
////////////////////////////////////////////////////////////////////////////
|
|
// BlockCGrQ implementation:
|
|
//--------------------------
|
|
// X is guess/Solution
|
|
// B is RHS
|
|
// Solve A X_i = B_i ; i refers to Nblock index
|
|
////////////////////////////////////////////////////////////////////////////
|
|
void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
|
|
{
|
|
int Orthog = blockDim; // First dimension is block dim; this is an assumption
|
|
Nblock = B._grid->_fdimensions[Orthog];
|
|
|
|
std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
|
|
|
|
X.checkerboard = B.checkerboard;
|
|
conformable(X, B);
|
|
|
|
Field tmp(B);
|
|
Field Q(B);
|
|
Field D(B);
|
|
Field Z(B);
|
|
Field AD(B);
|
|
|
|
Eigen::MatrixXcd m_DZ = Eigen::MatrixXcd::Identity(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_M = Eigen::MatrixXcd::Identity(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_rr = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
|
|
Eigen::MatrixXcd m_C = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_Cinv = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_S = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_Sinv = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
|
|
Eigen::MatrixXcd m_tmp = Eigen::MatrixXcd::Identity(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_tmp1 = Eigen::MatrixXcd::Identity(Nblock,Nblock);
|
|
|
|
// Initial residual computation & set up
|
|
std::vector<RealD> residuals(Nblock);
|
|
std::vector<RealD> ssq(Nblock);
|
|
|
|
sliceNorm(ssq,B,Orthog);
|
|
RealD sssum=0;
|
|
for(int b=0;b<Nblock;b++) sssum+=ssq[b];
|
|
|
|
sliceNorm(residuals,B,Orthog);
|
|
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
|
|
|
sliceNorm(residuals,X,Orthog);
|
|
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
|
|
|
/************************************************************************
|
|
* Block conjugate gradient rQ (Sebastien Birk Thesis, after Dubrulle 2001)
|
|
************************************************************************
|
|
* Dimensions:
|
|
*
|
|
* X,B==(Nferm x Nblock)
|
|
* A==(Nferm x Nferm)
|
|
*
|
|
* Nferm = Nspin x Ncolour x Ncomplex x Nlattice_site
|
|
*
|
|
* QC = R = B-AX, D = Q ; QC => Thin QR factorisation (google it)
|
|
* for k:
|
|
* Z = AD
|
|
* M = [D^dag Z]^{-1}
|
|
* X = X + D MC
|
|
* QS = Q - ZM
|
|
* D = Q + D S^dag
|
|
* C = S C
|
|
*/
|
|
///////////////////////////////////////
|
|
// Initial block: initial search dir is guess
|
|
///////////////////////////////////////
|
|
std::cout << GridLogMessage<<"BlockCGrQ algorithm initialisation " <<std::endl;
|
|
|
|
//1. QC = R = B-AX, D = Q ; QC => Thin QR factorisation (google it)
|
|
|
|
Linop.HermOp(X, AD);
|
|
tmp = B - AD;
|
|
ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);
|
|
D=Q;
|
|
|
|
std::cout << GridLogMessage<<"BlockCGrQ computed initial residual and QR fact " <<std::endl;
|
|
|
|
///////////////////////////////////////
|
|
// Timers
|
|
///////////////////////////////////////
|
|
GridStopWatch sliceInnerTimer;
|
|
GridStopWatch sliceMaddTimer;
|
|
GridStopWatch QRTimer;
|
|
GridStopWatch MatrixTimer;
|
|
GridStopWatch SolverTimer;
|
|
SolverTimer.Start();
|
|
|
|
int k;
|
|
for (k = 1; k <= MaxIterations; k++) {
|
|
|
|
//3. Z = AD
|
|
MatrixTimer.Start();
|
|
Linop.HermOp(D, Z);
|
|
MatrixTimer.Stop();
|
|
|
|
//4. M = [D^dag Z]^{-1}
|
|
sliceInnerTimer.Start();
|
|
sliceInnerProductMatrix(m_DZ,D,Z,Orthog);
|
|
sliceInnerTimer.Stop();
|
|
m_M = m_DZ.inverse();
|
|
|
|
//5. X = X + D MC
|
|
m_tmp = m_M * m_C;
|
|
sliceMaddTimer.Start();
|
|
sliceMaddMatrix(X,m_tmp, D,X,Orthog);
|
|
sliceMaddTimer.Stop();
|
|
|
|
//6. QS = Q - ZM
|
|
sliceMaddTimer.Start();
|
|
sliceMaddMatrix(tmp,m_M,Z,Q,Orthog,-1.0);
|
|
sliceMaddTimer.Stop();
|
|
QRTimer.Start();
|
|
ThinQRfact (m_rr, m_S, m_Sinv, Q, tmp);
|
|
QRTimer.Stop();
|
|
|
|
//7. D = Q + D S^dag
|
|
m_tmp = m_S.adjoint();
|
|
sliceMaddTimer.Start();
|
|
sliceMaddMatrix(D,m_tmp,D,Q,Orthog);
|
|
sliceMaddTimer.Stop();
|
|
|
|
//8. C = S C
|
|
m_C = m_S*m_C;
|
|
|
|
/*********************
|
|
* convergence monitor
|
|
*********************
|
|
*/
|
|
m_rr = m_C.adjoint() * m_C;
|
|
|
|
RealD max_resid=0;
|
|
RealD rrsum=0;
|
|
RealD rr;
|
|
|
|
for(int b=0;b<Nblock;b++) {
|
|
rrsum+=real(m_rr(b,b));
|
|
rr = real(m_rr(b,b))/ssq[b];
|
|
if ( rr > max_resid ) max_resid = rr;
|
|
}
|
|
|
|
std::cout << GridLogIterative << "\titeration "<<k<<" rr_sum "<<rrsum<<" ssq_sum "<< sssum
|
|
<<" ave "<<std::sqrt(rrsum/sssum) << " max "<< max_resid <<std::endl;
|
|
|
|
if ( max_resid < Tolerance*Tolerance ) {
|
|
|
|
SolverTimer.Stop();
|
|
|
|
std::cout << GridLogMessage<<"BlockCGrQ converged in "<<k<<" iterations"<<std::endl;
|
|
|
|
for(int b=0;b<Nblock;b++){
|
|
std::cout << GridLogMessage<< "\t\tblock "<<b<<" computed resid "
|
|
<< std::sqrt(real(m_rr(b,b))/ssq[b])<<std::endl;
|
|
}
|
|
std::cout << GridLogMessage<<"\tMax residual is "<<std::sqrt(max_resid)<<std::endl;
|
|
|
|
Linop.HermOp(X, AD);
|
|
AD = AD-B;
|
|
std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(norm2(AD)/norm2(B)) <<std::endl;
|
|
|
|
std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
|
|
std::cout << GridLogMessage << "\tElapsed " << SolverTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tInnerProd " << sliceInnerTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tMaddMatrix " << sliceMaddTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tThinQRfact " << QRTimer.Elapsed() <<std::endl;
|
|
|
|
IterationsToComplete = k;
|
|
return;
|
|
}
|
|
|
|
}
|
|
std::cout << GridLogMessage << "BlockConjugateGradient(rQ) did NOT converge" << std::endl;
|
|
|
|
if (ErrorOnNoConverge) assert(0);
|
|
IterationsToComplete = k;
|
|
}
|
|
//////////////////////////////////////////////////////////////////////////
|
|
// Block conjugate gradient; Original O'Leary Dimension zero should be the block direction
|
|
//////////////////////////////////////////////////////////////////////////
|
|
void BlockCGsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
|
|
{
|
|
int Orthog = blockDim; // First dimension is block dim; this is an assumption
|
|
Nblock = Src._grid->_fdimensions[Orthog];
|
|
|
|
std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
|
|
|
|
Psi.checkerboard = Src.checkerboard;
|
|
conformable(Psi, Src);
|
|
|
|
Field P(Src);
|
|
Field AP(Src);
|
|
Field R(Src);
|
|
|
|
Eigen::MatrixXcd m_pAp = Eigen::MatrixXcd::Identity(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_pAp_inv= Eigen::MatrixXcd::Identity(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_rr = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_rr_inv = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
|
|
Eigen::MatrixXcd m_alpha = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
Eigen::MatrixXcd m_beta = Eigen::MatrixXcd::Zero(Nblock,Nblock);
|
|
|
|
// Initial residual computation & set up
|
|
std::vector<RealD> residuals(Nblock);
|
|
std::vector<RealD> ssq(Nblock);
|
|
|
|
sliceNorm(ssq,Src,Orthog);
|
|
RealD sssum=0;
|
|
for(int b=0;b<Nblock;b++) sssum+=ssq[b];
|
|
|
|
sliceNorm(residuals,Src,Orthog);
|
|
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
|
|
|
sliceNorm(residuals,Psi,Orthog);
|
|
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
|
|
|
// Initial search dir is guess
|
|
Linop.HermOp(Psi, AP);
|
|
|
|
|
|
/************************************************************************
|
|
* Block conjugate gradient (Stephen Pickles, thesis 1995, pp 71, O Leary 1980)
|
|
************************************************************************
|
|
* O'Leary : R = B - A X
|
|
* O'Leary : P = M R ; preconditioner M = 1
|
|
* O'Leary : alpha = PAP^{-1} RMR
|
|
* O'Leary : beta = RMR^{-1}_old RMR_new
|
|
* O'Leary : X=X+Palpha
|
|
* O'Leary : R_new=R_old-AP alpha
|
|
* O'Leary : P=MR_new+P beta
|
|
*/
|
|
|
|
R = Src - AP;
|
|
P = R;
|
|
sliceInnerProductMatrix(m_rr,R,R,Orthog);
|
|
|
|
GridStopWatch sliceInnerTimer;
|
|
GridStopWatch sliceMaddTimer;
|
|
GridStopWatch MatrixTimer;
|
|
GridStopWatch SolverTimer;
|
|
SolverTimer.Start();
|
|
|
|
int k;
|
|
for (k = 1; k <= MaxIterations; k++) {
|
|
|
|
RealD rrsum=0;
|
|
for(int b=0;b<Nblock;b++) rrsum+=real(m_rr(b,b));
|
|
|
|
std::cout << GridLogIterative << "\titeration "<<k<<" rr_sum "<<rrsum<<" ssq_sum "<< sssum
|
|
<<" / "<<std::sqrt(rrsum/sssum) <<std::endl;
|
|
|
|
MatrixTimer.Start();
|
|
Linop.HermOp(P, AP);
|
|
MatrixTimer.Stop();
|
|
|
|
// Alpha
|
|
sliceInnerTimer.Start();
|
|
sliceInnerProductMatrix(m_pAp,P,AP,Orthog);
|
|
sliceInnerTimer.Stop();
|
|
m_pAp_inv = m_pAp.inverse();
|
|
m_alpha = m_pAp_inv * m_rr ;
|
|
|
|
// Psi, R update
|
|
sliceMaddTimer.Start();
|
|
sliceMaddMatrix(Psi,m_alpha, P,Psi,Orthog); // add alpha * P to psi
|
|
sliceMaddMatrix(R ,m_alpha,AP, R,Orthog,-1.0);// sub alpha * AP to resid
|
|
sliceMaddTimer.Stop();
|
|
|
|
// Beta
|
|
m_rr_inv = m_rr.inverse();
|
|
sliceInnerTimer.Start();
|
|
sliceInnerProductMatrix(m_rr,R,R,Orthog);
|
|
sliceInnerTimer.Stop();
|
|
m_beta = m_rr_inv *m_rr;
|
|
|
|
// Search update
|
|
sliceMaddTimer.Start();
|
|
sliceMaddMatrix(AP,m_beta,P,R,Orthog);
|
|
sliceMaddTimer.Stop();
|
|
P= AP;
|
|
|
|
/*********************
|
|
* convergence monitor
|
|
*********************
|
|
*/
|
|
RealD max_resid=0;
|
|
RealD rr;
|
|
for(int b=0;b<Nblock;b++){
|
|
rr = real(m_rr(b,b))/ssq[b];
|
|
if ( rr > max_resid ) max_resid = rr;
|
|
}
|
|
|
|
if ( max_resid < Tolerance*Tolerance ) {
|
|
|
|
SolverTimer.Stop();
|
|
|
|
std::cout << GridLogMessage<<"BlockCG converged in "<<k<<" iterations"<<std::endl;
|
|
for(int b=0;b<Nblock;b++){
|
|
std::cout << GridLogMessage<< "\t\tblock "<<b<<" computed resid "
|
|
<< std::sqrt(real(m_rr(b,b))/ssq[b])<<std::endl;
|
|
}
|
|
std::cout << GridLogMessage<<"\tMax residual is "<<std::sqrt(max_resid)<<std::endl;
|
|
|
|
Linop.HermOp(Psi, AP);
|
|
AP = AP-Src;
|
|
std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(norm2(AP)/norm2(Src)) <<std::endl;
|
|
|
|
std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
|
|
std::cout << GridLogMessage << "\tElapsed " << SolverTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tInnerProd " << sliceInnerTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tMaddMatrix " << sliceMaddTimer.Elapsed() <<std::endl;
|
|
|
|
IterationsToComplete = k;
|
|
return;
|
|
}
|
|
|
|
}
|
|
std::cout << GridLogMessage << "BlockConjugateGradient did NOT converge" << std::endl;
|
|
|
|
if (ErrorOnNoConverge) assert(0);
|
|
IterationsToComplete = k;
|
|
}
|
|
//////////////////////////////////////////////////////////////////////////
|
|
// multiRHS conjugate gradient. Dimension zero should be the block direction
|
|
// Use this for spread out across nodes
|
|
//////////////////////////////////////////////////////////////////////////
|
|
void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
|
|
{
|
|
int Orthog = blockDim; // First dimension is block dim
|
|
Nblock = Src._grid->_fdimensions[Orthog];
|
|
|
|
std::cout<<GridLogMessage<<"MultiRHS Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
|
|
|
|
Psi.checkerboard = Src.checkerboard;
|
|
conformable(Psi, Src);
|
|
|
|
Field P(Src);
|
|
Field AP(Src);
|
|
Field R(Src);
|
|
|
|
std::vector<ComplexD> v_pAp(Nblock);
|
|
std::vector<RealD> v_rr (Nblock);
|
|
std::vector<RealD> v_rr_inv(Nblock);
|
|
std::vector<RealD> v_alpha(Nblock);
|
|
std::vector<RealD> v_beta(Nblock);
|
|
|
|
// Initial residual computation & set up
|
|
std::vector<RealD> residuals(Nblock);
|
|
std::vector<RealD> ssq(Nblock);
|
|
|
|
sliceNorm(ssq,Src,Orthog);
|
|
RealD sssum=0;
|
|
for(int b=0;b<Nblock;b++) sssum+=ssq[b];
|
|
|
|
sliceNorm(residuals,Src,Orthog);
|
|
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
|
|
|
sliceNorm(residuals,Psi,Orthog);
|
|
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
|
|
|
// Initial search dir is guess
|
|
Linop.HermOp(Psi, AP);
|
|
|
|
R = Src - AP;
|
|
P = R;
|
|
sliceNorm(v_rr,R,Orthog);
|
|
|
|
GridStopWatch sliceInnerTimer;
|
|
GridStopWatch sliceMaddTimer;
|
|
GridStopWatch sliceNormTimer;
|
|
GridStopWatch MatrixTimer;
|
|
GridStopWatch SolverTimer;
|
|
|
|
SolverTimer.Start();
|
|
int k;
|
|
for (k = 1; k <= MaxIterations; k++) {
|
|
|
|
RealD rrsum=0;
|
|
for(int b=0;b<Nblock;b++) rrsum+=real(v_rr[b]);
|
|
|
|
std::cout << GridLogIterative << "\titeration "<<k<<" rr_sum "<<rrsum<<" ssq_sum "<< sssum
|
|
<<" / "<<std::sqrt(rrsum/sssum) <<std::endl;
|
|
|
|
MatrixTimer.Start();
|
|
Linop.HermOp(P, AP);
|
|
MatrixTimer.Stop();
|
|
|
|
// Alpha
|
|
sliceInnerTimer.Start();
|
|
sliceInnerProductVector(v_pAp,P,AP,Orthog);
|
|
sliceInnerTimer.Stop();
|
|
for(int b=0;b<Nblock;b++){
|
|
v_alpha[b] = v_rr[b]/real(v_pAp[b]);
|
|
}
|
|
|
|
// Psi, R update
|
|
sliceMaddTimer.Start();
|
|
sliceMaddVector(Psi,v_alpha, P,Psi,Orthog); // add alpha * P to psi
|
|
sliceMaddVector(R ,v_alpha,AP, R,Orthog,-1.0);// sub alpha * AP to resid
|
|
sliceMaddTimer.Stop();
|
|
|
|
// Beta
|
|
for(int b=0;b<Nblock;b++){
|
|
v_rr_inv[b] = 1.0/v_rr[b];
|
|
}
|
|
sliceNormTimer.Start();
|
|
sliceNorm(v_rr,R,Orthog);
|
|
sliceNormTimer.Stop();
|
|
for(int b=0;b<Nblock;b++){
|
|
v_beta[b] = v_rr_inv[b] *v_rr[b];
|
|
}
|
|
|
|
// Search update
|
|
sliceMaddTimer.Start();
|
|
sliceMaddVector(P,v_beta,P,R,Orthog);
|
|
sliceMaddTimer.Stop();
|
|
|
|
/*********************
|
|
* convergence monitor
|
|
*********************
|
|
*/
|
|
RealD max_resid=0;
|
|
for(int b=0;b<Nblock;b++){
|
|
RealD rr = v_rr[b]/ssq[b];
|
|
if ( rr > max_resid ) max_resid = rr;
|
|
}
|
|
|
|
if ( max_resid < Tolerance*Tolerance ) {
|
|
|
|
SolverTimer.Stop();
|
|
|
|
std::cout << GridLogMessage<<"MultiRHS solver converged in " <<k<<" iterations"<<std::endl;
|
|
for(int b=0;b<Nblock;b++){
|
|
std::cout << GridLogMessage<< "\t\tBlock "<<b<<" computed resid "<< std::sqrt(v_rr[b]/ssq[b])<<std::endl;
|
|
}
|
|
std::cout << GridLogMessage<<"\tMax residual is "<<std::sqrt(max_resid)<<std::endl;
|
|
|
|
Linop.HermOp(Psi, AP);
|
|
AP = AP-Src;
|
|
std::cout <<GridLogMessage << "\tTrue residual is " << std::sqrt(norm2(AP)/norm2(Src)) <<std::endl;
|
|
|
|
std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
|
|
std::cout << GridLogMessage << "\tElapsed " << SolverTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tInnerProd " << sliceInnerTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tNorm " << sliceNormTimer.Elapsed() <<std::endl;
|
|
std::cout << GridLogMessage << "\tMaddMatrix " << sliceMaddTimer.Elapsed() <<std::endl;
|
|
|
|
|
|
IterationsToComplete = k;
|
|
return;
|
|
}
|
|
|
|
}
|
|
std::cout << GridLogMessage << "MultiRHSConjugateGradient did NOT converge" << std::endl;
|
|
|
|
if (ErrorOnNoConverge) assert(0);
|
|
IterationsToComplete = k;
|
|
}
|
|
|
|
};
|
|
|
|
}
|
|
#endif
|