From a976bacfdf9513eccaf62f0432b6638e2b2514f3 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 28 Feb 2019 12:30:20 +0100
Subject: [PATCH] Davidson with sparse matrix

---
 CI/CI.ml          | 42 ++++++++++++++++++++++++++----------------
 Utils/Davidson.ml |  9 +++++++--
 Utils/Matrix.ml   |  4 ++++
 Utils/Matrix.mli  |  2 +-
 4 files changed, 38 insertions(+), 19 deletions(-)

diff --git a/CI/CI.ml b/CI/CI.ml
index 2fae19a..fd8e50e 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -5,8 +5,8 @@ module Ds = Determinant_space
 type t =
   {
     det_space   : Ds.t ;
-    m_H         : Mat.t lazy_t ;
-    m_S2        : Mat.t lazy_t ;
+    m_H         : Matrix.t lazy_t ;
+    m_S2        : Matrix.t lazy_t ;
     eigensystem : (Mat.t * Vec.t) lazy_t;
   }
 
@@ -53,18 +53,19 @@ let make det_space =
   let ndet = Ds.size det_space in
   let det  = Ds.determinants det_space in
   let mo_basis = Ds.mo_basis det_space in
-  (*
+
   let m_H = lazy (
-    Util.list_range 0 (ndet-1)
-    |> List.map (fun i -> let ki = det.(i) in
-                  Array.init ndet (fun j -> let kj = det.(j) in
+    Array.init ndet (fun i -> let ki = det.(i) in
+                  Vec.init ndet (fun j -> let kj = det.(j-1) in
                                     h_ij mo_basis ki kj
                                   )
-                  |> Vec.of_array)
-    |> Mat.of_col_vecs_list
+                  |> Vector.sparse_of_vec
+                  )
+    |> Matrix.sparse_of_vector_array
   )
   in
-  *)
+  
+
   (*
   let m_H = lazy (
     let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in
@@ -93,6 +94,8 @@ let make det_space =
     |> Mat.of_col_vecs_list
   ) in
   *)
+
+  (* Parallel
   let m_H = lazy (
     let h =
       if Parallel.master then
@@ -116,16 +119,19 @@ let make det_space =
           |> Util.list_some
           )
     |> Util.stream_to_list
-    |> List.iter (fun l ->  if Parallel.master then
+    |> List.iter (fun l -> if Parallel.master then
       List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l);
     Parallel.broadcast (lazy h)
   ) in
+  *)
   let m_S2 = lazy (
     Array.init ndet (fun i -> let ki = det.(i) in
-                      Array.init ndet (fun j -> let kj = det.(j) in
+                      Vec.init ndet (fun j -> let kj = det.(j-1) in
                                         CIMatrixElement.make_s2 ki kj
-                                      ))
-    |> Mat.of_array
+                                      )
+                      |> Vector.sparse_of_vec
+                      )
+    |> Matrix.sparse_of_vector_array
   )
   in
   let eigensystem = lazy (
@@ -134,9 +140,13 @@ let make det_space =
     Parallel.broadcast @@
       lazy (Util.diagonalize_symm m_H)
       *)
-    let diagonal = Vec.init (Mat.dim1 m_H) (fun i -> m_H.{i,i}) in
-    let matrix_vector psi = symv m_H psi in
-    Davidson.make diagonal matrix_vector
+    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
+    Davidson.make diagonal matrix_prod
   )
   in
   { det_space ; m_H ; m_S2 ; eigensystem }
diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml
index 67fcd57..bfca313 100644
--- a/Utils/Davidson.ml
+++ b/Utils/Davidson.ml
@@ -8,7 +8,7 @@ let make
     ?(n_iter=10)
     ?(threshold=1.e-8)
     diagonal
-    matrix_vector
+    matrix_prod
   =
 
   let n = Vec.dim diagonal in  (* Size of the matrix to diagonalize *)
@@ -70,7 +70,12 @@ let make
 
     (* Apply the operator the m last vectors *)
     let w_new =
-      List.map matrix_vector u_new_ortho
+      matrix_prod (
+        u_new_ortho
+        |> Mat.of_col_vecs_list
+        |> Matrix.dense_of_mat )
+      |> Matrix.to_mat
+      |> Mat.to_col_vecs_list
     in             
 
     (* Data for the next iteration *)
diff --git a/Utils/Matrix.ml b/Utils/Matrix.ml
index 94e7dff..690405f 100644
--- a/Utils/Matrix.ml
+++ b/Utils/Matrix.ml
@@ -60,6 +60,10 @@ let dense_of_sparse = function
 
 let dense_of_mat m = Dense m
 
+let rec to_vector_array ?(threshold=epsilon) = function
+  | Sparse {m ; n ; v} -> v
+  | Dense m -> to_vector_array (sparse_of_dense ~threshold (Dense m))
+
 
 let sparse_of_mat ?(threshold=epsilon) m =
   dense_of_mat m
diff --git a/Utils/Matrix.mli b/Utils/Matrix.mli
index 9fbf2a3..0d37233 100644
--- a/Utils/Matrix.mli
+++ b/Utils/Matrix.mli
@@ -27,7 +27,7 @@ val get : t -> int -> int -> float
 val to_mat : t -> Mat.t
 (** Convert into a Lacaml Mat. *)
 
-val to_vector_array : t -> Vector.t array
+val to_vector_array : ?threshold:float -> t -> Vector.t array
 (** Convert the matrix into an array of column vectors. *)
 
 val sparse_of_dense : ?threshold:float -> t -> t