From 37ccd24d3ed4a54b2852dd54a1eefca2efa458b7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 5 Apr 2019 16:54:38 +0200 Subject: [PATCH] Accelerated 4idx --- Utils/FourIdxStorage.ml | 163 +++++++++++++++++++++------------------ Utils/FourIdxStorage.mli | 16 ++-- Utils/Util.ml | 4 +- run_cas.ml | 3 + 4 files changed, 99 insertions(+), 87 deletions(-) diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 191f702..8a8d016 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -102,78 +102,75 @@ 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 - if i=k then - if j=l then - begin - unsafe_set t.two_index i j value; - unsafe_set t.two_index j i value; - end - else + 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; 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 - else if j=l then - begin - unsafe_set t.three_index (dense_index i k t.size) j value; - unsafe_set t.three_index (dense_index k i t.size) j value; - end + else if j=l then + begin + unsafe_set t.three_index (dense_index i k t.size) j value; + unsafe_set t.three_index (dense_index k i t.size) j value; + end - else if i=l then - if j=k then - begin - unsafe_set t.two_index_anti i j value; - unsafe_set t.two_index_anti j i value; - end - else + 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; 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 - else if j=k then - begin - unsafe_set t.three_index_anti (dense_index i l t.size) j value; - unsafe_set t.three_index_anti (dense_index l i t.size) j value; - end + else if j=k then + begin + unsafe_set t.three_index_anti (dense_index i l t.size) j value; + unsafe_set t.three_index_anti (dense_index l i t.size) j value; + end - else if i=j then - if k=l then - begin - unsafe_set t.two_index_anti i k value; - unsafe_set t.two_index_anti k i value; - end - else - (* *) + 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; 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 - else if k=l then - (* *) - begin - 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 + else if k=l then + (* *) + begin + 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) - and jl = (dense_index j l t.size) - and ki = (dense_index k i t.size) - and lj = (dense_index l j t.size) - and ik_s = (sym_index i k) - and jl_s = (sym_index j l) - in - begin - unsafe_set a ik jl_s value; - unsafe_set a ki jl_s value; - unsafe_set a jl ik_s value; - unsafe_set a lj ik_s value; - end - | Sparse a -> let key = key_of_indices ~r1 ~r2 in - Hashtbl.replace a key value + match t.four_index with + | Dense a -> let ik = (dense_index i k t.size) + and jl = (dense_index j l t.size) + and ki = (dense_index k i t.size) + and lj = (dense_index l j t.size) + and ik_s = (sym_index i k) + and jl_s = (sym_index j l) + in + begin + unsafe_set a ik jl_s value; + unsafe_set a ki jl_s value; + unsafe_set a jl ik_s value; + unsafe_set a lj ik_s value; + end + | Sparse a -> let key = key_of_indices ~r1 ~r2 in + Hashtbl.replace a key value let set_four_index ~r1 ~r2 ~value t = @@ -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 = - Array.init d.size (fun j -> get_phys_all_i d ~j:(j+1) ~k ~l) +let get_phys_all_ij d ~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 = - Array.init d.size (fun j -> get_chem_all_i d ~j:(j+1) ~k ~l) +let get_chem_all_ij d ~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 -> - incr jk; - get_chem_all_i source ~j ~k ~l - |> Array.iteri (fun i x -> o.{i+1,!jk} <- x) - (* - lacpy ~bc:!jk ~b:o - (Mat.of_col_vecs [| Vec.of_array (get_chem_all_i source ~j ~k ~l) |] ) - |> ignore - *) - ) range_ao + get_chem_all_ij source ~k ~l + |> lacpy ~b:o ~bc:!jk + |> ignore; + jk := !jk + ao_num; ) range_ao; (* o_i_jk *) diff --git a/Utils/FourIdxStorage.mli b/Utils/FourIdxStorage.mli index 2df9ccd..9bddd54 100644 --- a/Utils/FourIdxStorage.mli +++ b/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 -(** Get all integrals in an array [a.(i-1) =] {% $(\cdot j|kl)$ %} . *) +val get_chem_all_i : t -> j:int -> k:int -> l:int -> Vec.t +(** 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 -(** Get all integrals in an array [a.(i-1) =] {% $\langle \cdot j|kl \rangle$ %} . *) +val get_phys_all_i : t -> j:int -> k:int -> l:int -> Vec.t +(** 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 -(** Get all integrals in an array [a.(j-1).(i-1) =] {% $(\cdot \cdot|kl)$ %} . *) +val get_chem_all_ij : t -> k:int -> l:int -> Mat.t +(** 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 -(** Get all integrals in an array [a.(j-1).(i-1) =] {% $\langle \cdot \cdot|kl \rangle$ %} . *) +val get_phys_all_ij : t -> k:int -> l:int -> Mat.t +(** 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. *) diff --git a/Utils/Util.ml b/Utils/Util.ml index 591c70d..aae6faf 100644 --- a/Utils/Util.ml +++ b/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) ;; diff --git a/run_cas.ml b/run_cas.ml index 89a91e2..436ff94 100644 --- a/run_cas.ml +++ b/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