From b89b1280d5d384ce10ea8fc77eae1c804d2e9569 Mon Sep 17 00:00:00 2001
From: Yong-Chull Jang <integration.field@gmail.com>
Date: Tue, 31 Mar 2020 05:39:31 -0400
Subject: [PATCH] use gemm twice to complete the Gram Schmidt

---
 .../ImplicitlyRestartedBlockLanczos.h         | 23 ++++++++++++++++---
 1 file changed, 20 insertions(+), 3 deletions(-)

diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
index cd6bd682..95fedbca 100644
--- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
@@ -303,7 +303,7 @@ public:
       //Glog << stat << std::endl;
       
       grid->GlobalSumVector((double*)c_acc,2*Nu*Nevec_acc);
-      
+#if 0      
       for (int i=0; i<Nu; ++i) {
         for (size_t j=0; j<Nevec_acc; ++j) {
           cuDoubleComplex z = c_acc[i*Nevec_acc+j];
@@ -314,10 +314,27 @@ public:
           }
           w[i] = w[i] - ip * evec[b*Nevec_acc+j];
         }
-        //assert(normalize(w[i],do_print)!=0);
+      }
+#else
+      alpha = make_cuDoubleComplex(-1.0,0.0);
+      beta = make_cuDoubleComplex(1.0,0.0);
+      stat = cublasZgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, 12*sites, Nu, Nevec_acc,
+                         &alpha, 
+                         evec_acc, 12*sites, c_acc, Nevec_acc,  
+                         &beta, 
+                         w_acc, 12*sites);
+      //Glog << stat << std::endl;
+#endif
+    }
+#if 1    
+    for (int col=0; col<Nu; ++col) {
+      auto w_v = w[col].View();
+      cuDoubleComplex *z = reinterpret_cast<cuDoubleComplex*>(&w_v[0]);
+      for (size_t row=0; row<sites*12; ++row) {
+        z[row] = w_acc[col*sites*12+row];
       }
     }
-
+#endif
     for (int i=0; i<Nu; ++i) {
       assert(normalize(w[i],do_print)!=0);
     }