Updates to Pull Request #9 - Oblique projection

include/eigsolve_mugiq.h

@@ -74,7 +74,9 @@ 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
 
   //- This switch is required so that the dirac object is NOT
   //- deleted when NOT created within Eigsolve.
@@ -151,50 +153,23 @@ 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);
 
   /** @brief Compute eigenvalues
    */
   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<ColorSpinorField *> &getEvecs(){ return eVecs;}
 
-  /** @brief Accessor to get the Quda eigenvalues outside of the class
-   */
-  std::vector<Complex>* 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<Complex> *getEvals(){ return eVals;}
 
-  /** @brief Accessor to get the Eigsolve_Mugiq singular values outside of the class
-   */
-  std::vector<double> *getEvalsSigma(){ return eVals_sigma;}
-
-  /** @brief Accessor to get the residual of the computed eigenvalues
-   */
-  std::vector<double>* 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 
 
 

lib/eigsolve_mugiq.cpp

@@ -14,7 +14,9 @@ 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),
   eVals_quda(nullptr),
   eVals(nullptr),
@@ -31,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();
@@ -62,7 +67,9 @@ 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),
   eVals_quda(nullptr),
   eVals(nullptr),
@@ -94,6 +101,7 @@ Eigsolve_Mugiq::~Eigsolve_Mugiq(){
     delete eVals_sigma;
 
   if(mat) delete mat;
+  if(matFine) delete matFine;
   mat = nullptr;
 
   if(useMGenv){
@@ -204,6 +212,9 @@ 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__);
+
+  if(computeCoarse) matFine = new DiracM(*diracFine);
+  else matFine = new DiracM(*dirac);
 }
 
 
@@ -283,6 +294,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; i<eigParams->nEv; i++){
+      gamma5(*eVecs[i], *eVecs[i]);
+    }
+  }
+
   delete eigSolve;
 }
 
@@ -296,21 +314,22 @@ void Eigsolve_Mugiq::computeEvals(){
   std::vector<double> &r = *evals_res;
 
   double kappa = invParams->kappa;
+  DiracMatrix *mat_direct = new DiracM(*dirac);	  
 
   for(int i=0; i<eigParams->nEv; 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||
   }
 
-  if(eigParams->diracType == MUGIQ_EIG_OPERATOR_MdagM || eigParams->diracType == MUGIQ_EIG_OPERATOR_MMdag){
-    std::vector<double> &sigma = *eVals_sigma;    
-    for(int i=0; i<eigParams->nEv; i++) sigma[i] = sqrt(lambda[i].real());
-  }
+  std::vector<double> &sigma = *eVals_sigma;    
+  for(int i=0; i<eigParams->nEv; i++) sigma[i] = lambda[i].real();
 
+  delete mat_direct;
   delete w;
 }
 
@@ -335,14 +354,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; i<eigParams->nEv; 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; lev<mg_env->nCoarseLevels; 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<Complex> s(eigParams->nEv);
+  std::vector<ColorSpinorField *> 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<ColorSpinorField *> 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;
 }

lib/loop_mugiq.cpp

@@ -281,7 +281,7 @@ void Loop_Mugiq<Float, fieldOrder>::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