Davidson works

CI/CI.ml

@@ -107,7 +107,7 @@ let make det_space =
     |> Farm.run ~ordered:false
       ~f:(fun i ->
           let ki = det.(i) in
-          Printf.eprintf "%d / %d\n%!" i ndet;
+          Printf.eprintf "%8d / %8d\r%!" i ndet;
           List.init ndet (fun j ->
             let kj = det.(j) in
             let x = h_ij mo_basis ki kj in
@@ -130,9 +130,15 @@ let make det_space =
   in
   let eigensystem = lazy (
     let m_H = Lazy.force m_H in
+    (*
     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
   )
   in
   { det_space ; m_H ; m_S2 ; eigensystem }
 
+

Utils/Constants.ml

@@ -3,6 +3,7 @@ let epsilon = 1.e-20
 
 (** Constants *)
 let pi = acos (-1.)
+let two_pi = 2. *. pi
 let sq_pi = sqrt pi
 let sq_pi_over_two = sq_pi *. 0.5
 let pi_inv = 1. /. pi

Utils/Constants.mli

@@ -14,6 +14,9 @@ val integrals_cutoff : float
 val pi : float
 (** {% $\pi$ %} = 3.141_592_653_589_793_12 *)
 
+val two_pi : float
+(** {% $2\pi$ %} *)
+
 val sq_pi : float
 (** {% $\sqrt{\pi}$ %} *)
 

Utils/Davidson.ml

@@ -0,0 +1,133 @@
+open Lacaml.D
+
+type t
+
+let make
+    ?guess
+    ?(n_states=8)
+    ?(n_iter=10)
+    ?(threshold=1.e-10)
+    diagonal
+    matrix_vector
+  =
+
+  let n = Vec.dim diagonal in  (* Size of the matrix to diagonalize *)
+
+
+  let m = (* Number of requested states *)
+    match guess with
+    | Some vectors -> (Mat.dim2 vectors) * n_states
+    | None -> n_states
+  in
+
+  (* Create guess vectors u, with randomly initialized unknown vectors. *)
+  let random_vectors =
+    let random_vector k =
+      Vec.init n (fun i ->
+          let r1 = Random.float 1.
+          and r2 = Random.float 1.
+          in
+          let a = sqrt (-2. *. log r1)
+          and b = Constants.two_pi *. r2
+          in
+          if i<k then 0.
+          else if i>k then
+            a *. cos b
+          else 100.0
+        )
+      |> Util.normalize
+    in
+    List.init m (fun i -> random_vector i) 
+  in
+
+  let pick_new u =
+    Mat.to_col_vecs_list u
+    |> Util.list_pack m
+    |> List.rev
+    |> List.hd
+  in
+
+  let u_new =
+   match guess with
+   | Some vectors -> Mat.to_col_vecs_list vectors 
+   | None -> random_vectors
+  in
+
+  let rec iteration u u_new w iter =
+    (* u is a list of orthonormal vectors, on which the operator has
+       been applied : w = op.u
+       u_new is a list of vector which will increase the size of the
+       space.
+    *)
+
+    (* Orthonormalize input vectors u_new *)
+    let u_new_ortho =
+      List.concat [u ; u_new]
+      |> Mat.of_col_vecs_list
+      |> Util.qr_ortho 
+      |> pick_new
+    in
+
+    (* Apply the operator the m last vectors *)
+    let w_new =
+      List.map matrix_vector u_new_ortho
+    in             
+
+    (* Data for the next iteration *)
+    let u_next =
+      List.concat [ u ; u_new_ortho ]
+    and w_next = 
+      List.concat [ w ; w_new ]
+    in
+
+    (* Build the small matrix h = <U_k | W_l> *)
+    let m_U =
+      Mat.of_col_vecs_list u_next
+    and m_W =
+      Mat.of_col_vecs_list w_next
+    in
+    let m_h =
+      gemm ~transa:`T m_U m_W
+    in
+
+    (* Diagonalize h *)
+    let y, lambda =
+      Util.diagonalize_symm m_h
+    in
+
+    (* Express m lowest eigenvectors of h in the large basis *)
+    let m_new_U =
+      gemm  ~n:m  m_U  y
+    and m_new_W =
+      gemm  ~n:m  m_W  y
+    in
+
+    (* Compute the residual as proposed new vectors *)
+    let u_proposed = 
+      Mat.init_cols n m (fun i k -> (lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. 
+                              (max (diagonal.{i} -. lambda.{k}) 0.01) )
+      |> Mat.to_col_vecs_list
+    in
+
+    let residual_norms = List.map nrm2 u_proposed in
+    let residual_norm  = List.fold_left (fun accu i -> max accu i) 0. residual_norms in
+
+    Printf.printf "%3d %16.10f  %16.8e%!\n" iter lambda.{1} residual_norm;
+
+    if residual_norm > threshold then
+      let u_next, w_next, iter =
+        if iter = n_iter then
+          m_new_U |> pick_new, 
+          m_new_W |> pick_new,
+          0
+        else
+          u_next, w_next, iter
+      in
+      iteration u_next u_proposed w_next (iter+1)
+    else
+      (Mat.of_col_vecs_list u_next |> pick_new |> Mat.of_col_vecs_list), lambda
+
+  in
+  iteration [] u_new [] 1
+
+

Utils/Util.ml

@@ -178,6 +178,29 @@ let list_some l =
   |> List.map (function Some x -> x | _ -> assert false)
 
 
+let list_range first last = 
+  if last < first then [] else
+  let rec aux accu = function
+  | 0 -> first :: accu
+  | i -> aux ( (first+i)::accu ) (i-1)
+  in
+  aux [] (last-first)
+
+
+let list_pack n l =
+  let rec aux i accu1 accu2 = function
+  | [] -> if accu1 = [] then
+             List.rev accu2
+          else 
+             List.rev ((List.rev accu1) :: accu2)
+  | a :: rest ->
+    match i with
+    | 0 -> aux (n-1)  [] ((List.rev (a::accu1)) :: accu2) rest
+    | _ -> aux (i-1)  (a::accu1) accu2 rest
+  in
+  aux (n-1) [] [] l
+
+
 (** {2 Stream functions} *)
 
 let stream_range first last = 
@@ -195,16 +218,6 @@ let stream_to_list stream =
   in aux []
 
 
-let list_range first last = 
-  if last < first then [] else
-  let rec aux accu = function
-  | 0 -> first :: accu
-  | i -> aux ( (first+i)::accu ) (i-1)
-  in
-  aux [] (last-first)
-
-
-
 
 (** {2 Linear algebra} *)
 
@@ -250,6 +263,27 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c =
   gemm c u
 
 
+let qr_ortho m =
+  (** Performed twice for precision *)
+  let result = lacpy m in
+  let tau = geqrf result in
+  orgqr ~tau result;
+  let tau = geqrf result in
+  orgqr ~tau result;
+  result
+
+
+let normalize v =
+  let result = copy v in
+  scal (1. /. (nrm2 v)) result;
+  result
+
+
+let normalize_mat m =
+  Mat.to_col_vecs m 
+  |> Array.map (fun v -> normalize v) 
+  |> Mat.of_col_vecs
+
 
 (** {2 Bitstring functions} *)
 let bit_permtutations m n =

Utils/Util.mli

@@ -1,3 +1,5 @@
+open Lacaml.D
+
 (** All utilities which should be included in all source files are defined here *)
 
 (** {2 Functions from libm} *)
@@ -69,6 +71,15 @@ val list_some : 'a option list -> 'a list
 val list_range : int -> int -> int list
 (** [list_range first last] returns a list [first; first+1 ; ... ; last-1 ; last ]. *)
 
+val list_pack : int -> 'a list -> 'a list list
+(** Example:
+{[
+    list_pack 3 [ 1; 2; 3; ...; 18; 19; 20 ] =
+      [[1; 2; 3]; [4; 5; 6]; [7; 8; 9]; [10; 11; 12]; [13; 14; 15];
+      [16; 17; 18]; [19; 20]]
+]}
+*)
+
 (** {2 Useful streams} *)
 val stream_range : int -> int -> int Stream.t
 (** [stream_range first last] returns a stream <first ; first+1 ; ... ; last-1 ; last>. *)
@@ -78,16 +89,25 @@ val stream_to_list : 'a Stream.t -> 'a list
 
 (** {2 Linear algebra } *)
 
-val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec
+val normalize : Vec.t -> Vec.t
+(** Returns a the vector normalized *)
+
+val normalize_mat : Mat.t -> Mat.t
+(** Returns the matrix with all the column vectors normalized *)
+
+val diagonalize_symm : Mat.t -> Mat.t * Vec.t
 (** Diagonalize a symmetric matrix. Returns the eigenvectors and the eigenvalues. *)
 
-val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
+val xt_o_x : o:Mat.t -> x:Mat.t -> Mat.t
 (** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *)
 
-val x_o_xt : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
+val x_o_xt : o:Mat.t -> x:Mat.t -> Mat.t
 (** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *)
 
-val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat
+val qr_ortho : Mat.t -> Mat.t
+(** QR orthogonalization of the input matrix *)
+
+val canonical_ortho: ?thresh:float -> overlap:Mat.t -> Mat.t -> Mat.t
 (** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %},
     and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize.
 
@@ -100,14 +120,14 @@ $$
     
     *)
 
-val debug_matrix: string -> Lacaml.D.mat -> unit
+val debug_matrix: string -> Mat.t -> unit
 (** Prints a matrix in stdout for debug *)
 
-val matrix_of_file : string -> Lacaml.D.mat
+val matrix_of_file : string -> Mat.t
 (** Reads a matrix from a file with format "%d %d %f" corresponding to
     [i, j, A.{i,j}]. *)
 
-val sym_matrix_of_file : string -> Lacaml.D.mat
+val sym_matrix_of_file : string -> Mat.t
 (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to
     [i, j, A.{i,j}]. *)
 
@@ -164,5 +184,5 @@ val pp_float_2darray : Format.formatter -> float array array -> unit
 val pp_bitstring : int -> Format.formatter -> Z.t -> unit
 (** Example: [ pp_bitstring 14 ppf z -> +++++------+-- ] *)
 
-val pp_matrix : Format.formatter -> Lacaml.D.mat -> unit
+val pp_matrix : Format.formatter -> Mat.t -> unit