Changed batchedInnerProduct for portability

Grid/algorithms/iterative/RestartedLanczosBidiagonalization.h

@@ -153,7 +153,7 @@ public:
       // Step 1: extend from pStart steps to Nm steps
       // ----------------------------------------------------------------
       extendBasis(pStart, Nm, betaRestart);
-      verify();
+//      verify();
 
       // ----------------------------------------------------------------
       // Step 2: SVD of the Nm x Nm B matrix.
@@ -208,7 +208,7 @@ public:
       // Step 4: implicit restart — compress to Nk steps
       // ----------------------------------------------------------------
       implicitRestart(Nm, Nk, sigma, X, Y, order, betaK, betaRestart);
-      verify();
+//      verify();
 
       // Lucky breakdown: exact invariant subspace found; convergence is exact.
       // B_p^+ = diag(alpha[0..Nk-1]); extract directly from restart basis.
@@ -360,7 +360,7 @@ private:
     if (restart_col >= 0 && restart_col < m && (int)fvec.size() > 0) {
       for (int j = 0; j < restart_col && j < (int)fvec.size(); ++j){
         B(j, restart_col) = fvec[j];
-        std::cout << GridLogMessage << "buildFullB: B  " <<j<<" "<<restart_col<<B(j, restart_col);
+        std::cout << GridLogDebug << "buildFullB: B  " <<j<<" "<<restart_col<<B(j, restart_col)<<std::endl;
       }
     }
     return B;
@@ -464,20 +464,150 @@ private:
     }
   }
 
+public:
   // ------------------------------------------------------------------
-  // Full reorthogonalization of vec against the vectors in basis.
-  // Subtracts projections only — does NOT normalize.
+  // Block reorthogonalization helpers.
+  // Declared public because CUDA extended lambdas cannot live inside
+  // private/protected member functions.
+  //
+  // batchInnerProducts: computes c[j] = <basis[j], vec> for all j
+  //   in a single GPU pass (one accelerator_barrier instead of n).
+  //   Queues n pairs of (per-site kernel, reduceKernel) to computeStream
+  //   without intermediate CPU syncs, then syncs once at the end.
+  //
+  // batchUpdate: computes vec -= sum_j c[j]*basis[j] in one GPU kernel.
+  //
+  // reorthogonalize: two-pass Classical Gram-Schmidt (CGS2) using the
+  //   two helpers above.  Each pass costs 2 GPU syncs (1 IP + 1 update)
+  //   instead of 2n syncs per pass in the old sequential MGS.
+  // ------------------------------------------------------------------
+
+  void batchInnerProducts(const Field &vec,
+                          const std::vector<Field> &basis,
+                          std::vector<ComplexD> &c)
+  {
+    int n = (int)basis.size();
+    c.resize(n);
+    if (n == 0) return;
+
+    typedef typename Field::vector_object         vobj;
+    typedef decltype(innerProduct(vobj(), vobj())) inner_t;
+    typedef decltype(basis[0].View(AcceleratorRead)) View;
+
+    GridBase *grid = vec.Grid();
+    uint64_t oSites = grid->oSites();
+    uint64_t nsimd  = grid->Nsimd();
+
+    // all_ip[j * oSites + ss] = per-site inner product of basis[j] and vec at site ss.
+    // Layout: n contiguous blocks of oSites each.
+    deviceVector<inner_t> all_ip((uint64_t)n * oSites);
+    inner_t *all_ip_p = &all_ip[0];
+
+    hostVector<View>   h_basis_v(n);
+    deviceVector<View> d_basis_v(n);
+    for (int j = 0; j < n; ++j) {
+      h_basis_v[j] = basis[j].View(AcceleratorRead);
+      acceleratorPut(d_basis_v[j], h_basis_v[j]);
+    }
+    View *basis_vp = &d_basis_v[0];
+
+    // Queue n per-site kernels to the accelerator stream — no intermediate barriers.
+    autoView(vec_v, vec, AcceleratorRead);
+    for (int j = 0; j < n; ++j) {
+      int      jj      = j;
+      uint64_t oSites_ = oSites;
+      accelerator_for(ss, oSites, nsimd, {
+        auto x = coalescedRead(basis_vp[jj][ss]);
+        auto y = coalescedRead(vec_v[ss]);
+        coalescedWrite(all_ip_p[jj * oSites_ + ss], innerProduct(x, y));
+      });
+    }
+
+    // ONE sync after all n kernels
+    accelerator_barrier();
+
+    // Copy all per-site results to host
+    hostVector<inner_t> all_ip_h((uint64_t)n * oSites);
+    acceleratorCopyFromDevice(all_ip_p, &all_ip_h[0], (uint64_t)n * oSites * sizeof(inner_t));
+
+    // Reduce on host: sum over oSites, then collapse SIMD lanes via Reduce(TensorRemove(...))
+    // TensorRemove strips the iSinglet tensor wrapper to expose the SIMD scalar type.
+    // Reduce sums all nsimd lanes and returns a plain scalar (RealD or ComplexD).
+    std::vector<ComplexD> raw(n);
+    for (int j = 0; j < n; ++j) {
+      inner_t sum = Zero();
+      for (uint64_t ss = 0; ss < oSites; ++ss)
+        sum += all_ip_h[(uint64_t)j * oSites + ss];
+      raw[j] = ComplexD(Reduce(TensorRemove(sum)));
+    }
+    grid->GlobalSumVector(&raw[0], n);
+    for (int j = 0; j < n; ++j) c[j] = raw[j];
+
+    for (int j = 0; j < n; ++j) h_basis_v[j].ViewClose();
+  }
+
+  void batchUpdate(Field &vec,
+                   const std::vector<Field> &basis,
+                   const std::vector<ComplexD> &c)
+  {
+    int n = (int)basis.size();
+    if (n == 0) return;
+
+    typedef typename Field::vector_object vobj;
+    typedef decltype(basis[0].View(AcceleratorRead)) View;
+
+    GridBase *grid = vec.Grid();
+    uint64_t oSites = grid->oSites();
+    uint64_t nsimd  = grid->Nsimd();
+
+    // Split complex coefficients into real/imag double arrays on device.
+    // Using doubles avoids potential ComplexD-device-code compatibility issues.
+    hostVector<double>   h_re(n), h_im(n);
+    deviceVector<double> d_re(n), d_im(n);
+    for (int k = 0; k < n; ++k) {
+      h_re[k] = c[k].real();
+      h_im[k] = c[k].imag();
+    }
+    acceleratorCopyToDevice(&h_re[0], &d_re[0], n * sizeof(double));
+    acceleratorCopyToDevice(&h_im[0], &d_im[0], n * sizeof(double));
+    double *re_p = &d_re[0];
+    double *im_p = &d_im[0];
+
+    // Basis views
+    hostVector<View>   h_basis_v(n);
+    deviceVector<View> d_basis_v(n);
+    for (int k = 0; k < n; ++k) {
+      h_basis_v[k] = basis[k].View(AcceleratorRead);
+      acceleratorPut(d_basis_v[k], h_basis_v[k]);
+    }
+    View *basis_vp = &d_basis_v[0];
+
+    // Single kernel: vec[ss] -= sum_k (re[k] + i*im[k]) * basis[k][ss]
+    autoView(vec_v, vec, AcceleratorWrite);
+    accelerator_for(ss, oSites, nsimd, {
+      auto v = coalescedRead(vec_v[ss]);
+      for (int k = 0; k < n; ++k) {
+        auto b = coalescedRead(basis_vp[k][ss]);
+        v = v - re_p[k] * b - timesI(im_p[k] * b);
+      }
+      coalescedWrite(vec_v[ss], v);
+    });
+
+    for (int k = 0; k < n; ++k) h_basis_v[k].ViewClose();
+  }
+
+  // ------------------------------------------------------------------
+  // Full reorthogonalization using two-pass Classical Gram-Schmidt (CGS2).
+  // Each pass calls batchInnerProducts (1 GPU sync) + batchUpdate (1 sync),
+  // replacing the old 2n GPU syncs per pass from sequential MGS.
   // ------------------------------------------------------------------
   void reorthogonalize(Field &vec, const std::vector<Field> &basis)
   {
-    for (int j = 0; j < (int)basis.size(); ++j) {
-      ComplexD ip = innerProduct(basis[j], vec);
-      vec = vec - ip * basis[j];
-    }
-    // Second pass for numerical stability
-    for (int j = 0; j < (int)basis.size(); ++j) {
-      ComplexD ip = innerProduct(basis[j], vec);
-      vec = vec - ip * basis[j];
+    if (basis.empty()) return;
+    std::vector<ComplexD> c;
+    for (int pass = 0; pass < 2; ++pass) {
+      batchInnerProducts(vec, basis, c);
+      batchUpdate(vec, basis, c);
     }
   }
 

tests/lanczos/Test_wilson_specflow.cc

@@ -200,8 +200,8 @@ int main(int argc, char** argv) {
     RD.close();
   }
 
-//  std::vector<Complex> boundary = {1,1,1,-1};
-  std::vector<Complex> boundary = {1,1,1,1};
+  std::vector<Complex> boundary = {1,1,1,-1};
+//  std::vector<Complex> boundary = {1,1,1,1};
   FermionOp::ImplParams Params(boundary);