From 8fd16686dcae1db1f314a7b50ddb2fed577e8cc2 Mon Sep 17 00:00:00 2001
From: Chulwoo Jung <chulwoo@bnl.gov>
Date: Sun, 5 Sep 2021 23:04:41 -0400
Subject: [PATCH] Checking block lanczos deleted some diag outputs

---
 .../ImplicitlyRestartedBlockLanczos.h         | 42 ++++++++---
 Grid/allocator/MemoryManager.cc               | 10 +--
 tests/lanczos/Test_dwf_block_lanczos.README   | 73 +++++++++++++++++++
 3 files changed, 108 insertions(+), 17 deletions(-)
 create mode 100644 tests/lanczos/Test_dwf_block_lanczos.README

diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
index 809d64d9..1b60f89f 100644
--- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
+++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
@@ -721,7 +721,7 @@ cudaStat = cudaMallocManaged((void **)&evec_acc, Nevec_acc*sites*12*sizeof(CUDA_
     Eigen::MatrixXcd    Q = Eigen::MatrixXcd::Zero(Nm,Nm);
 
     std::vector<int>   Iconv(Nm);
-    int Ntest=Nu;
+//    int Ntest=Nu;
 //    std::vector<Field>  B(Nm,grid); // waste of space replicating
     std::vector<Field>  B(1,grid); // waste of space replicating
     
@@ -731,7 +731,7 @@ cudaStat = cudaMallocManaged((void **)&evec_acc, Nevec_acc*sites*12*sizeof(CUDA_
     
     Nconv = 0;
     
-    RealD beta_k;
+//    RealD beta_k;
   
     // set initial vector
     for (int i=0; i<Nu; ++i) {
@@ -838,22 +838,40 @@ cudaStat = cudaMallocManaged((void **)&evec_acc, Nevec_acc*sites*12*sizeof(CUDA_
     } else {
       Glog << fname + " CONVERGED ; Summary :\n";
       // Sort convered eigenpairs.
-      std::vector<Field>  Btmp(Nconv,grid); // waste of space replicating
+      std::vector<Field>  Btmp(Nstop,grid); // waste of space replicating
+#if 0
       for(int i=0; i<Nconv; ++i) Btmp[i]=0;
       for(int i=0; i<Nconv; ++i)
       for(int k = 0; k<Nr; ++k){
          Btmp[i].Checkerboard() = evec[k].Checkerboard();
-         Btmp[i] += evec[k]*Qt(k,Iconv[i]);
+         Btmp[i] += evec[k]*Qt(k,i);
       }
+#endif
 
-      eval.resize(Nconv);
-      evec.resize(Nconv,grid);
-
-      for(int i=0; i<Nconv; ++i){
-        eval[i] = eval2[Iconv[i]];
-        evec[i] = Btmp[i];
+      for(int i=0; i<Nstop; ++i){
+	  Btmp[i]=0.;
+          for(int k = 0; k<Nr; ++k){
+             Btmp[i].Checkerboard() = evec[k].Checkerboard();
+             Btmp[i] += evec[k]*Qt(k,i);
+          }
+          _Linop.HermOp(Btmp[i],v);
+          RealD vnum = real(innerProduct(Btmp[i],v)); // HermOp.
+          RealD vden = norm2(Btmp[i]);
+//          eval2[j] = vnum/vden;
+          v -= vnum/vden*Btmp[i];
+          RealD vv = norm2(v);
+//          resid[j] = vv;
+          
+          std::cout.precision(13);
+          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+          std::cout << "eval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< vnum/vden<<" "<<eval2[i];;
+          std::cout << "   resid^2 = "<< std::setw(20)<< std::setiosflags(std::ios_base::right)<< vv<< std::endl;
+        eval[i] = vnum/vden;
       }
-      _sort.push(eval,evec,Nconv);
+      for(int i=0; i<Nstop; ++i) evec[i] = Btmp[i];
+      eval.resize(Nstop);
+      evec.resize(Nstop,grid);
+      _sort.push(eval,evec,Nstop);
     }
     Glog << std::string(74,'*') << std::endl;
     Glog << " -- Iterations    = "<< iter   << "\n";
@@ -972,7 +990,7 @@ if(split_test){
     }
 
     // re-orthogonalization for numerical stability
-#if 0
+#if 1
     Glog << "Gram Schmidt"<< std::endl;
     orthogonalize(w,Nu,evec,R);
 #else
diff --git a/Grid/allocator/MemoryManager.cc b/Grid/allocator/MemoryManager.cc
index bfa06938..32eecbec 100644
--- a/Grid/allocator/MemoryManager.cc
+++ b/Grid/allocator/MemoryManager.cc
@@ -44,7 +44,7 @@ void  MemoryManager::AcceleratorFree    (void *ptr,size_t bytes)
   if ( __freeme ) {
     acceleratorFreeDevice(__freeme);
     total_device-=bytes;
-       PrintBytes();
+//       PrintBytes();
   }
 }
 void *MemoryManager::SharedAllocate(size_t bytes)
@@ -53,8 +53,8 @@ void *MemoryManager::SharedAllocate(size_t bytes)
   if ( ptr == (void *) NULL ) {
     ptr = (void *) acceleratorAllocShared(bytes);
     total_shared+=bytes;
-        std::cout <<"AcceleratorAllocate: allocated Shared pointer "<<std::hex<<ptr<<std::dec<<std::endl;
-        PrintBytes();
+        std::cout <<GridLogMessage<<"AcceleratorAllocate: allocated Shared pointer "<<std::hex<<ptr<<std::dec<<std::endl;
+//        PrintBytes();
   }
   return ptr;
 }
@@ -74,7 +74,7 @@ void *MemoryManager::CpuAllocate(size_t bytes)
   if ( ptr == (void *) NULL ) {
     ptr = (void *) acceleratorAllocShared(bytes);
     total_host+=bytes;
-    std::cout << GridLogMessage<< "MemoryManager:: CpuAllocate  total_host= "<<total_host<<" "<< ptr << std::endl;
+//    std::cout << GridLogMessage<< "MemoryManager:: CpuAllocate  total_host= "<<total_host<<" "<< ptr << std::endl;
   }
   return ptr;
 }
@@ -84,7 +84,7 @@ void  MemoryManager::CpuFree    (void *_ptr,size_t bytes)
   void *__freeme = Insert(_ptr,bytes,Cpu);
   if ( __freeme ) { 
     acceleratorFreeShared(__freeme);
-    std::cout << GridLogMessage<< "MemoryManager:: CpuFree  total_host= "<<total_host<<" "<< __freeme << std::endl;
+//    std::cout << GridLogMessage<< "MemoryManager:: CpuFree  total_host= "<<total_host<<" "<< __freeme << std::endl;
     total_host-=bytes;
   }
 }
diff --git a/tests/lanczos/Test_dwf_block_lanczos.README b/tests/lanczos/Test_dwf_block_lanczos.README
new file mode 100644
index 00000000..179f9037
--- /dev/null
+++ b/tests/lanczos/Test_dwf_block_lanczos.README
@@ -0,0 +1,73 @@
+#Example script 
+DIR=/gpfs/alpine/phy157/proj-shared/phy157dwf/chulwoo/Grid/BL/build/tests/lanczos
+BIN=${DIR}/Test_dwf_block_lanczos
+
+VOL='--grid 16.16.16.32 '
+GRID='--mpi 1.1.1.4 ' 
+CONF='--gconf ckpoint_lat.IEEE64BIG.2000 '
+OPT='--mass 0.01 --M5 1.8 --phase in.params --omega in.params --shm 4096' 
+#BL='--rbl 16.1024.128.1000.10 --split 1.1.4.4 --check_int 100 --resid 1.0e-5 --cheby_l 0.007 --cheby_u 7 --cheby_n 51'
+BL='--rbl 4.128.16.100.10 --split 1.1.1.4 --check_int 25 --resid 1.0e-5 --cheby_l 0.007 --cheby_u 7 --cheby_n 51'
+
+ARGS=${CONF}" "${OPT}" "${BL}" "${VOL}" "${GRID}
+export APP="${BIN}  ${ARGS}"
+echo APP=${APP}
+#export JS="jsrun --nrs 32 -a4 -g4 -c42 -dpacked -b packed:7 --smpiargs="-gpu" "
+export JS="jsrun --nrs 1 -a4 -g4 -c42 -dpacked -b  packed:10  --smpiargs="-gpu" "
+$JS  $APP
+
+#sample in.param 
+
+boundary_phase 0 1 0
+boundary_phase 1 1 0
+boundary_phase 2 1 0
+boundary_phase 3 -1 0
+
+omega 0 0.5 0
+omega 1 0.5 0
+omega 2 0.5 0
+omega 3 0.5 0
+omega 4 0.5 0
+omega 5 0.5 0
+omega 6 0.5 0
+omega 7 0.5 0
+omega 8 0.5 0
+omega 9 0.5 0
+omega 10 0.5 0
+omega 11 0.5 0
+
+
+#output 
+
+Grid : Message : 1.717474 s :  Gauge Configuration ckpoint_lat.IEEE64BIG.2000
+Grid : Message : 1.717478 s :  boundary_phase[0] = (1,0)
+Grid : Message : 1.717497 s :  boundary_phase[1] = (1,0)
+Grid : Message : 1.717500 s :  boundary_phase[2] = (1,0)
+Grid : Message : 1.717503 s :  boundary_phase[3] = (-1,0)
+Grid : Message : 1.717506 s :  Ls 12
+Grid : Message : 1.717507 s :  mass 0.01
+Grid : Message : 1.717510 s :  M5 1.8
+Grid : Message : 1.717512 s :  mob_b 1.5
+Grid : Message : 1.717514 s :  omega[0] = (0.5,0)
+Grid : Message : 1.717517 s :  omega[1] = (0.5,0)
+Grid : Message : 1.717520 s :  omega[2] = (0.5,0)
+Grid : Message : 1.717523 s :  omega[3] = (0.5,0)
+Grid : Message : 1.717526 s :  omega[4] = (0.5,0)
+Grid : Message : 1.717529 s :  omega[5] = (0.5,0)
+Grid : Message : 1.717532 s :  omega[6] = (0.5,0)
+Grid : Message : 1.717535 s :  omega[7] = (0.5,0)
+Grid : Message : 1.717538 s :  omega[8] = (0.5,0)
+Grid : Message : 1.717541 s :  omega[9] = (0.5,0)
+Grid : Message : 1.717544 s :  omega[10] = (0.5,0)
+Grid : Message : 1.717547 s :  omega[11] = (0.5,0)
+Grid : Message : 1.717550 s :  Nu 4
+Grid : Message : 1.717551 s :  Nk 128
+Grid : Message : 1.717552 s :  Np 16
+Grid : Message : 1.717553 s :  Nm 288
+Grid : Message : 1.717554 s :  Nstop 100
+Grid : Message : 1.717555 s :  Ntest 25
+Grid : Message : 1.717557 s :  MaxIter 10
+Grid : Message : 1.717558 s :  resid 1e-05
+Grid : Message : 1.717560 s :  Cheby Poly 0.007,7,51
+
+