From ca4d1f01786a312da106a9844a948764f26fa833 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E2=80=9Chborchert=E2=80=9D?=
 <“hannes.borchert@stonybrook.edu”>
Date: Wed, 14 Aug 2024 17:21:15 -0400
Subject: [PATCH] added tiling of apply BSH in SCF.cc

---
 src/madness/chem/SCF.cc             | 47 ++++++++++++++++++++++++-----
 src/madness/chem/exchangeoperator.h |  4 +--
 2 files changed, 41 insertions(+), 10 deletions(-)

diff --git a/src/madness/chem/SCF.cc b/src/madness/chem/SCF.cc
index 653118cc8e2..805b312fc9b 100644
--- a/src/madness/chem/SCF.cc
+++ b/src/madness/chem/SCF.cc
@@ -1549,26 +1549,57 @@ vecfuncT SCF::compute_residual(World& world, tensorT& occ, tensorT& fock,
     fpsi.clear();
     std::vector<double> fac(nmo, -2.0);
     scale(world, Vpsi, fac);
-    std::vector<poperatorT> ops = make_bsh_operators(world, eps);
-    set_thresh(world, Vpsi, FunctionDefaults<3>::get_thresh());
+    //std::vector<poperatorT> ops = make_bsh_operators(world, eps);
+    //set_thresh(world, Vpsi, FunctionDefaults<3>::get_thresh());
     END_TIMER(world, "Compute residual stuff");
 
     START_TIMER(world);
 
-    //TODO: tile apply
-    vecfuncT new_psi = apply(world, ops, Vpsi);
+    //vecfuncT new_psi = apply(world, ops, Vpsi);
+
+    // TODO: tile apply
+    size_t min_tile = 10;
+    size_t ntile = std::min(amo.size(), min_tile);
+    vecfuncT new_psi = zero_functions<double,3>(world, Vpsi.size());
+
+    for (size_t ilo=0; ilo<Vpsi.size(); ilo+=ntile) {
+        size_t iend = std::min(ilo+ntile,Vpsi.size());
+        vecfuncT tmp_Vpsi(Vpsi.begin()+ilo,Vpsi.begin()+iend);
+
+        int tmp_nmo = tmp_Vpsi.size();
+        tensorT tmp_eps(tmp_nmo);
+        for (int i = 0; i < tmp_nmo; ++i) {
+            tmp_eps(i) = std::min(-0.05, fock(i+ilo, i+ilo));
+        }
+
+        std::vector<poperatorT> ops = make_bsh_operators(world, tmp_eps);
+        set_thresh(world, tmp_Vpsi, FunctionDefaults<3>::get_thresh());
+
+        vecfuncT tmp_new_psi = apply(world, ops, tmp_Vpsi);
+        print_size(world, tmp_new_psi, "tmp_new_psi before truncation");
+
+        //truncate tmp_new_psi
+        truncate(world, tmp_new_psi);
+        print_size(world, tmp_new_psi, "tmp_new_psi after truncation");
+
+        //put the results into their final home
+        for (size_t i = ilo; i<iend; ++i){
+            new_psi[i] += tmp_new_psi[i-ilo];
+        }
+        ops.clear();
+    }
 
     END_TIMER(world, "Apply BSH");
-    ops.clear();
+    //ops.clear();
     Vpsi.clear();
     world.gop.fence();
 
     // Thought it was a bad idea to truncate *before* computing the residual
     // but simple tests suggest otherwise ... no more iterations and
     // reduced iteration time from truncating.
-    START_TIMER(world);
-    truncate(world, new_psi);
-    END_TIMER(world, "Truncate new psi");
+    //START_TIMER(world);
+    //truncate(world, new_psi);
+    //END_TIMER(world, "Truncate new psi");
 
     START_TIMER(world);
     vecfuncT r = sub(world, psi, new_psi);
diff --git a/src/madness/chem/exchangeoperator.h b/src/madness/chem/exchangeoperator.h
index 563e3f1c9c4..9be188eb3d9 100644
--- a/src/madness/chem/exchangeoperator.h
+++ b/src/madness/chem/exchangeoperator.h
@@ -412,10 +412,10 @@ class Exchange<T,NDIM>::ExchangeImpl {
                 //}
                 Kf[0] += tmp_Kf;
                 print_size(world, Kf, "Kf before truncation");
-                truncate(world, Kf, mul_tol);
+                truncate(world, Kf);
                 print_size(world, Kf, "Kf after truncation");
             }
-            Kf[0] +=dot(world, mo_ket, psif);
+            //Kf[0] +=dot(world, mo_ket, psif);
             //if (world.rank()==0 && printdebug()) Kf[0].print_size("Kf after dot");
             
             //truncate(world, Kf);