From 346db83439b8f8ca04f7f656a7e3c8f0bc2885c3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 10 Feb 2018 22:19:13 +0100 Subject: [PATCH] Less allocations --- Basis/TwoElectronRRVectorized.ml | 153 ++++++++++++++++++------------- 1 file changed, 89 insertions(+), 64 deletions(-) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 375e5dc..f552385 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -24,9 +24,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) map_1d map_2d np nq = - let empty = - Array.make nq 0. - in + let empty = Array.make nq 0. in + let tmp_1 = Array.make nq 0. in let totAngMom_a = Angular_momentum.to_int totAngMom_a @@ -54,59 +53,78 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) if amxyz < 0 then None else - let v1_top, p1_top = - if abs_float (Coordinate.coord center_ab xyz) < cutoff then - None, - vrr0_v (m+1) am (totAngMom_a-1) - else - vrr0_v m am (totAngMom_a-1), - vrr0_v (m+1) am (totAngMom_a-1) - in - let v1_top2, p1_top2 = - if amxyz < 1 then (None,None) else - vrr0_v m amm (totAngMom_a-2), - vrr0_v (m+1) amm (totAngMom_a-2) - in + begin + let cab = Coordinate.coord center_ab xyz in + let v1_top, p1_top = + if abs_float cab < cutoff then + None, + vrr0_v (m+1) am (totAngMom_a-1) + else + vrr0_v m am (totAngMom_a-1), + vrr0_v (m+1) am (totAngMom_a-1) + in + let v1_top2, p1_top2 = + if amxyz < 1 then (None,None) else + vrr0_v m amm (totAngMom_a-2), + vrr0_v (m+1) amm (totAngMom_a-2) + in - Some ( - Array.init np (fun l -> - let v1 = - let f = - -. expo_b.(l) *. expo_inv_p.(l) *. (Coordinate.coord center_ab xyz) - in - match v1_top with - | Some v1_top -> - v1_top.(l) |> Array.map (fun x -> f *. x) - | None -> empty - in - let p1 = - match p1_top with - | Some p1_top -> p1_top.(l) - | _ -> assert false - in - let p1 = - Array.init nq (fun k -> - v1.(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p1.(k)) - in - if amxyz < 1 then p1 else - let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5 - in - let v1 = - match v1_top2 with - | Some v1_top2 -> v1_top2.(l) - | None -> assert false - in - let v2 = - match p1_top2 with - | Some p1_top2 -> p1_top2.(l) - | None -> assert false - in - Array.init nq (fun k -> - p1.(k) +. f *. (v1.(k) +. v2.(k) *. expo_inv_p.(l) ) ) - ) - ) - in Zmap.add map_1d.(m) key result; - result + let result = Array.make_matrix np nq 0. in + for l=0 to np-1 + do + let p0 = + match p1_top with + | Some p1_top -> p1_top.(l) + | _ -> assert false + in + let f0 = + -. expo_b.(l) *. expo_inv_p.(l) *. cab + in + let v0 = + match v1_top with + | Some v1_top -> v1_top.(l) + | None -> empty + in + for k=0 to nq-1 + do + tmp_1.(k) <- v0.(k) *. f0 + done; + let v0 = tmp_1 in + if amxyz < 1 then + begin + for k=0 to nq-1 + do + result.(l).(k) <- v0.(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(k) + done + end + else + begin + let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5 + in + let v1 = + match v1_top2 with + | Some v1_top2 -> v1_top2.(l) + | None -> assert false + in + let v2 = + match p1_top2 with + | Some p1_top2 -> p1_top2.(l) + | None -> assert false + in + for k=0 to nq-1 + do + result.(l).(k) <- + v0.(k) + +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(k) + +. f *. (v1.(k) +. v2.(k) *. expo_inv_p.(l) ) + done + end + done; + Some result + end + in + Zmap.add map_1d.(m) key result; + result and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = @@ -140,6 +158,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (acx-2,acy,acz), aax,acx, 0 in + (* + let result = Array.make_matrix np nq 0. in + *) let f1 = let f = (Coordinate.coord center_cd xyz) in Array.init nq (fun k -> expo_d.(k) *. expo_inv_q.(k) *. f) @@ -484,19 +505,23 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.init (maxm+1) (fun _ -> Array.init np (fun _ -> Array.make nq 0. ) ) in + let empty = Array.make (maxm+1) 0. in Array.iteri (fun ab shell_ab -> let zero_m_array_tmp = Array.mapi (fun cd shell_cd -> - let expo_pq_inv = - expo_inv_p.(ab) +. expo_inv_q.(cd) - in - let norm_pq_sq = - center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +. - center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +. - center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd) - in + if (abs_float coef.(ab).(cd) < cutoff) then + empty + else + let expo_pq_inv = + expo_inv_p.(ab) +. expo_inv_q.(cd) + in + let norm_pq_sq = + center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +. + center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +. + center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd) + in - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq ) sq (* |> Array.to_list