Accelerated 4idx

Utils/FourIdxStorage.ml

@@ -102,14 +102,15 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
   let open Bigarray.Array2 in
   let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in 
 
+  let () = 
+
     if i=k then
+      begin
         if j=l then
           begin
             unsafe_set t.two_index i j value;
             unsafe_set t.two_index j i value;
-      end
+          end;
-    else
-      begin
         unsafe_set t.three_index (dense_index j l t.size) i value;
         unsafe_set t.three_index (dense_index l j t.size) i value;
       end
@@ -120,13 +121,12 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
         end
 
     else if i=l then
+      begin
         if j=k then
           begin
             unsafe_set t.two_index_anti i j value;
             unsafe_set t.two_index_anti j i value;
-      end
+          end;
-    else
-      begin
         unsafe_set t.three_index_anti (dense_index j k t.size) i value;
         unsafe_set t.three_index_anti (dense_index k j t.size) i value;
       end
@@ -137,14 +137,12 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
         end
 
     else if i=j then 
+      begin
         if k=l then
           begin
             unsafe_set t.two_index_anti i k value;
             unsafe_set t.two_index_anti k i value;
-      end
+          end;
-    else 
-      (* <ii|kl> *)
-      begin
         unsafe_set t.three_index_anti (dense_index k l t.size) i value;
         unsafe_set t.three_index_anti (dense_index l k t.size) i value;
       end
@@ -155,8 +153,7 @@ let unsafe_set_four_index ~r1 ~r2 ~value t =
           unsafe_set t.three_index_anti (dense_index i j t.size) k value;
           unsafe_set t.three_index_anti (dense_index j i t.size) k value;
         end
-
+  in
-  else
 
   match t.four_index with
   | Dense  a -> let ik = (dense_index i k t.size)
@@ -268,19 +265,37 @@ type element = (** Element for the stream *)
 
 
 let get_phys_all_i d ~j ~k ~l  =
-  Array.init d.size (fun i -> get_phys d (i+1) j k l)
+  Vec.init d.size (fun i -> get_phys d i j k l)
     
 
 let get_chem_all_i d ~j ~k ~l  =
-  Array.init d.size (fun i -> get_chem d (i+1) j k l)
+  Vec.init d.size (fun i -> get_chem d i j k l)
 
 
-let get_phys_all_ji d ~k ~l  =
+let get_phys_all_ij d ~k ~l  =
-  Array.init d.size (fun j -> get_phys_all_i d ~j:(j+1) ~k ~l)
+  Mat.init_cols d.size d.size (fun i j -> get_phys d i j k l)
     
 
-let get_chem_all_ji d ~k ~l  =
+let get_chem_all_ij d ~k ~l  =
-  Array.init d.size (fun j -> get_chem_all_i d ~j:(j+1) ~k ~l)
+  (*
+  if k = l then
+      let result = 
+        Mat.col d.three_index k
+        |> Bigarray.genarray_of_array1 
+      in
+      Bigarray.reshape_2 result d.size d.size
+  else
+    *)
+  match d.four_index with
+  | Dense  a -> 
+      let kl = sym_index k l in
+      let result = 
+        Mat.col a kl
+        |> Bigarray.genarray_of_array1
+      in
+      Bigarray.reshape_2 result d.size d.size
+  | Sparse a ->
+      Mat.init_cols d.size d.size (fun i j -> get_chem d i j k l)
     
 
     
@@ -292,16 +307,16 @@ let to_stream d =
   and l = ref 1
   in
   let rec f_dense _ =
-    i := !i+1;
+    incr i;
     if !i > !k then begin
       i := 1;
-      j := !j + 1;
+      incr j;
       if !j > !l then begin
         j := 1;
-        k := !k + 1;
+        incr k;
         if !k > !l then begin
           k := 1;
-          l := !l + 1;
+          incr l;
         end;
       end;
      end;
@@ -430,18 +445,12 @@ let four_index_transform coef source =
       List.iter (fun l ->
           if abs_float coef.{l,delta} > epsilon then
             begin
-              let jk = ref 0 in
+              let jk = ref 1 in
               List.iter (fun k ->
-                  List.iter (fun j ->
+                  get_chem_all_ij source ~k ~l 
-                      incr jk;
+                  |> lacpy ~b:o ~bc:!jk 
-                      get_chem_all_i source ~j ~k ~l
+                  |> ignore;
-                      |> Array.iteri (fun i x -> o.{i+1,!jk} <- x)
+                  jk := !jk + ao_num;
-                      (*
-                      lacpy ~bc:!jk ~b:o
-                      (Mat.of_col_vecs [| Vec.of_array (get_chem_all_i source ~j ~k ~l) |] ) 
-                      |> ignore
-                      *)
-                    ) range_ao
                 ) range_ao;
               (* o_i_jk *)
 

Utils/FourIdxStorage.mli

@@ -38,17 +38,17 @@ val set_chem : t -> int -> int -> int -> int -> float -> unit
 val set_phys : t -> int -> int -> int -> int -> float -> unit
 (** Set an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *)
 
-val get_chem_all_i : t -> j:int -> k:int -> l:int -> float array
+val get_chem_all_i : t -> j:int -> k:int -> l:int -> Vec.t
-(** Get all integrals in an array [a.(i-1) =] {% $(\cdot j|kl)$ %} . *)
+(** Get all integrals in an array [a.{i} =] {% $(\cdot j|kl)$ %} . *)
 
-val get_phys_all_i : t -> j:int -> k:int -> l:int -> float array
+val get_phys_all_i : t -> j:int -> k:int -> l:int -> Vec.t
-(** Get all integrals in an array [a.(i-1) =] {% $\langle \cdot j|kl \rangle$ %} . *)
+(** Get all integrals in an array [a.{i} =] {% $\langle \cdot j|kl \rangle$ %} . *)
 
-val get_chem_all_ji : t -> k:int -> l:int -> float array array
+val get_chem_all_ij : t -> k:int -> l:int -> Mat.t
-(** Get all integrals in an array [a.(j-1).(i-1) =] {% $(\cdot \cdot|kl)$ %} . *)
+(** Get all integrals in an array [a.{i,j} =] {% $(\cdot \cdot|kl)$ %} . *)
 
-val get_phys_all_ji : t -> k:int -> l:int -> float array array
+val get_phys_all_ij : t -> k:int -> l:int -> Mat.t
-(** Get all integrals in an array [a.(j-1).(i-1) =] {% $\langle \cdot \cdot|kl \rangle$ %} . *)
+(** Get all integrals in an array [a.{i,j} =] {% $\langle \cdot \cdot|kl \rangle$ %} . *)
 
 val to_stream : t -> element Stream.t
 (** Retrun the data structure as a stream. *)

Utils/Util.ml

@@ -35,10 +35,10 @@ external leadz : int64 -> int32 = "leadz_bytecode" "leadz"
 let leadz i = leadz i |> Int32.to_int 
 
 
-exception Ctrl_C
+exception SIGTERM
 
 let () = 
-  let f _ = raise Ctrl_C in
+  let f _ = raise SIGTERM in
   Sys.set_signal Sys.sigint (Sys.Signal_handle f) 
   ;;
 

run_cas.ml

@@ -77,11 +77,14 @@ let () =
   Format.fprintf ppf "CAS-MP2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
 *)
 
+(*
+
   let pt2 = CI.pt2_en ci in
   Format.fprintf ppf "CAS-EN2 energy             : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
 
   let pt2 = CI.pt2_en_reference ci in
   Format.fprintf ppf "CAS-EN2 energy (reference) : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
+  *)
 
 (*
   let variance = CI.variance ci in