From 3c3b2f14abf5024082b149fef02009d8f78dcb1f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 2 Apr 2019 11:54:04 +0200
Subject: [PATCH] Working on CI

---
 CI/CI.ml          | 73 +++++++++++++++++++++++++++++++++--------------
 Utils/Davidson.ml | 27 +++++-------------
 2 files changed, 59 insertions(+), 41 deletions(-)

diff --git a/CI/CI.ml b/CI/CI.ml
index d7476be..ab0b0b6 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -323,10 +323,7 @@ let create_matrix_spin f det_space =
     in
 
     let result =
-      if Parallel.master then
-        Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
-      else
-        Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
+      Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
     in
 
     List.mapi (fun i i_alfa -> i*n_beta, i_alfa) a
@@ -334,7 +331,7 @@ let create_matrix_spin f det_space =
     |> Farm.run ~ordered:false ~f:task
     |> Stream.iter (fun (k, r) ->
         Array.iteri (fun j r_j -> result.(k+j) <- r_j) r;
-        Printf.eprintf "%8d / %8d\r%!" (k+1) ndet;
+        Printf.eprintf "%8d / %8d\r%!" (k+Array.length r) ndet;
       ) ;
     Matrix.sparse_of_vector_array result
   )
@@ -343,8 +340,7 @@ let create_matrix_spin f det_space =
 
 
 
-
-let make ?(n_states=1) det_space = 
+let make ?(n_states=1) ?(algo=`Direct) det_space = 
 
   let mo_basis = Ds.mo_basis det_space in
 
@@ -372,7 +368,7 @@ let make ?(n_states=1) det_space =
          if ki <> kj then
             h_ij mo_basis ki kj
          else
-            h_ij mo_basis ki kj -. e_shift
+            h_ij mo_basis ki ki -. e_shift
       ) det_space 
   in
 
@@ -386,22 +382,57 @@ let make ?(n_states=1) det_space =
   in
 
   let eigensystem = lazy (
-    let m_H =
-      Lazy.force m_H
+
+    let eigensystem_incore () = 
+      let m_H =
+        Lazy.force m_H
+      in
+      let diagonal =
+        Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
+      in
+      let matrix_prod psi = 
+        Matrix.mm ~transa:`T m_H psi
+      in
+      let eigenvectors, eigenvalues = 
+        let result = lazy (
+          Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
+        ) in
+        Parallel.broadcast result
+      in
+      let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
+      eigenvectors, eigenvalues
     in
-    let diagonal =
-      Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
+
+    let eigensystem_direct () = 
+      eigensystem_incore ()
     in
-    let matrix_prod psi = 
-      Matrix.mm ~transa:`T m_H psi
+(*
+      let diagonal =
+        let stream = Ds.determinant_stream det_space in
+        Vec.init (Ds.size det_space) (fun _ -> 
+          let ki = Stream.next stream in
+          h_ij mo_basis ki ki -. e_shift)
+      in
+
+
+      let matrix_prod psi = 
+        (*TODO*)
+      in
+
+      let eigenvectors, eigenvalues = 
+        let result = 
+          Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
+        in
+        Parallel.broadcast (lazy result)
+      in
+      let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
+      eigenvectors, eigenvalues
     in
-    let eigenvectors, eigenvalues = 
-      Parallel.broadcast (lazy (
-      Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
-      ))
-    in
-    let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
-    eigenvectors, eigenvalues
+*)
+
+    match algo with
+    | `Direct -> eigensystem_direct ()
+    | `InCore -> eigensystem_incore ()
   )
   in
   { det_space ; e_shift ; m_H ; m_S2 ; eigensystem ; n_states }
diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml
index a6b313a..2f2a827 100644
--- a/Utils/Davidson.ml
+++ b/Utils/Davidson.ml
@@ -22,26 +22,13 @@ let make
   in
 
   (* Create guess vectors u, with randomly initialized unknown vectors. *)
-  let random_vectors =
-    let random_vector k =
-      Vec.init n (fun i ->
-          if i<k then 0.
-          else if i=k then 1.e5
-          else
-            0.
-            (*
-            let r1 = Random.float 1.
-            and r2 = Random.float 1.
-            in
-            let a = sqrt (-2. *. log r1)
-            and b = Constants.two_pi *. r2 in
-            let c = a *. cos b in
-            if abs_float c > 1.e-1 then c else 0.
-              *)
-        )
-      |> Util.normalize
+  let init_vectors =
+    let init_vector k =
+      Vector.sparse_of_assoc_list n [ (k,1.0) ]
     in
-    List.init m (fun i -> random_vector (i+1)) 
+    Array.init m (fun i -> init_vector (i+1)) 
+    |> Matrix.sparse_of_vector_array
+    |> Matrix.to_mat
   in
 
   let pick_new u =
@@ -54,7 +41,7 @@ let make
   let u_new =
    match guess with
    | Some vectors -> Mat.to_col_vecs_list vectors 
-   | None -> random_vectors
+   | None -> Mat.to_col_vecs_list init_vectors
   in
 
   let rec iteration u u_new w iter =