From b4a14224c3d1b590cf3ee8b1c42fcf622be071e9 Mon Sep 17 00:00:00 2001 From: Eloy Romero Date: Tue, 16 Feb 2021 10:35:28 -0500 Subject: [PATCH 1/6] Reimplement projectVector to do oblique projection. --- include/eigsolve_mugiq.h | 16 +++------ lib/eigsolve_mugiq.cpp | 72 +++++++++++++++++++++++++++++++++++----- 2 files changed, 68 insertions(+), 20 deletions(-) diff --git a/include/eigsolve_mugiq.h b/include/eigsolve_mugiq.h index 30f58ab..8f223fe 100644 --- a/include/eigsolve_mugiq.h +++ b/include/eigsolve_mugiq.h @@ -75,6 +75,7 @@ class Eigsolve_Mugiq { const Dirac *dirac; DiracMatrix *mat; // The Dirac operator whose eigenpairs we are computing + DiracMatrix *matFine; // The Dirac operator whose eigenpairs we are computing //- This switch is required so that the dirac object is NOT //- deleted when NOT created within Eigsolve. @@ -151,7 +152,7 @@ class Eigsolve_Mugiq { */ void computeEvals(); - /** @brief Perform the projection: out = \sum_i v_i dot(v*_i,in) + /** @brief Perform the projection: out = \sum_i evecs_i * dot(evecs_i*,\gamma_5 * fine_op * in) / eval_i */ void projectVector(ColorSpinorField &out, ColorSpinorField &in); @@ -159,22 +160,15 @@ class Eigsolve_Mugiq { */ void printEvals(); - /** @brief Accessor to get the eigenvectors outside of the class + /** @brief Accessor to get approx. right singular vector of the given fine/coarse operator */ std::vector &getEvecs(){ return eVecs;} - /** @brief Accessor to get the Quda eigenvalues outside of the class - */ - std::vector* getEvalsQuda(){ return eVals_quda;} - - /** @brief Accessor to get the Eigsolve_Mugiq eigenvalues outside of the class + /** @brief Accessor to get the approx. eigenvalues of \gamma_5 * the given fine/coarse operator; + * that is, evals[i] = dot(evecs[i], \gamma_5 * fine/coarse operator * evecs[i]) */ std::vector *getEvals(){ return eVals;} - /** @brief Accessor to get the Eigsolve_Mugiq singular values outside of the class - */ - std::vector *getEvalsSigma(){ return eVals_sigma;} - /** @brief Accessor to get the residual of the computed eigenvalues */ std::vector* getEvalsRes(){ return evals_res;} diff --git a/lib/eigsolve_mugiq.cpp b/lib/eigsolve_mugiq.cpp index 73d3a14..31ee0d2 100644 --- a/lib/eigsolve_mugiq.cpp +++ b/lib/eigsolve_mugiq.cpp @@ -15,6 +15,7 @@ Eigsolve_Mugiq::Eigsolve_Mugiq(MugiqEigParam *eigParams_, eigProfile(eigProfile_), dirac(nullptr), mat(nullptr), + matFine(nullptr), diracCreated(MUGIQ_BOOL_FALSE), eVals_quda(nullptr), eVals(nullptr), @@ -63,6 +64,7 @@ Eigsolve_Mugiq::Eigsolve_Mugiq(MugiqEigParam *eigParams_, eigProfile(eigProfile_), dirac(nullptr), mat(nullptr), + matFine(nullptr), diracCreated(MUGIQ_BOOL_FALSE), eVals_quda(nullptr), eVals(nullptr), @@ -94,6 +96,7 @@ Eigsolve_Mugiq::~Eigsolve_Mugiq(){ delete eVals_sigma; if(mat) delete mat; + if(matFine) delete mat; mat = nullptr; if(useMGenv){ @@ -204,6 +207,8 @@ void Eigsolve_Mugiq::createNewDiracMatrix(){ else if (eigParams->diracType == MUGIQ_EIG_OPERATOR_MdagM) mat = new DiracMdagM(*dirac); else if (eigParams->diracType == MUGIQ_EIG_OPERATOR_MMdag) mat = new DiracMMdag(*dirac); else errorQuda("%s: Unsupported Dirac operator type\n", __func__); + + matFine = new DiracM(*mg_env->mg_solver->d); } @@ -283,6 +288,13 @@ void Eigsolve_Mugiq::computeEvecs(){ EigenSolver *eigSolve = EigenSolver::create(eigParams->QudaEigParams, *mat, *eigProfile); (*eigSolve)(eVecs, *eVals_quda); + // Get the right singular vectors if the solver returns the left singular vectors + if(eigParams->diracType == MUGIQ_EIG_OPERATOR_Mdag || eigParams->diracType == MUGIQ_EIG_OPERATOR_MMdag){ + for(int i=0; inEv; i++){ + gamma5(*eVecs[i], *eVecs[i]); + } + } + delete eigSolve; } @@ -296,13 +308,15 @@ void Eigsolve_Mugiq::computeEvals(){ std::vector &r = *evals_res; double kappa = invParams->kappa; + DiracMatrix *mat_direct = new DiracM(*dirac); for(int i=0; inEv; i++){ - (*mat)(*w,*eVecs[i]); //- w = M*v_i + (*mat_direct)(*w,*eVecs[i]); //- w = M*v_i if(invParams->mass_normalization == QUDA_MASS_NORMALIZATION) blas::ax(0.25/(kappa*kappa), *w); + gamma5(*w, *w); lambda[i] = blas::cDotProduct(*eVecs[i], *w) / sqrt(blas::norm2(*eVecs[i])); // lambda_i = (v_i^dag M v_i) / ||v_i|| Complex Cm1(-1.0, 0.0); - blas::caxpby(lambda[i], *eVecs[i], Cm1, *w); // w = lambda_i*v_i - A*v_i + blas::caxpby(lambda[i], *eVecs[i], Cm1, *w); // w = lambda_i*v_i - \gamma_5*A*v_i r[i] = sqrt(blas::norm2(*w)); // r = ||w|| } @@ -311,6 +325,7 @@ void Eigsolve_Mugiq::computeEvals(){ for(int i=0; inEv; i++) sigma[i] = sqrt(lambda[i].real()); } + delete mat_direct; delete w; } @@ -335,14 +350,53 @@ void Eigsolve_Mugiq::printEvals(){ } /** - * Perform the projection: out = \sum_i v_i dot(v*_i,in) + * Perform the projection: out = \sum_i evecs_i * dot(evecs_i*,\gamma_5 * fine_op * in) / eval_i */ void Eigsolve_Mugiq::projectVector(ColorSpinorField &out, ColorSpinorField &in){ - - blas::zero(out); - for(int i=0; inEv; i++){ - Complex dp = blas::cDotProduct(*eVecs[i], in); // dp = dot(v_i*,in) - blas::caxpy(dp,*eVecs[i],out); // out = dp*v_i + out + // Min = gamma_5 * matFine * in + ColorSpinorParam csParam(*tmpCSF[0]); + ColorSpinorField *Min = ColorSpinorField::Create(csParam); + (*matFine)(*Min,in); + gamma5(*Min, *Min); + + // Transfer Min to the coarsest level + ColorSpinorField *Min_coarse = nullptr, *Mout_coarse = nullptr; + if (computeCoarse && mg_env->nCoarseLevels > 0) { + blas::zero(*tmpCSF[1]); + mg_env->transfer[0]->R(*tmpCSF[1], *Min); + for(int lev=1; levnCoarseLevels; lev++){ + blas::zero(*tmpCSF[lev+1]); + if(!mg_env->transfer[lev]) errorQuda("%s: Transfer operator for level %d does not exist!\n", __func__, lev); + mg_env->transfer[lev]->R(*tmpCSF[lev+1], *tmpCSF[lev]); + } + Mout_coarse = Min_coarse = tmpCSF[mg_env->nCoarseLevels]; + } else { + Min_coarse = Min; + Mout_coarse = &out; } - + + // s[i] = dot(eVecs[i], Min_coarse) + std::vector s(eigParams->nEv); + std::vector Minc_(1, Min_coarse); + blas::cDotProduct(s.data(), eVecs, Minc_); + + // Mout_coarse = Sum_i s[i] / evals[i] * eVecs[i] + blas::zero(*Mout_coarse); + for (int i = 0; i < eigParams->nEv; i++) { + s[i] /= (*eVals)[i]; + } + std::vector Moutc_(1, Mout_coarse); + blas::caxpy(s.data(), eVecs, Moutc_); + + // Transfer Mout_coarse to the fine level + if (computeCoarse && mg_env->nCoarseLevels > 0) { + for(int lev=mg_env->nCoarseLevels; lev>1; lev--){ + blas::zero(*tmpCSF[lev-1]); + mg_env->transfer[lev-1]->P(*tmpCSF[lev-1], *tmpCSF[lev]); + } + blas::zero(out); + mg_env->transfer[0]->P(out, *tmpCSF[1]); + } + + delete Min; } From 4d9ee05bc5707ae13172831bc3180d80d2bca715 Mon Sep 17 00:00:00 2001 From: ckallidonis Date: Wed, 17 Feb 2021 15:04:51 -0500 Subject: [PATCH 2/6] Fixed bug on sigma-eigenvalue computation --- lib/eigsolve_mugiq.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/lib/eigsolve_mugiq.cpp b/lib/eigsolve_mugiq.cpp index 31ee0d2..5b051d6 100644 --- a/lib/eigsolve_mugiq.cpp +++ b/lib/eigsolve_mugiq.cpp @@ -320,10 +320,8 @@ void Eigsolve_Mugiq::computeEvals(){ r[i] = sqrt(blas::norm2(*w)); // r = ||w|| } - if(eigParams->diracType == MUGIQ_EIG_OPERATOR_MdagM || eigParams->diracType == MUGIQ_EIG_OPERATOR_MMdag){ - std::vector &sigma = *eVals_sigma; - for(int i=0; inEv; i++) sigma[i] = sqrt(lambda[i].real()); - } + std::vector &sigma = *eVals_sigma; + for(int i=0; inEv; i++) sigma[i] = lambda[i].real(); delete mat_direct; delete w; From ca56c933df47ccad8f2dfbf91da85cd8bf105aaa Mon Sep 17 00:00:00 2001 From: ckallidonis Date: Wed, 17 Feb 2021 15:09:55 -0500 Subject: [PATCH 3/6] Fixed bug when creating the fine Dirac operator in Eigsolve class --- lib/eigsolve_mugiq.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/eigsolve_mugiq.cpp b/lib/eigsolve_mugiq.cpp index 5b051d6..55f3b1c 100644 --- a/lib/eigsolve_mugiq.cpp +++ b/lib/eigsolve_mugiq.cpp @@ -208,7 +208,7 @@ void Eigsolve_Mugiq::createNewDiracMatrix(){ else if (eigParams->diracType == MUGIQ_EIG_OPERATOR_MMdag) mat = new DiracMMdag(*dirac); else errorQuda("%s: Unsupported Dirac operator type\n", __func__); - matFine = new DiracM(*mg_env->mg_solver->d); + matFine = new DiracM(*dirac); } From 91bedec7666cf349fdb11abef3acf636d70147de Mon Sep 17 00:00:00 2001 From: ckallidonis Date: Wed, 17 Feb 2021 15:20:02 -0500 Subject: [PATCH 4/6] Removed unneccesary get-er functions from Eigsolve class --- include/eigsolve_mugiq.h | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/include/eigsolve_mugiq.h b/include/eigsolve_mugiq.h index 8f223fe..9856b1b 100644 --- a/include/eigsolve_mugiq.h +++ b/include/eigsolve_mugiq.h @@ -169,26 +169,10 @@ class Eigsolve_Mugiq { */ std::vector *getEvals(){ return eVals;} - /** @brief Accessor to get the residual of the computed eigenvalues - */ - std::vector* getEvalsRes(){ return evals_res;} - /** @brief Accessor to get the Multigrid environment structure */ MG_Mugiq* getMGEnv(){ return mg_env;} - /** @brief Accessor to get the Mugiq eigsolve parameter structure - */ - MugiqEigParam* getMugiqEigParams(){ return eigParams;} - - /** @brief Accessor to get the Quda eigsolve parameter structure - */ - QudaEigParam* getQudaEigParams(){ return eigParams->QudaEigParams;} - - /** @brief Accessor to get the invert parameter structure - */ - QudaInvertParam* getInvParams(){ return invParams;} - }; // class Eigsolve_Mugiq From ed47c059eb2e3e6023f6b134decd8ba5fbd7f7a1 Mon Sep 17 00:00:00 2001 From: ckallidonis Date: Wed, 17 Feb 2021 15:53:43 -0500 Subject: [PATCH 5/6] Removed get-er function for MG environment from Eigsolve class --- include/eigsolve_mugiq.h | 4 ---- lib/loop_mugiq.cpp | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/include/eigsolve_mugiq.h b/include/eigsolve_mugiq.h index 9856b1b..f59bc20 100644 --- a/include/eigsolve_mugiq.h +++ b/include/eigsolve_mugiq.h @@ -169,10 +169,6 @@ class Eigsolve_Mugiq { */ std::vector *getEvals(){ return eVals;} - /** @brief Accessor to get the Multigrid environment structure - */ - MG_Mugiq* getMGEnv(){ return mg_env;} - }; // class Eigsolve_Mugiq diff --git a/lib/loop_mugiq.cpp b/lib/loop_mugiq.cpp index b0acd04..d5a2061 100644 --- a/lib/loop_mugiq.cpp +++ b/lib/loop_mugiq.cpp @@ -281,7 +281,7 @@ void Loop_Mugiq::prolongateEvec(ColorSpinorField *fineEvec, C if(!eigsolve->computeCoarse) errorQuda("%s: Not supposed to be called when computeCoarse is False\n", __func__); if(fieldOrder != QUDA_FLOAT2_FIELD_ORDER) errorQuda("%s: Vector prolongation requires fieldOrder = FLOAT2\n", __func__); - MG_Mugiq &mg_env = *(eigsolve->getMGEnv()); + MG_Mugiq &mg_env = *(eigsolve->mg_env); //- Create one fine and N_coarse temporary coarse fields //- Will be used for prolongating the coarse eigenvectors back to the fine lattice From 80d4f401ac60c0de0d5cdf3854c55c2a34443780 Mon Sep 17 00:00:00 2001 From: ckallidonis Date: Thu, 18 Feb 2021 09:49:13 -0500 Subject: [PATCH 6/6] Fixed issue with fine Dirac operator creation --- include/eigsolve_mugiq.h | 1 + lib/eigsolve_mugiq.cpp | 10 ++++++++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/include/eigsolve_mugiq.h b/include/eigsolve_mugiq.h index f59bc20..9ce6c31 100644 --- a/include/eigsolve_mugiq.h +++ b/include/eigsolve_mugiq.h @@ -74,6 +74,7 @@ class Eigsolve_Mugiq { TimeProfile *eigProfile; // Used for profiling const Dirac *dirac; + const Dirac *diracFine; DiracMatrix *mat; // The Dirac operator whose eigenpairs we are computing DiracMatrix *matFine; // The Dirac operator whose eigenpairs we are computing diff --git a/lib/eigsolve_mugiq.cpp b/lib/eigsolve_mugiq.cpp index 55f3b1c..771b4af 100644 --- a/lib/eigsolve_mugiq.cpp +++ b/lib/eigsolve_mugiq.cpp @@ -14,6 +14,7 @@ Eigsolve_Mugiq::Eigsolve_Mugiq(MugiqEigParam *eigParams_, invParams(eigParams->QudaEigParams->invert_param), eigProfile(eigProfile_), dirac(nullptr), + diracFine(nullptr), mat(nullptr), matFine(nullptr), diracCreated(MUGIQ_BOOL_FALSE), @@ -32,6 +33,9 @@ Eigsolve_Mugiq::Eigsolve_Mugiq(MugiqEigParam *eigParams_, //- This is diracCoarseResidual of the QUDA MG class dirac = mg_env->diracCoarse; if(typeid(*dirac) != typeid(DiracCoarse)) errorQuda("The Coarse Dirac operator must not be preconditioned!\n"); + + //-The Fine Dirac operator + diracFine = mg_env->mg_solver->d; } else{ allocateFineEvecs(); @@ -63,6 +67,7 @@ Eigsolve_Mugiq::Eigsolve_Mugiq(MugiqEigParam *eigParams_, invParams(eigParams->QudaEigParams->invert_param), eigProfile(eigProfile_), dirac(nullptr), + diracFine(nullptr), mat(nullptr), matFine(nullptr), diracCreated(MUGIQ_BOOL_FALSE), @@ -96,7 +101,7 @@ Eigsolve_Mugiq::~Eigsolve_Mugiq(){ delete eVals_sigma; if(mat) delete mat; - if(matFine) delete mat; + if(matFine) delete matFine; mat = nullptr; if(useMGenv){ @@ -208,7 +213,8 @@ void Eigsolve_Mugiq::createNewDiracMatrix(){ else if (eigParams->diracType == MUGIQ_EIG_OPERATOR_MMdag) mat = new DiracMMdag(*dirac); else errorQuda("%s: Unsupported Dirac operator type\n", __func__); - matFine = new DiracM(*dirac); + if(computeCoarse) matFine = new DiracM(*diracFine); + else matFine = new DiracM(*dirac); }