From bb677d4f3a034350cad195ba539a1c561bc3d137 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 28 Feb 2019 18:18:26 +0100
Subject: [PATCH] Accelerated CI

---
 CI/CI.ml              | 109 +++++++++++++++++++++++++++++++++---------
 CI/Spindeterminant.ml |   5 +-
 2 files changed, 90 insertions(+), 24 deletions(-)

diff --git a/CI/CI.ml b/CI/CI.ml
index 91c4e45..ef35da9 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -70,30 +70,95 @@ let make ?(n_states=1) det_space =
     in
 
     let m_H_spin a b = lazy (
-      let v = Vec.make0 ndet in
+      let n_alfa = Array.length a in
+      let n_beta = Array.length b in
+      let result = Array.init ndet (fun _ -> []) in
+
+      (** Update function when ki and kj are connected *)
+      let update i j ki kj =
+        let x = h_ij mo_basis ki kj in
+        if x <> 0. then
+          result.(i) <- (j, x) :: result.(i) ;
+      in
+
+      (** Array of (list of singles, list of doubles) in the beta spin *)
+      let degree_bb = 
+        Array.map (fun det_i ->
+          let deg = Spindeterminant.degree det_i in
+          let doubles =
+            Array.mapi (fun i det_j -> 
+              let d = deg det_j in
+              if d < 3 then
+                Some (i,d,det_j)
+              else
+                None
+              ) b
+            |> Array.to_list
+            |> Util.list_some
+          in
+          let singles = 
+            List.filter (fun (i,d,det_j) -> d < 2) doubles
+            |> List.map (fun (i,_,det_j) -> (i,det_j))
+          in
+          let doubles = 
+            List.map (fun (i,_,det_j) -> (i,det_j)) doubles
+          in
+          (singles, doubles)
+        ) b
+      in
+      let a = Array.to_list a
+      and b = Array.to_list b
+      in
       let i = ref 0 in
-      let result = Array.make ndet (Vector.sparse_of_vec @@ Vec.make0 1) in
-      Array.iteri (fun ia i_alfa ->
-          Array.iteri (fun ib i_beta ->
-              Printf.eprintf "%8d / %8d\r%!" !i ndet;
-              let ki =
-                Determinant.of_spindeterminants i_alfa i_beta
-              in
-              let j = ref 0 in
-              Array.iteri (fun ja j_alfa ->
-                  Array.iteri (fun jb j_beta ->
-                      let kj =
-                        Determinant.of_spindeterminants j_alfa j_beta
-                      in
-                      incr j;
-                      v.{!j} <- h_ij mo_basis ki kj
-                    ) b;
-                ) a;
-              result.(!i) <- Vector.sparse_of_vec v;
-              incr i;
-            ) b
+      List.iteri (fun ia i_alfa ->
+          Printf.eprintf "%8d / %8d\r%!" ia n_alfa;
+          let j = ref 1 in
+          let deg_a = Spindeterminant.degree i_alfa in
+          List.iter (fun j_alfa ->
+              let degree_a = deg_a j_alfa in
+              begin
+                match degree_a with
+                  | 2 ->
+                      let i' = ref !i in
+                      List.iteri (fun ib i_beta ->
+                          let ki = Determinant.of_spindeterminants i_alfa i_beta in
+                          let kj = Determinant.of_spindeterminants j_alfa i_beta in
+                          update !i' (ib + !j) ki kj;
+                          incr i';
+                        ) b;
+                  | 1 -> 
+                      let i' = ref !i in
+                      List.iteri (fun ib i_beta ->
+                          let ki = Determinant.of_spindeterminants i_alfa i_beta in
+                          let singles, _ = degree_bb.(ib) in
+                          List.iter (fun (j', j_beta) ->
+                              let kj = Determinant.of_spindeterminants j_alfa j_beta in
+                              update !i' (j' + !j) ki kj
+                            ) singles;
+                          incr i';
+                        ) b;
+                  | 0 ->
+                      let i' = ref !i in
+                      List.iteri (fun ib i_beta ->
+                          let ki = Determinant.of_spindeterminants i_alfa i_beta in
+                          let _singles, doubles = degree_bb.(ib) in
+                          List.iter (fun (j', j_beta) ->
+                                let kj = Determinant.of_spindeterminants j_alfa j_beta in
+                                update !i' (j' + !j) ki kj
+                            ) doubles;
+                          incr i';
+                        ) b;
+                  | _ -> ();
+              end;
+              j := !j + n_beta
+            ) a;
+          i := !i + n_beta
         ) a;
-      Matrix.sparse_of_vector_array result
+
+      Array.map (fun l ->
+        List.sort compare l
+        |> Vector.sparse_of_assoc_list ndet ) result
+      |> Matrix.sparse_of_vector_array
     )
     in
 
diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml
index bbcab53..429f333 100644
--- a/CI/Spindeterminant.ml
+++ b/CI/Spindeterminant.ml
@@ -80,8 +80,9 @@ let rec bits_to_list accu = function
          in
          bits_to_list newlist Z.(logand t (t-one))
 
-let degree t t' =
-  Z.hamdist (bitstring t) (bitstring t') / 2
+let degree t = 
+  let bt = bitstring t in
+  fun t' -> Z.hamdist bt (bitstring t') / 2
 
 let holes_of t t' =
   Z.logand (bitstring t) (Z.logxor (bitstring t) (bitstring t'))