merge hippo and aplus cuda pcg kernels

Author: zhi-wang

include/ff/aplus/induce.h

@@ -6,8 +6,5 @@ void dfieldAplus(real (*field)[3]);
 
 void ufieldAplus(const real (*uind)[3], real (*field)[3]);
 
-void diagPrecond3(const real (*rsd)[3], real (*zrsd)[3]);
-
-void sparsePrecondBuild3();
 void sparsePrecondApply3(const real (*rsd)[3], real (*zrsd)[3]);
 }

src/acc/aplus/induce.cpp

@@ -1,239 +1,10 @@
-#include "ff/amoeba/induce.h"
 #include "ff/amoebamod.h"
-#include "ff/aplus/induce.h"
 #include "ff/atom.h"
 #include "ff/image.h"
 #include "ff/nblist.h"
-#include "ff/switch.h"
 #include "seq/add.h"
 #include "seq/damp.h"
-#include "tool/error.h"
 #include "tool/gpucard.h"
-#include "tool/ioprint.h"
-#include <tinker/detail/inform.hh>
-#include <tinker/detail/polpcg.hh>
-#include <tinker/detail/polpot.hh>
-#include <tinker/detail/units.hh>
-
-namespace tinker {
-void induceMutualPcg3_acc(real (*uind)[3])
-{
-   auto* field = work01_;
-   auto* rsd = work02_;
-   auto* zrsd = work03_;
-   auto* conj = work04_;
-   auto* vec = work05_;
-
-   bool dirguess = polpcg::pcgguess;
-   // use sparse matrix preconditioner
-   // or just use diagonal matrix preconditioner
-   const bool sparse_prec = polpcg::pcgprec and (switchOff(Switch::USOLVE) > 0);
-
-   // zero out the induced dipoles at each site
-
-   bool predict = polpred != UPred::NONE;
-   if (predict and nualt < maxualt) {
-      predict = false;
-      dirguess = true;
-   }
-   // get the electrostatic field due to permanent multipoles
-   dfieldAplus(field);
-   // direct induced dipoles
-
-   #pragma acc parallel loop independent async\
-               deviceptr(polarity,udir,field)
-   for (int i = 0; i < n; ++i) {
-      real poli = polarity[i];
-      #pragma acc loop seq
-      for (int j = 0; j < 3; ++j)
-         udir[i][j] = poli * field[i][j];
-   }
-
-   // initial induced dipole
-   if (predict) {
-      ulspredSum(uind, nullptr);
-   } else if (dirguess) {
-      darray::copy(g::q0, n, uind, udir);
-   } else {
-      darray::zero(g::q0, n, uind);
-   }
-   // initial residual r(0)
-
-   // if do not use pcgguess, r(0) = E - T Zero = E
-   // if use pcgguess, r(0) = E - (inv_alpha + Tu) alpha E
-   //                       = E - E -Tu udir
-   //                       = -Tu udir
-
-   if (predict) {
-      ufieldAplus(uind, field);
-      #pragma acc parallel loop independent async\
-              deviceptr(polarity_inv,rsd,udir,uind,field)
-      for (int i = 0; i < n; ++i) {
-         real pol = polarity_inv[i];
-         #pragma acc loop seq
-         for (int j = 0; j < 3; ++j)
-            rsd[i][j] = (udir[i][j] - uind[i][j]) * pol + field[i][j];
-      }
-   } else if (dirguess) {
-      ufieldAplus(uind, rsd);
-   } else {
-      darray::copy(g::q0, n, rsd, field);
-   }
-   #pragma acc parallel loop independent async deviceptr(polarity,rsd)
-   for (int i = 0; i < n; ++i) {
-      if (polarity[i] == 0) {
-         rsd[i][0] = 0;
-         rsd[i][1] = 0;
-         rsd[i][2] = 0;
-      }
-   }
-
-   // initial M r(0) and p(0)
-
-   if (sparse_prec) {
-      sparsePrecondBuild3();
-      sparsePrecondApply3(rsd, zrsd);
-   } else {
-      diagPrecond3(rsd, zrsd);
-   }
-   darray::copy(g::q0, n, conj, zrsd);
-
-   // initial r(0) M r(0)
-
-   real sum;
-   sum = darray::dotThenReturn(g::q0, n, rsd, zrsd);
-
-   // conjugate gradient iteration of the mutual induced dipoles
-
-   const bool debug = inform::debug;
-   const int politer = polpot::politer;
-   const real poleps = polpot::poleps;
-   const real debye = units::debye;
-   const real pcgpeek = polpcg::pcgpeek;
-   const int maxiter = 100; // see also subroutine induce0a in induce.f
-   const int miniter = std::min(3, n);
-
-   bool done = false;
-   int iter = 0;
-   real eps = 100;
-   // real epsold;
-
-   while (not done) {
-      ++iter;
-
-      // vec = (inv_alpha + Tu) conj, field = -Tu conj
-      // vec = inv_alpha * conj - field
-      ufieldAplus(conj, field);
-      #pragma acc parallel loop independent async\
-                  deviceptr(polarity_inv,vec,conj,field)
-      for (int i = 0; i < n; ++i) {
-         real poli_inv = polarity_inv[i];
-         #pragma acc loop seq
-         for (int j = 0; j < 3; ++j)
-            vec[i][j] = poli_inv * conj[i][j] - field[i][j];
-      }
-
-      // a <- p T p
-      real a;
-      a = darray::dotThenReturn(g::q0, n, conj, vec);
-      // a <- r M r / p T p
-      if (a != 0)
-         a = sum / a;
-
-      // u <- u + a p
-      // r <- r - a T p
-      #pragma acc parallel loop independent async\
-                  deviceptr(polarity,uind,conj,rsd,vec)
-      for (int i = 0; i < n; ++i) {
-         #pragma acc loop seq
-         for (int j = 0; j < 3; ++j) {
-            uind[i][j] += a * conj[i][j];
-            rsd[i][j] -= a * vec[i][j];
-         }
-         if (polarity[i] == 0) {
-            rsd[i][0] = 0;
-            rsd[i][1] = 0;
-            rsd[i][2] = 0;
-         }
-      }
-
-      // calculate/update M r
-      if (sparse_prec)
-         sparsePrecondApply3(rsd, zrsd);
-      else
-         diagPrecond3(rsd, zrsd);
-
-      real b;
-      real sum1;
-      sum1 = darray::dotThenReturn(g::q0, n, rsd, zrsd);
-      b = sum1 / sum;
-      if (sum == 0)
-         b = 0;
-
-      // calculate/update p
-      #pragma acc parallel loop independent async\
-                  deviceptr(conj,zrsd)
-      for (int i = 0; i < n; ++i) {
-         #pragma acc loop seq
-         for (int j = 0; j < 3; ++j)
-            conj[i][j] = zrsd[i][j] + b * conj[i][j];
-      }
-
-      sum = sum1;
-
-      real epsd;
-      epsd = darray::dotThenReturn(g::q0, n, rsd, rsd);
-      // epsold = eps;
-      eps = epsd;
-      eps = debye * REAL_SQRT(eps / n);
-
-      if (debug) {
-         if (iter == 1) {
-            print(stdout,
-               "\n Determination of SCF Induced Dipole Moments\n\n"
-               "    Iter    RMS Residual (Debye)\n\n");
-         }
-         print(stdout, " %8d       %-16.10f\n", iter, eps);
-      }
-
-      if (eps < poleps)
-         done = true;
-      // if (eps > epsold) done = true;
-      if (iter < miniter) done = false;
-      if (iter >= politer)
-         done = true;
-
-      // apply a "peek" iteration to the mutual induced dipoles
-
-      if (done) {
-         #pragma acc parallel loop independent async\
-                     deviceptr(polarity,uind,rsd)
-         for (int i = 0; i < n; ++i) {
-            real term = pcgpeek * polarity[i];
-            #pragma acc loop seq
-            for (int j = 0; j < 3; ++j)
-               uind[i][j] += term * rsd[i][j];
-         }
-      }
-   }
-
-   // print the results from the conjugate gradient iteration
-
-   if (debug) {
-      print(stdout,
-         " Induced Dipoles :    Iterations %4d      RMS "
-         "Residual %14.10f\n",
-         iter, eps);
-   }
-
-   // terminate the calculation if dipoles failed to converge
-
-   if (iter >= maxiter) {
-      printError();
-      TINKER_THROW("INDUCE  --  Warning, Induced Dipoles are not Converged");
-   }
-}
-}
 
 namespace tinker {
 #define APPLY_DPTRS rsd, zrsd, x, y, z, polarity, pdamp, thole

src/acc/hippo/induce.cpp

@@ -1,5 +1,6 @@
 #include "ff/amoeba/induce.h"
 #include "ff/amoebamod.h"
+#include "ff/aplus/induce.h"
 #include "ff/atom.h"
 #include "ff/hippo/induce.h"
 #include "ff/hippomod.h"
@@ -17,26 +18,7 @@
 #include <tinker/detail/polpot.hh>
 #include <tinker/detail/units.hh>
 
-namespace tinker { // PCG
-//
-// M = preconditioner
-// T = inv_alpha + Tu, -Tu = ufield_chgpen
-//
-// r(0) = E - T u(0)
-// p(0) (or conj(0)) = M r(0)
-// ------------------------------
-// gamma (or a) = r(i) M r(i) / p(i) T p(i)
-// u(i+1) = u(i) + a p(i)
-// r(i+1) = r(i) - a T p(i)
-// beta (or b(i+1)) = r(i+1) M r(i+1) / r(i) M r(i)
-// p(i+1) = M r(r+1) + b(i+1) p(i)
-// ------------------------------
-//
-// subroutine induce0a in induce.f
-// rsd = r
-// zrsd = M r
-// conj = p
-// vec = T P
+namespace tinker {
 void induceMutualPcg2_acc(real (*uind)[3])
 {
    auto* field = work01_;
@@ -45,6 +27,10 @@ void induceMutualPcg2_acc(real (*uind)[3])
    auto* conj = work04_;
    auto* vec = work05_;
 
+   auto pcg_dfield = polpot::use_tholed ? dfieldAplus : dfieldChgpen;
+   auto pcg_ufield = polpot::use_tholed ? ufieldAplus : ufieldChgpen;
+   auto pcg_sparse_apply = polpot::use_tholed ? sparsePrecondApply3 : sparsePrecondApply2;
+
    bool dirguess = polpcg::pcgguess;
    // use sparse matrix preconditioner
    // or just use diagonal matrix preconditioner
@@ -58,7 +44,7 @@ void induceMutualPcg2_acc(real (*uind)[3])
       dirguess = true;
    }
    // get the electrostatic field due to permanent multipoles
-   dfieldChgpen(field);
+   pcg_dfield(field);
    // direct induced dipoles
 
    #pragma acc parallel loop independent async\
@@ -86,7 +72,7 @@ void induceMutualPcg2_acc(real (*uind)[3])
    //                       = -Tu udir
 
    if (predict) {
-      ufieldChgpen(uind, field);
+      pcg_ufield(uind, field);
       #pragma acc parallel loop independent async\
               deviceptr(polarity_inv,rsd,udir,uind,field)
       for (int i = 0; i < n; ++i) {
@@ -96,7 +82,7 @@ void induceMutualPcg2_acc(real (*uind)[3])
             rsd[i][j] = (udir[i][j] - uind[i][j]) * pol + field[i][j];
       }
    } else if (dirguess) {
-      ufieldChgpen(uind, rsd);
+      pcg_ufield(uind, rsd);
    } else {
       darray::copy(g::q0, n, rsd, field);
    }
@@ -113,7 +99,7 @@ void induceMutualPcg2_acc(real (*uind)[3])
 
    if (sparse_prec) {
       sparsePrecondBuild2();
-      sparsePrecondApply2(rsd, zrsd);
+      pcg_sparse_apply(rsd, zrsd);
    } else {
       diagPrecond2(rsd, zrsd);
    }
@@ -144,7 +130,7 @@ void induceMutualPcg2_acc(real (*uind)[3])
 
       // vec = (inv_alpha + Tu) conj, field = -Tu conj
       // vec = inv_alpha * conj - field
-      ufieldChgpen(conj, field);
+      pcg_ufield(conj, field);
       #pragma acc parallel loop independent async\
                   deviceptr(polarity_inv,vec,conj,field)
       for (int i = 0; i < n; ++i) {
@@ -180,7 +166,7 @@ void induceMutualPcg2_acc(real (*uind)[3])
 
       // calculate/update M r
       if (sparse_prec)
-         sparsePrecondApply2(rsd, zrsd);
+         pcg_sparse_apply(rsd, zrsd);
       else
          diagPrecond2(rsd, zrsd);
 
@@ -220,7 +206,8 @@ void induceMutualPcg2_acc(real (*uind)[3])
       if (eps < poleps)
          done = true;
       // if (eps > epsold) done = true;
-      if (iter < miniter) done = false;
+      if (iter < miniter)
+         done = false;
       if (iter >= politer)
          done = true;
 

src/aplus/induce.cpp

@@ -4,13 +4,6 @@
 #include "tool/externfunc.h"
 
 namespace tinker {
-void diagPrecond3(const real (*rsd)[3], real (*zrsd)[3])
-{
-   diagPrecond2(rsd, zrsd);
-}
-
-void sparsePrecondBuild3() {}
-
 TINKER_FVOID2(acc1, cu1, sparsePrecondApply3, const real (*)[3], real (*)[3]);
 void sparsePrecondApply3(const real (*rsd)[3], real (*zrsd)[3])
 {