mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Block CG improvements to develop
This commit is contained in:
Azusa Yamaguchi
parent
8f514ae550
commit
4a47b11876
@ -33,7 +33,7 @@ directory
|
||||
|
||||
namespace Grid {
|
||||
|
||||
enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS };
|
||||
enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS, BlockCGVec, BlockCGrQVec };
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////
|
||||
// Block conjugate gradient. Dimension zero should be the block direction
|
||||
@ -42,7 +42,6 @@ template <class Field>
|
||||
class BlockConjugateGradient : public OperatorFunction<Field> {
|
||||
public:
|
||||
|
||||
|
||||
typedef typename Field::scalar_type scomplex;
|
||||
|
||||
int blockDim ;
|
||||
@ -54,21 +53,15 @@ class BlockConjugateGradient : public OperatorFunction<Field> {
|
||||
RealD Tolerance;
|
||||
Integer MaxIterations;
|
||||
Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
|
||||
Integer PrintInterval; //GridLogMessages or Iterative
|
||||
|
||||
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)
|
||||
: Tolerance(tol), CGtype(cgtype), blockDim(_Orthog), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv),PrintInterval(100)
|
||||
{};
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// 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
|
||||
@ -85,22 +78,20 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr,
|
||||
// Cdag C = Rdag R ; passes.
|
||||
// QdagQ = 1 ; passes
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
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
|
||||
sliceInnerProductMatrix(m_rr,R,R,Orthog);
|
||||
|
||||
// Force manifest hermitian to avoid rounding related
|
||||
m_rr = 0.5*(m_rr+m_rr.adjoint());
|
||||
|
||||
#if 0
|
||||
std::cout << " Calling Cholesky ldlt on m_rr " << m_rr <<std::endl;
|
||||
Eigen::MatrixXcd L_ldlt = m_rr.ldlt().matrixL();
|
||||
std::cout << " Called Cholesky ldlt on m_rr " << L_ldlt <<std::endl;
|
||||
auto D_ldlt = m_rr.ldlt().vectorD();
|
||||
std::cout << " Called Cholesky ldlt on m_rr " << D_ldlt <<std::endl;
|
||||
#endif
|
||||
|
||||
// std::cout << " Calling Cholesky llt on m_rr " <<std::endl;
|
||||
Eigen::MatrixXcd L = m_rr.llt().matrixL();
|
||||
// std::cout << " Called Cholesky llt on m_rr " << L <<std::endl;
|
||||
|
||||
C = L.adjoint();
|
||||
Cinv = C.inverse();
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
@ -112,6 +103,25 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr,
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
sliceMulMatrix(Q,Cinv,R,Orthog);
|
||||
}
|
||||
// see comments above
|
||||
void ThinQRfact (Eigen::MatrixXcd &m_rr,
|
||||
Eigen::MatrixXcd &C,
|
||||
Eigen::MatrixXcd &Cinv,
|
||||
std::vector<Field> & Q,
|
||||
const std::vector<Field> & R)
|
||||
{
|
||||
InnerProductMatrix(m_rr,R,R);
|
||||
|
||||
m_rr = 0.5*(m_rr+m_rr.adjoint());
|
||||
|
||||
Eigen::MatrixXcd L = m_rr.llt().matrixL();
|
||||
|
||||
C = L.adjoint();
|
||||
Cinv = C.inverse();
|
||||
|
||||
MulMatrix(Q,Cinv,R);
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Call one of several implementations
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
@ -119,14 +129,20 @@ 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);
|
||||
}
|
||||
}
|
||||
void operator()(LinearOperatorBase<Field> &Linop, const std::vector<Field> &Src, std::vector<Field> &Psi)
|
||||
{
|
||||
if ( CGtype == BlockCGrQVec ) {
|
||||
BlockCGrQsolveVec(Linop,Src,Psi);
|
||||
} else {
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
// BlockCGrQ implementation:
|
||||
@ -139,7 +155,8 @@ 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];
|
||||
|
||||
/* FAKE */
|
||||
Nblock=8;
|
||||
std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
|
||||
|
||||
X.checkerboard = B.checkerboard;
|
||||
@ -202,15 +219,10 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
|
||||
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;
|
||||
//std::cout << GridLogMessage << " initial tmp " << norm2(tmp)<< std::endl;
|
||||
|
||||
ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);
|
||||
//std::cout << GridLogMessage << " initial Q " << norm2(Q)<< std::endl;
|
||||
//std::cout << GridLogMessage << " m_rr " << m_rr<<std::endl;
|
||||
//std::cout << GridLogMessage << " m_C " << m_C<<std::endl;
|
||||
//std::cout << GridLogMessage << " m_Cinv " << m_Cinv<<std::endl;
|
||||
D=Q;
|
||||
|
||||
std::cout << GridLogMessage<<"BlockCGrQ computed initial residual and QR fact " <<std::endl;
|
||||
@ -232,14 +244,12 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
|
||||
MatrixTimer.Start();
|
||||
Linop.HermOp(D, Z);
|
||||
MatrixTimer.Stop();
|
||||
//std::cout << GridLogMessage << " norm2 Z " <<norm2(Z)<<std::endl;
|
||||
|
||||
//4. M = [D^dag Z]^{-1}
|
||||
sliceInnerTimer.Start();
|
||||
sliceInnerProductMatrix(m_DZ,D,Z,Orthog);
|
||||
sliceInnerTimer.Stop();
|
||||
m_M = m_DZ.inverse();
|
||||
//std::cout << GridLogMessage << " m_DZ " <<m_DZ<<std::endl;
|
||||
|
||||
//5. X = X + D MC
|
||||
m_tmp = m_M * m_C;
|
||||
@ -257,6 +267,7 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
|
||||
|
||||
//7. D = Q + D S^dag
|
||||
m_tmp = m_S.adjoint();
|
||||
|
||||
sliceMaddTimer.Start();
|
||||
sliceMaddMatrix(D,m_tmp,D,Q,Orthog);
|
||||
sliceMaddTimer.Stop();
|
||||
@ -317,152 +328,6 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
|
||||
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
|
||||
//////////////////////////////////////////////////////////////////////////
|
||||
@ -600,6 +465,233 @@ void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &
|
||||
IterationsToComplete = k;
|
||||
}
|
||||
|
||||
void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector<Field> &X, const std::vector<Field> &Y){
|
||||
for(int b=0;b<Nblock;b++){
|
||||
for(int bp=0;bp<Nblock;bp++) {
|
||||
m(b,bp) = innerProduct(X[b],Y[bp]);
|
||||
}}
|
||||
}
|
||||
void MaddMatrix(std::vector<Field> &AP, Eigen::MatrixXcd &m , const std::vector<Field> &X,const std::vector<Field> &Y,RealD scale=1.0){
|
||||
// Should make this cache friendly with site outermost, parallel_for
|
||||
// Deal with case AP aliases with either Y or X
|
||||
std::vector<Field> tmp(Nblock,X[0]);
|
||||
for(int b=0;b<Nblock;b++){
|
||||
tmp[b] = Y[b];
|
||||
for(int bp=0;bp<Nblock;bp++) {
|
||||
tmp[b] = tmp[b] + (scale*m(bp,b))*X[bp];
|
||||
}
|
||||
}
|
||||
for(int b=0;b<Nblock;b++){
|
||||
AP[b] = tmp[b];
|
||||
}
|
||||
}
|
||||
void MulMatrix(std::vector<Field> &AP, Eigen::MatrixXcd &m , const std::vector<Field> &X){
|
||||
// Should make this cache friendly with site outermost, parallel_for
|
||||
for(int b=0;b<Nblock;b++){
|
||||
AP[b] = zero;
|
||||
for(int bp=0;bp<Nblock;bp++) {
|
||||
AP[b] += (m(bp,b))*X[bp];
|
||||
}
|
||||
}
|
||||
}
|
||||
double normv(const std::vector<Field> &P){
|
||||
double nn = 0.0;
|
||||
for(int b=0;b<Nblock;b++) {
|
||||
nn+=norm2(P[b]);
|
||||
}
|
||||
return nn;
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
// BlockCGrQvec implementation:
|
||||
//--------------------------
|
||||
// X is guess/Solution
|
||||
// B is RHS
|
||||
// Solve A X_i = B_i ; i refers to Nblock index
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field> &B, std::vector<Field> &X)
|
||||
{
|
||||
Nblock = B.size();
|
||||
assert(Nblock == X.size());
|
||||
|
||||
std::cout<<GridLogMessage<<" Block Conjugate Gradient Vec rQ : Nblock "<<Nblock<<std::endl;
|
||||
|
||||
for(int b=0;b<Nblock;b++){
|
||||
X[b].checkerboard = B[b].checkerboard;
|
||||
conformable(X[b], B[b]);
|
||||
conformable(X[b], X[0]);
|
||||
}
|
||||
|
||||
Field Fake(B[0]);
|
||||
|
||||
std::vector<Field> tmp(Nblock,Fake);
|
||||
std::vector<Field> Q(Nblock,Fake);
|
||||
std::vector<Field> D(Nblock,Fake);
|
||||
std::vector<Field> Z(Nblock,Fake);
|
||||
std::vector<Field> AD(Nblock,Fake);
|
||||
|
||||
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);
|
||||
|
||||
RealD sssum=0;
|
||||
for(int b=0;b<Nblock;b++){ ssq[b] = norm2(B[b]);}
|
||||
for(int b=0;b<Nblock;b++) sssum+=ssq[b];
|
||||
|
||||
for(int b=0;b<Nblock;b++){ residuals[b] = norm2(B[b]);}
|
||||
for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
|
||||
|
||||
for(int b=0;b<Nblock;b++){ residuals[b] = norm2(X[b]);}
|
||||
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<<"BlockCGrQvec algorithm initialisation " <<std::endl;
|
||||
|
||||
//1. QC = R = B-AX, D = Q ; QC => Thin QR factorisation (google it)
|
||||
for(int b=0;b<Nblock;b++) {
|
||||
Linop.HermOp(X[b], AD[b]);
|
||||
tmp[b] = B[b] - AD[b];
|
||||
}
|
||||
|
||||
ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);
|
||||
|
||||
for(int b=0;b<Nblock;b++) D[b]=Q[b];
|
||||
|
||||
std::cout << GridLogMessage<<"BlockCGrQ vec 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();
|
||||
for(int b=0;b<Nblock;b++) Linop.HermOp(D[b], Z[b]);
|
||||
MatrixTimer.Stop();
|
||||
|
||||
//4. M = [D^dag Z]^{-1}
|
||||
sliceInnerTimer.Start();
|
||||
InnerProductMatrix(m_DZ,D,Z);
|
||||
sliceInnerTimer.Stop();
|
||||
m_M = m_DZ.inverse();
|
||||
|
||||
//5. X = X + D MC
|
||||
m_tmp = m_M * m_C;
|
||||
sliceMaddTimer.Start();
|
||||
MaddMatrix(X,m_tmp, D,X);
|
||||
sliceMaddTimer.Stop();
|
||||
|
||||
//6. QS = Q - ZM
|
||||
sliceMaddTimer.Start();
|
||||
MaddMatrix(tmp,m_M,Z,Q,-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();
|
||||
MaddMatrix(D,m_tmp,D,Q);
|
||||
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 << "\t Block Iteration "<<k<<" ave resid "<< sqrt(rrsum/sssum) << " max "<< sqrt(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;
|
||||
|
||||
for(int b=0;b<Nblock;b++) Linop.HermOp(X[b], AD[b]);
|
||||
for(int b=0;b<Nblock;b++) AD[b] = AD[b]-B[b];
|
||||
std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(normv(AD)/normv(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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -133,7 +133,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
|
||||
LinalgTimer.Stop();
|
||||
|
||||
std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k
|
||||
<< " residual " << cp << " target " << rsq << std::endl;
|
||||
<< " residual^2 " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
|
||||
|
||||
// Stopping condition
|
||||
if (cp <= rsq) {
|
||||
@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction<Field> {
|
||||
std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
|
||||
std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
|
||||
|
||||
std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
|
||||
std::cout << GridLogPerformance << "\tElapsed " << SolverTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tLinalg " << LinalgTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tInner " << InnerTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tAxpyNorm " << AxpyNormTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<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 << "\tLinalg " << LinalgTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tInner " << InnerTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tAxpyNorm " << AxpyNormTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
|
||||
|
||||
if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
|
||||
|
||||
|
||||
@ -38,6 +38,7 @@ int main (int argc, char ** argv)
|
||||
typedef typename DomainWallFermionR::ComplexField ComplexField;
|
||||
typename DomainWallFermionR::ImplParams params;
|
||||
|
||||
double stp=1.0e-5;
|
||||
const int Ls=4;
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
@ -197,7 +198,7 @@ int main (int argc, char ** argv)
|
||||
|
||||
MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
|
||||
MdagMLinearOperator<DomainWallFermionR,FermionField> HermOpCk(Dchk);
|
||||
ConjugateGradient<FermionField> CG((1.0e-2),10000);
|
||||
ConjugateGradient<FermionField> CG((stp),10000);
|
||||
s_res = zero;
|
||||
CG(HermOp,s_src,s_res);
|
||||
|
||||
@ -227,5 +228,11 @@ int main (int argc, char ** argv)
|
||||
std::cout << GridLogMessage<<" resid["<<n<<"] "<< norm2(tmp)/norm2(src[n])<<std::endl;
|
||||
}
|
||||
|
||||
for(int s=0;s<nrhs;s++) result[s]=zero;
|
||||
int blockDim = 0;//not used for BlockCGVec
|
||||
BlockConjugateGradient<FermionField> BCGV (BlockCGVec,blockDim,stp,10000);
|
||||
BCGV.PrintInterval=10;
|
||||
BCGV(HermOpCk,src,result);
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
220
tests/solver/Test_mobius_bcg.cc
Normal file
220
tests/solver/Test_mobius_bcg.cc
Normal file
@ -0,0 +1,220 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_mrhs_cg.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
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 */
|
||||
#include <Grid/Grid.h>
|
||||
#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
typedef typename MobiusFermionR::FermionField FermionField;
|
||||
typedef typename MobiusFermionR::ComplexField ComplexField;
|
||||
typename MobiusFermionR::ImplParams params;
|
||||
|
||||
const int Ls=12;
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
std::vector<int> mpi_split (mpi_layout.size(),1);
|
||||
std::vector<int> split_coor (mpi_layout.size(),1);
|
||||
std::vector<int> split_dim (mpi_layout.size(),1);
|
||||
|
||||
std::vector<ComplexD> boundary_phases(Nd,1.);
|
||||
boundary_phases[Nd-1]=-1.;
|
||||
params.boundary_phases = boundary_phases;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Split into 1^4 mpi communicators
|
||||
/////////////////////////////////////////////
|
||||
|
||||
for(int i=0;i<argc;i++){
|
||||
if(std::string(argv[i]) == "--split"){
|
||||
for(int k=0;k<mpi_layout.size();k++){
|
||||
std::stringstream ss;
|
||||
ss << argv[i+1+k];
|
||||
ss >> mpi_split[k];
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double stp = 1.e-8;
|
||||
int nrhs = 1;
|
||||
int me;
|
||||
for(int i=0;i<mpi_layout.size();i++){
|
||||
// split_dim[i] = (mpi_layout[i]/mpi_split[i]);
|
||||
nrhs *= (mpi_layout[i]/mpi_split[i]);
|
||||
// split_coor[i] = FGrid._processor_coor[i]/mpi_split[i];
|
||||
}
|
||||
std::cout << GridLogMessage << "Creating split grids " <<std::endl;
|
||||
GridCartesian * SGrid = new GridCartesian(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
mpi_split,
|
||||
*UGrid,me);
|
||||
std::cout << GridLogMessage <<"Creating split ferm grids " <<std::endl;
|
||||
|
||||
GridCartesian * SFGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,SGrid);
|
||||
std::cout << GridLogMessage <<"Creating split rb grids " <<std::endl;
|
||||
GridRedBlackCartesian * SrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(SGrid);
|
||||
std::cout << GridLogMessage <<"Creating split ferm rb grids " <<std::endl;
|
||||
GridRedBlackCartesian * SFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,SGrid);
|
||||
std::cout << GridLogMessage << "Made the grids"<<std::endl;
|
||||
///////////////////////////////////////////////
|
||||
// Set up the problem as a 4d spreadout job
|
||||
///////////////////////////////////////////////
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
|
||||
std::vector<FermionField> src(nrhs,FGrid);
|
||||
std::vector<FermionField> src_chk(nrhs,FGrid);
|
||||
std::vector<FermionField> result(nrhs,FGrid);
|
||||
FermionField tmp(FGrid);
|
||||
std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
|
||||
|
||||
for(int s=0;s<nrhs;s++) result[s]=zero;
|
||||
GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds);
|
||||
for(int s=0;s<nrhs;s++) {
|
||||
random(pRNG5,src[s]);
|
||||
std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
if(0) {
|
||||
FieldMetaData header;
|
||||
std::string file("./lat.in");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
std::cout << GridLogMessage << " "<<file<<" successfully read" <<std::endl;
|
||||
} else {
|
||||
GridParallelRNG pRNG(UGrid );
|
||||
std::cout << GridLogMessage << "Intialising 4D RNG "<<std::endl;
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
std::cout << GridLogMessage << "Intialised 4D RNG "<<std::endl;
|
||||
SU3::HotConfiguration(pRNG,Umu);
|
||||
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
|
||||
std::cout << " Site zero "<< Umu._odata[0] <<std::endl;
|
||||
}
|
||||
|
||||
/////////////////
|
||||
// MPI only sends
|
||||
/////////////////
|
||||
LatticeGaugeField s_Umu(SGrid);
|
||||
FermionField s_src(SFGrid);
|
||||
FermionField s_tmp(SFGrid);
|
||||
FermionField s_res(SFGrid);
|
||||
|
||||
std::cout << GridLogMessage << "Made the split grid fields"<<std::endl;
|
||||
///////////////////////////////////////////////////////////////
|
||||
// split the source out using MPI instead of I/O
|
||||
///////////////////////////////////////////////////////////////
|
||||
Grid_split (Umu,s_Umu);
|
||||
Grid_split (src,s_src);
|
||||
std::cout << GridLogMessage << " split rank " <<me << " s_src "<<norm2(s_src)<<std::endl;
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Set up N-solvers as trivially parallel
|
||||
///////////////////////////////////////////////////////////////
|
||||
std::cout << GridLogMessage << " Building the solvers"<<std::endl;
|
||||
// RealD mass=0.00107;
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.8;
|
||||
RealD mobius_factor=32./12.;
|
||||
RealD mobius_b=0.5*(mobius_factor+1.);
|
||||
RealD mobius_c=0.5*(mobius_factor-1.);
|
||||
MobiusFermionR Dchk(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,mobius_b,mobius_c,params);
|
||||
MobiusFermionR Ddwf(s_Umu,*SFGrid,*SFrbGrid,*SGrid,*SrbGrid,mass,M5,mobius_b,mobius_c,params);
|
||||
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
|
||||
MdagMLinearOperator<MobiusFermionR,FermionField> HermOp(Ddwf);
|
||||
MdagMLinearOperator<MobiusFermionR,FermionField> HermOpCk(Dchk);
|
||||
ConjugateGradient<FermionField> CG((stp),100000);
|
||||
s_res = zero;
|
||||
|
||||
CG(HermOp,s_src,s_res);
|
||||
|
||||
std::cout << GridLogMessage << " split residual norm "<<norm2(s_res)<<std::endl;
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Report how long they all took
|
||||
/////////////////////////////////////////////////////////////
|
||||
std::vector<uint32_t> iterations(nrhs,0);
|
||||
iterations[me] = CG.IterationsToComplete;
|
||||
|
||||
for(int n=0;n<nrhs;n++){
|
||||
UGrid->GlobalSum(iterations[n]);
|
||||
std::cout << GridLogMessage<<" Rank "<<n<<" "<< iterations[n]<<" CG iterations"<<std::endl;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Gather and residual check on the results
|
||||
/////////////////////////////////////////////////////////////
|
||||
std::cout << GridLogMessage<< "Unsplitting the result"<<std::endl;
|
||||
Grid_unsplit(result,s_res);
|
||||
|
||||
|
||||
std::cout << GridLogMessage<< "Checking the residuals"<<std::endl;
|
||||
for(int n=0;n<nrhs;n++){
|
||||
std::cout << GridLogMessage<< " res["<<n<<"] norm "<<norm2(result[n])<<std::endl;
|
||||
HermOpCk.HermOp(result[n],tmp); tmp = tmp - src[n];
|
||||
std::cout << GridLogMessage<<" resid["<<n<<"] "<< std::sqrt(norm2(tmp)/norm2(src[n]))<<std::endl;
|
||||
}
|
||||
|
||||
|
||||
for(int s=0;s<nrhs;s++){
|
||||
result[s]=zero;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Try block CG
|
||||
/////////////////////////////////////////////////////////////
|
||||
int blockDim = 0;//not used for BlockCGVec
|
||||
BlockConjugateGradient<FermionField> BCGV (BlockCGrQVec,blockDim,stp,100000);
|
||||
{
|
||||
BCGV(HermOpCk,src,result);
|
||||
}
|
||||
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
144
tests/solver/Test_mobius_bcg_nosplit.cc
Normal file
144
tests/solver/Test_mobius_bcg_nosplit.cc
Normal file
@ -0,0 +1,144 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_mrhs_cg.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
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 */
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
typedef typename DomainWallFermionR::FermionField FermionField;
|
||||
typedef typename DomainWallFermionR::ComplexField ComplexField;
|
||||
typename DomainWallFermionR::ImplParams params;
|
||||
|
||||
const int Ls=16;
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
|
||||
std::vector<ComplexD> boundary_phases(Nd,1.);
|
||||
boundary_phases[Nd-1]=-1.;
|
||||
params.boundary_phases = boundary_phases;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
double stp = 1.e-8;
|
||||
int nrhs = 2;
|
||||
|
||||
///////////////////////////////////////////////
|
||||
// Set up the problem as a 4d spreadout job
|
||||
///////////////////////////////////////////////
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
|
||||
std::vector<FermionField> src(nrhs,FGrid);
|
||||
std::vector<FermionField> src_chk(nrhs,FGrid);
|
||||
std::vector<FermionField> result(nrhs,FGrid);
|
||||
FermionField tmp(FGrid);
|
||||
std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
|
||||
|
||||
for(int s=0;s<nrhs;s++) result[s]=zero;
|
||||
GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds);
|
||||
for(int s=0;s<nrhs;s++) {
|
||||
random(pRNG5,src[s]);
|
||||
std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
int conf = 0;
|
||||
if(conf==0) {
|
||||
FieldMetaData header;
|
||||
std::string file("./lat.in");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
std::cout << GridLogMessage << " Config "<<file<<" successfully read" <<std::endl;
|
||||
} else if (conf==1){
|
||||
GridParallelRNG pRNG(UGrid );
|
||||
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
SU3::HotConfiguration(pRNG,Umu);
|
||||
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
|
||||
} else {
|
||||
SU3::ColdConfiguration(Umu);
|
||||
std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Set up N-solvers as trivially parallel
|
||||
///////////////////////////////////////////////////////////////
|
||||
std::cout << GridLogMessage << " Building the solvers"<<std::endl;
|
||||
RealD mass=0.01;
|
||||
RealD M5=1.8;
|
||||
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,params);
|
||||
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
|
||||
MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
|
||||
ConjugateGradient<FermionField> CG((stp),100000);
|
||||
|
||||
for(int rhs=0;rhs<1;rhs++){
|
||||
result[rhs] = zero;
|
||||
CG(HermOp,src[rhs],result[rhs]);
|
||||
}
|
||||
|
||||
for(int rhs=0;rhs<1;rhs++){
|
||||
std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Try block CG
|
||||
/////////////////////////////////////////////////////////////
|
||||
int blockDim = 0;//not used for BlockCGVec
|
||||
for(int s=0;s<nrhs;s++){
|
||||
result[s]=zero;
|
||||
}
|
||||
BlockConjugateGradient<FermionField> BCGV (BlockCGrQVec,blockDim,stp,100000);
|
||||
{
|
||||
BCGV(HermOp,src,result);
|
||||
}
|
||||
|
||||
for(int rhs=0;rhs<nrhs;rhs++){
|
||||
std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
148
tests/solver/Test_mobius_bcg_phys_nosplit.cc
Normal file
148
tests/solver/Test_mobius_bcg_phys_nosplit.cc
Normal file
@ -0,0 +1,148 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_mrhs_cg.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
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 */
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
typedef typename DomainWallFermionR::FermionField FermionField;
|
||||
typedef typename DomainWallFermionR::ComplexField ComplexField;
|
||||
typename DomainWallFermionR::ImplParams params;
|
||||
|
||||
const int Ls=16;
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
|
||||
std::vector<ComplexD> boundary_phases(Nd,1.);
|
||||
boundary_phases[Nd-1]=-1.;
|
||||
params.boundary_phases = boundary_phases;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
double stp = 1.e-8;
|
||||
int nrhs = 2;
|
||||
|
||||
///////////////////////////////////////////////
|
||||
// Set up the problem as a 4d spreadout job
|
||||
///////////////////////////////////////////////
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
|
||||
std::vector<FermionField> src4(nrhs,UGrid);
|
||||
std::vector<FermionField> src(nrhs,FGrid);
|
||||
std::vector<FermionField> src_chk(nrhs,FGrid);
|
||||
std::vector<FermionField> result(nrhs,FGrid);
|
||||
FermionField tmp(FGrid);
|
||||
std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
|
||||
|
||||
for(int s=0;s<nrhs;s++) result[s]=zero;
|
||||
GridParallelRNG pRNG4(UGrid); pRNG4.SeedFixedIntegers(seeds);
|
||||
for(int s=0;s<nrhs;s++) {
|
||||
random(pRNG4,src4[s]);
|
||||
std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
int conf = 0;
|
||||
if(conf==0) {
|
||||
FieldMetaData header;
|
||||
std::string file("./lat.in");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
std::cout << GridLogMessage << " Config "<<file<<" successfully read" <<std::endl;
|
||||
} else if (conf==1){
|
||||
GridParallelRNG pRNG(UGrid );
|
||||
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
SU3::HotConfiguration(pRNG,Umu);
|
||||
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
|
||||
} else {
|
||||
SU3::ColdConfiguration(Umu);
|
||||
std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Set up N-solvers as trivially parallel
|
||||
///////////////////////////////////////////////////////////////
|
||||
std::cout << GridLogMessage << " Building the solvers"<<std::endl;
|
||||
RealD mass=0.01;
|
||||
RealD M5=1.8;
|
||||
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,params);
|
||||
for(int s=0;s<nrhs;s++) {
|
||||
Ddwf.ImportPhysicalFermionSource(src4[s],src[s]);
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
|
||||
MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
|
||||
ConjugateGradient<FermionField> CG((stp),100000);
|
||||
|
||||
for(int rhs=0;rhs<1;rhs++){
|
||||
result[rhs] = zero;
|
||||
// CG(HermOp,src[rhs],result[rhs]);
|
||||
}
|
||||
|
||||
for(int rhs=0;rhs<1;rhs++){
|
||||
std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Try block CG
|
||||
/////////////////////////////////////////////////////////////
|
||||
int blockDim = 0;//not used for BlockCGVec
|
||||
for(int s=0;s<nrhs;s++){
|
||||
result[s]=zero;
|
||||
}
|
||||
BlockConjugateGradient<FermionField> BCGV (BlockCGrQVec,blockDim,stp,100000);
|
||||
{
|
||||
BCGV(HermOp,src,result);
|
||||
}
|
||||
|
||||
for(int rhs=0;rhs<nrhs;rhs++){
|
||||
std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
144
tests/solver/Test_mobius_bcg_prec_nosplit.cc
Normal file
144
tests/solver/Test_mobius_bcg_prec_nosplit.cc
Normal file
@ -0,0 +1,144 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_mrhs_cg.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
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 */
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
typedef typename DomainWallFermionR::FermionField FermionField;
|
||||
typedef typename DomainWallFermionR::ComplexField ComplexField;
|
||||
typename DomainWallFermionR::ImplParams params;
|
||||
|
||||
const int Ls=16;
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||
|
||||
std::vector<ComplexD> boundary_phases(Nd,1.);
|
||||
boundary_phases[Nd-1]=-1.;
|
||||
params.boundary_phases = boundary_phases;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
double stp = 1.e-8;
|
||||
int nrhs = 4;
|
||||
|
||||
///////////////////////////////////////////////
|
||||
// Set up the problem as a 4d spreadout job
|
||||
///////////////////////////////////////////////
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
|
||||
std::vector<FermionField> src(nrhs,FGrid);
|
||||
std::vector<FermionField> src_chk(nrhs,FGrid);
|
||||
std::vector<FermionField> result(nrhs,FGrid);
|
||||
FermionField tmp(FGrid);
|
||||
std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
|
||||
|
||||
for(int s=0;s<nrhs;s++) result[s]=zero;
|
||||
GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds);
|
||||
for(int s=0;s<nrhs;s++) {
|
||||
random(pRNG5,src[s]);
|
||||
std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
int conf = 0;
|
||||
if(conf==0) {
|
||||
FieldMetaData header;
|
||||
std::string file("./lat.in");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
std::cout << GridLogMessage << " Config "<<file<<" successfully read" <<std::endl;
|
||||
} else if (conf==1){
|
||||
GridParallelRNG pRNG(UGrid );
|
||||
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
SU3::HotConfiguration(pRNG,Umu);
|
||||
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
|
||||
} else {
|
||||
SU3::ColdConfiguration(Umu);
|
||||
std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Set up N-solvers as trivially parallel
|
||||
///////////////////////////////////////////////////////////////
|
||||
std::cout << GridLogMessage << " Building the solvers"<<std::endl;
|
||||
RealD mass=0.01;
|
||||
RealD M5=1.8;
|
||||
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,params);
|
||||
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
|
||||
MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
|
||||
ConjugateGradient<FermionField> CG((stp),100000);
|
||||
|
||||
for(int rhs=0;rhs<1;rhs++){
|
||||
result[rhs] = zero;
|
||||
CG(HermOp,src[rhs],result[rhs]);
|
||||
}
|
||||
|
||||
for(int rhs=0;rhs<1;rhs++){
|
||||
std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Try block CG
|
||||
/////////////////////////////////////////////////////////////
|
||||
int blockDim = 0;//not used for BlockCGVec
|
||||
for(int s=0;s<nrhs;s++){
|
||||
result[s]=zero;
|
||||
}
|
||||
BlockConjugateGradient<FermionField> BCGV (BlockCGrQVec,blockDim,stp,100000);
|
||||
{
|
||||
BCGV(HermOp,src,result);
|
||||
}
|
||||
|
||||
for(int rhs=0;rhs<nrhs;rhs++){
|
||||
std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -67,7 +67,22 @@ int main (int argc, char ** argv)
|
||||
GridParallelRNG pRNG(UGrid ); pRNG.SeedFixedIntegers(seeds);
|
||||
GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds);
|
||||
|
||||
FermionField src(FGrid); random(pRNG5,src);
|
||||
FermionField src(FGrid);
|
||||
FermionField tt(FGrid);
|
||||
#if 1
|
||||
random(pRNG5,src);
|
||||
#else
|
||||
src=zero;
|
||||
ComplexField coor(FGrid);
|
||||
LatticeCoordinate(coor,0);
|
||||
for(int ss=0;ss<FGrid->oSites();ss++){
|
||||
src._odata[ss]()()(0)=coor._odata[ss]()()();
|
||||
}
|
||||
LatticeCoordinate(coor,1);
|
||||
for(int ss=0;ss<FGrid->oSites();ss++){
|
||||
src._odata[ss]()()(0)+=coor._odata[ss]()()();
|
||||
}
|
||||
#endif
|
||||
FermionField src_o(FrbGrid); pickCheckerboard(Odd,src_o,src);
|
||||
FermionField result_o(FrbGrid); result_o=zero;
|
||||
RealD nrm = norm2(src);
|
||||
@ -89,7 +104,8 @@ int main (int argc, char ** argv)
|
||||
ConjugateGradient<FermionField> CG(1.0e-8,10000);
|
||||
int blockDim = 0;
|
||||
BlockConjugateGradient<FermionField> BCGrQ(BlockCGrQ,blockDim,1.0e-8,10000);
|
||||
BlockConjugateGradient<FermionField> BCG (BlockCG,blockDim,1.0e-8,10000);
|
||||
BlockConjugateGradient<FermionField> BCG (BlockCGrQ,blockDim,1.0e-8,10000);
|
||||
BlockConjugateGradient<FermionField> BCGv (BlockCGrQVec,blockDim,1.0e-8,10000);
|
||||
BlockConjugateGradient<FermionField> mCG (CGmultiRHS,blockDim,1.0e-8,10000);
|
||||
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
@ -158,7 +174,7 @@ int main (int argc, char ** argv)
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling Block CG for "<<Ls <<" right hand sides" <<std::endl;
|
||||
std::cout << GridLogMessage << " Calling Block CGrQ for "<<Ls <<" right hand sides" <<std::endl;
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
Ds.ZeroCounters();
|
||||
result_o=zero;
|
||||
@ -176,6 +192,49 @@ int main (int argc, char ** argv)
|
||||
Ds.Report();
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling Block CG for "<<Ls <<" right hand sides" <<std::endl;
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
Ds.ZeroCounters();
|
||||
result_o=zero;
|
||||
{
|
||||
double t1=usecond();
|
||||
BCG(HermOp,src_o,result_o);
|
||||
double t2=usecond();
|
||||
double ncall=BCGrQ.IterationsToComplete*Ls;
|
||||
double flops = deodoe_flops * ncall;
|
||||
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
|
||||
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
|
||||
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
|
||||
HermOp.Report();
|
||||
}
|
||||
Ds.Report();
|
||||
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
std::cout << GridLogMessage << " Calling BCGvec "<<std::endl;
|
||||
std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
|
||||
std::vector<FermionField> src_v (Ls,UrbGrid);
|
||||
std::vector<FermionField> result_v(Ls,UrbGrid);
|
||||
for(int s=0;s<Ls;s++) result_v[s] = zero;
|
||||
for(int s=0;s<Ls;s++) {
|
||||
FermionField src4(UGrid);
|
||||
ExtractSlice(src4,src,s,0);
|
||||
pickCheckerboard(Odd,src_v[s],src4);
|
||||
}
|
||||
|
||||
{
|
||||
double t1=usecond();
|
||||
BCGv(HermOp4d,src_v,result_v);
|
||||
double t2=usecond();
|
||||
double ncall=BCGv.IterationsToComplete*Ls;
|
||||
double flops = deodoe_flops * ncall;
|
||||
std::cout<<GridLogMessage << "usec = "<< (t2-t1)<<std::endl;
|
||||
std::cout<<GridLogMessage << "flops = "<< flops<<std::endl;
|
||||
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t2-t1)<<std::endl;
|
||||
// HermOp4d.Report();
|
||||
}
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user