Accelerated direct FCI

Basis/TwoElectronIntegralsNonSeparable.ml

@@ -283,7 +283,7 @@ module Make(Zero_m : Zero_mType) = struct
       else
         Fis.create ~size:n `Dense
     in
-    Farm.run ~ordered:false ~f input_stream
+    Farm.run ~ordered:true ~f input_stream
     |> Stream.iter (fun l ->
         Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
             set_chem eri_array i_c j_c k_c l_c value) l);

Basis/TwoElectronIntegralsSeparable.ml

@@ -307,7 +307,7 @@ let of_basis_parallel basis =
     else
       Fis.create ~size:0 `Dense
   in
-  Farm.run ~ordered:false ~f input_stream
+  Farm.run ~ordered:true ~f input_stream
   |> Stream.iter (fun l ->
       Array.iter (fun (i_c,j_c,k_c,l_c,value) ->
           set_chem eri_array i_c j_c k_c l_c value) l);

CI/CI.ml

@@ -358,7 +358,6 @@ let create_matrix_spin_computed f det_space =
       | _ -> assert false
     in
     let n_beta = Array.length b in
-    let n_alfa = Array.length a in
 
     let h i_alfa j_alfa =
         let deg_a = Spindeterminant.degree a.(i_alfa) a.(j_alfa) in
@@ -401,40 +400,42 @@ let create_matrix_spin_computed f det_space =
     let i_prev = ref (-10)
     and result = ref (fun _ -> 0.)
     in
+    let h123 = ref (fun _ -> 0.)  in
 
     let g i =
       if i <> !i_prev then
         begin
           i_prev := i;
           let i_a = (i-1)/n_beta in
-          let i_alfa = i_a + 1 in
-          let h1 = 
-            h (i_alfa-1) 
-          in
-          let i_beta = i - i_a*n_beta  in
-          let bi = (i_beta-1) in
-          let h123_prev = ref (fun _ -> 0.) in
-          let j_a = ref (-n_alfa) in
-          let j_alfa_prev = ref (-10) in
-          let j0 = ref (!j_a * n_beta) in
+          let h1 = h i_a in
+          let i_b = i - i_a*n_beta - 1  in
+          let j0 = ref max_int in
+          let j1 = ref min_int in
+          let j_a = ref 0 in
           result := fun j ->
-              if  j > !j0 + n_beta
-              ||  j < !j0 
-              then begin
-                j_a := (j-1)/n_beta;
-                j0 := !j_a * n_beta
-              end;
-              let j_alfa = !j_a + 1 in
-              let h123 = 
-                if j_alfa <> !j_alfa_prev then
-                  begin
-                    j_alfa_prev := j_alfa ;
-                    h123_prev := (h1 (j_alfa-1) bi)
-                  end;
-                !h123_prev
-              in 
-              let j_beta = j - !j0 in
-              h123 (j_beta-1)
+              (* The following test checks if !j0 < j < !j1 *)
+              if  (!j1 - j) lor (j - !j0) > 0 then
+                begin
+                  let j_b = j - !j0 - 1 in
+                  !h123 j_b
+                end
+              else
+                begin
+                  if j < !j1 + n_beta then
+                    begin
+                      incr j_a;
+                      j0 := !j1;
+                    end
+                  else
+                    begin
+                      j_a := (j-1)/n_beta;
+                      j0 := !j_a * n_beta;
+                    end;
+                  j1 := !j0 + n_beta;
+                  h123 := h1 !j_a i_b;
+                  let j_b = j - !j0 - 1 in
+                  !h123 j_b
+                end
         end;
       !result
     in

CI/Spindeterminant.ml

@@ -119,6 +119,9 @@ let rec to_list = function
       else List.rev accu
     in aux [] spindet.bitstring
 
+let rec to_array t = 
+  to_list t
+  |> Array.of_list
 
 let n_electrons = function
   | Some t -> Bitstring.popcount t.bitstring

CI/Spindeterminant.mli

@@ -76,6 +76,9 @@ val of_list : int -> int list -> t
 val to_list : t -> int list
 (** Transforms a spin-determinant into a list of orbital indices. *)
 
+val to_array : t -> int array
+(** Transforms a spin-determinant into an array of orbital indices. *)
+
 (** {1 Printers}. *)
 
 val pp_spindet : int -> Format.formatter -> t -> unit

Utils/FourIdxStorage.ml

@@ -386,7 +386,7 @@ let four_index_transform coef source =
 
   let n = ref 0 in
   Stream.of_list range_mo
-  |> Farm.run ~f:task ~ordered:false ~comm:Parallel.Node.comm
+  |> Farm.run ~f:task ~ordered:true ~comm:Parallel.Node.comm
   |> Stream.iter (fun l ->
      if Parallel.master then (incr n ; Printf.eprintf "\r%d / %d%!" !n mo_num);
      Array.iter (fun (alpha, beta, gamma, delta, x) ->

Utils/Matrix.ml

@@ -358,10 +358,17 @@ let rec mm ?(transa=`N) ?(transb=`N) ?(threshold=epsilon) a b =
             (Mat.to_col_vecs b).(j)
         in
         let accu = Vec.make0 m' in
-        Vec.iteri (fun k a ->
+        let v = Vec.make0 m' in
+        let bj = Vec.to_array bj in
+        Array.iteri (fun k a ->
             if a <> 0. then
-              Vec.iteri (fun i vi ->  accu.{i} <- vi +. (f' i k) *. a) accu
-                ) bj;
+              begin
+                for i = 1 to m' do
+                  Bigarray.Array1.unsafe_set v i  (f' i (k+1));
+                done;
+                axpy ~alpha:a v accu
+              end
+          ) bj;
         accu
       )
       |> Mat.of_col_vecs