From bdccd66b1abb1a8fb481576d3670a0eccd349bb1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 29 Jan 2018 13:39:59 +0100 Subject: [PATCH] Vectorizing --- Basis/TwoElectronRRVectorized.ml | 93 +++++++++++++++++++------------- 1 file changed, 57 insertions(+), 36 deletions(-) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 85f0a19..bf6dcd8 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -26,13 +26,11 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0 m angMom_a = function | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod - |> getk | 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2 in Array.mapi (fun k c -> c *. expo_inv_p *. ( (Coordinate.coord center_pq.(k) i) *. zero_m_array.(k).(m+1) -. expo_b *. (Coordinate.coord center_ab i) *. zero_m_array.(k).(m) ) ) coef_prod - |> getk | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three @@ -52,19 +50,28 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) in am.(xyz) <- am.(xyz) - 1; amm.(xyz) <- amm.(xyz) - 2; -if am.(xyz) < 0 then 0. else - Array.mapi (fun k _ -> - (-. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)) - *. (vrr0 m am (totAngMom_a-1) ) - +. (expo_inv_p *. (Coordinate.coord center_pq.(k) xyz)) - *.(vrr0 (m+1) am (totAngMom_a-1) ) - +. (if amm.(xyz) < 0 then 0. else - ((float_of_int am.(xyz)) *. expo_inv_p *. 0.5) - *. (vrr0 m amm (totAngMom_a-2) - +. expo_inv_p *. (vrr0 (m+1) amm (totAngMom_a-2) ) ) ) -) coef_prod - |> getk -) + if am.(xyz) < 0 then + Array.map (fun _ -> 0.) coef_prod + else + let v1 = + vrr0 m am (totAngMom_a-1) + and v2 = + vrr0 (m+1) am (totAngMom_a-1) + in + let p1 = + Array.mapi (fun k _ -> + (-. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)) *. v1.(k) + +. (expo_inv_p *. (Coordinate.coord center_pq.(k) xyz)) *. v2.(k) + ) coef_prod + in + if amm.(xyz) < 0 then p1 else + let v1 = vrr0 m amm (totAngMom_a-2) + and v2 = vrr0 (m+1) amm (totAngMom_a-2) + in + Array.mapi (fun k _ -> p1.(k) +. + (float_of_int am.(xyz)) *. expo_inv_p *. 0.5 + *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) coef_prod + ) in if not found then Zmap.add map.(m) key result; @@ -74,7 +81,6 @@ if am.(xyz) < 0 then 0. else match (totAngMom_a, totAngMom_c) with | (0,0) -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod - |> getk | (_,0) -> vrr0 m angMom_a totAngMom_a | (_,_) -> @@ -99,20 +105,28 @@ if am.(xyz) < 0 then 0. else am.(xyz) <- am.(xyz) - 1; cm.(xyz) <- cm.(xyz) - 1; cmm.(xyz) <- cmm.(xyz) - 2; -if cm.(xyz) < 0 then 0. else Array.mapi (fun k _ -> - (-. expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) ) - *.(vrr m angMom_a cm totAngMom_a (totAngMom_c-1) ) - -. (expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz)) - *.(vrr (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) ) +if cm.(xyz) < 0 then 0. else + let v1 = + vrr m angMom_a cm totAngMom_a (totAngMom_c-1) + and v2 = + vrr (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + in + -. v1.(k) *. expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) + -. v2.(k) *. (expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz)) +. (if cmm.(xyz) < 0 then 0. else - ((float_of_int cm.(xyz)) *. expo_inv_q.(k) *. 0.5 ) - *.(vrr m angMom_a cmm totAngMom_a (totAngMom_c-2) - +. expo_inv_q.(k) *. (vrr (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) ) ) ) + let v1 = + vrr m angMom_a cmm totAngMom_a (totAngMom_c-2) + and v2 = + vrr (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) + in + ((float_of_int cm.(xyz)) *. expo_inv_q.(k) *. 0.5 ) *.(v1.(k) +. expo_inv_q.(k) *. v2.(k) ) ) -. (if am.(xyz) lor cm.(xyz) < 0 then 0. else - ((float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q.(k) *. 0.5 ) - *.(vrr (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) ) )) coef_prod - |> getk + let v = + vrr (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) + in + (float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q.(k) *. 0.5 *.v.(k) ) + ) coef_prod ) in if not found then @@ -128,11 +142,17 @@ if cm.(xyz) < 0 then 0. else match totAngMom_b with | 0 -> vrr m angMom_a angMom_c totAngMom_a totAngMom_c + |> getk | 1 -> let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in ap.(xyz) <- ap.(xyz) + 1; - vrr m ap angMom_c (totAngMom_a+1) totAngMom_c - +. (Coordinate.coord center_ab xyz) *. (vrr m angMom_a angMom_c totAngMom_a totAngMom_c) + let v1 = + vrr m ap angMom_c (totAngMom_a+1) totAngMom_c + |> getk + and v2 = + vrr m angMom_a angMom_c totAngMom_a totAngMom_c + |> getk + in v1 +. v2 *. (Coordinate.coord center_ab xyz) | _ -> let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |] @@ -144,18 +164,19 @@ if cm.(xyz) < 0 then 0. else in ap.(xyz) <- ap.(xyz) + 1; bm.(xyz) <- bm.(xyz) - 1; -if (bm.(xyz) < 0) then 0. else - hrr0 m ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c - +. (Coordinate.coord center_ab xyz) - *.(hrr0 m angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) - totAngMom_c ) + if (bm.(xyz) < 0) then 0. else + let h1 = + hrr0 m ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c + and h2 = + hrr0 m angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + in h1 +. h2 *. (Coordinate.coord center_ab xyz) and hrr m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = match (totAngMom_b, totAngMom_d) with - | (0,0) -> - vrr m angMom_a angMom_c totAngMom_a totAngMom_c + | (0,0) -> vrr m angMom_a angMom_c totAngMom_a totAngMom_c + |> getk | (_,0) -> hrr0 m angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |]