From a0b84249ec5c85a82ad20b3a572314d558d1fb56 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 5 Apr 2019 09:46:23 +0200 Subject: [PATCH] Accelerated direct FCI --- Basis/TwoElectronIntegralsNonSeparable.ml | 2 +- Basis/TwoElectronIntegralsSeparable.ml | 2 +- CI/CI.ml | 57 ++++++++++++----------- CI/Spindeterminant.ml | 3 ++ CI/Spindeterminant.mli | 3 ++ Utils/FourIdxStorage.ml | 2 +- Utils/Matrix.ml | 13 ++++-- 7 files changed, 48 insertions(+), 34 deletions(-) diff --git a/Basis/TwoElectronIntegralsNonSeparable.ml b/Basis/TwoElectronIntegralsNonSeparable.ml index 5d73675..83995ae 100644 --- a/Basis/TwoElectronIntegralsNonSeparable.ml +++ b/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); diff --git a/Basis/TwoElectronIntegralsSeparable.ml b/Basis/TwoElectronIntegralsSeparable.ml index b5fe2b8..a3b9156 100644 --- a/Basis/TwoElectronIntegralsSeparable.ml +++ b/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); diff --git a/CI/CI.ml b/CI/CI.ml index 1cf7487..c192bec 100644 --- a/CI/CI.ml +++ b/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 diff --git a/CI/Spindeterminant.ml b/CI/Spindeterminant.ml index 2a254b9..3bce41e 100644 --- a/CI/Spindeterminant.ml +++ b/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 diff --git a/CI/Spindeterminant.mli b/CI/Spindeterminant.mli index d87df4f..48651ef 100644 --- a/CI/Spindeterminant.mli +++ b/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 diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 9ba2996..ed38b84 100644 --- a/Utils/FourIdxStorage.ml +++ b/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) -> diff --git a/Utils/Matrix.ml b/Utils/Matrix.ml index 0fbb228..b11dec5 100644 --- a/Utils/Matrix.ml +++ b/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