From 03380374a1c0edb650a7e5a85e7498e007e0a9ac Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 11 Feb 2018 01:05:30 +0100 Subject: [PATCH] Better vectorization --- Basis/TwoElectronRRVectorized.ml | 99 +++++++++++++++++++++----------- 1 file changed, 65 insertions(+), 34 deletions(-) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 928e807..dca9853 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -24,7 +24,6 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) map_1d map_2d np nq = - let totAngMom_a = Angular_momentum.to_int totAngMom_a and totAngMom_b = Angular_momentum.to_int totAngMom_b and totAngMom_c = Angular_momentum.to_int totAngMom_c @@ -82,15 +81,15 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab in - for k=0 to nq-1 - do + for k=0 to nq-1 do result.(l).(k) <- v0.(l).(k) *. f0 done done end; for l=0 to np-1 do for k=0 to nq-1 do - result.(l).(k) <- result.(l).(k) +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) + result.(l).(k) <- result.(l).(k) + +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) done done end @@ -106,11 +105,8 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | None -> () | Some v0 -> for l=0 to np-1 do - let f0 = - -. expo_b.(l) *. expo_inv_p.(l) *. cab - in - for k=0 to nq-1 - do + let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab in + for k=0 to nq-1 do result.(l).(k) <- v0.(l).(k) *. f0 done done @@ -129,9 +125,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let f = (float_of_int amxyz) *. expo_inv_p.(l) *. 0.5 in for k=0 to nq-1 do - result.(l).(k) <- result.(l).(k) - +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) - +. f *. (v1.(l).(k) +. v2.(l).(k) *. expo_inv_p.(l) ) + result.(l).(k) <- result.(l).(k) +. + expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) +. + f *. (v1.(l).(k) +. v2.(l).(k) *. expo_inv_p.(l)) done done end; @@ -179,12 +175,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let do_compute = ref false in let v1 = let f = -. (Coordinate.coord center_cd xyz) in - let f1 = - Array.init nq (fun k -> + + let f1 = + Array.init nq (fun k -> let x = expo_d.(k) *. expo_inv_q.(k) *. f in - if (!do_compute) then x - else (if abs_float x > cutoff then do_compute := true ; x) - ) + if ( (not !do_compute) && (abs_float x > cutoff) ) then + do_compute := true; + x) in if (!do_compute) then match vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) with @@ -194,7 +191,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = Array.make_matrix np nq 0. in for l=0 to np-1 do for k=0 to nq-1 do -(* v1 rec *) result.(l).(k) <- v1.(l).(k) *. f1.(k) + result.(l).(k) <- v1.(l).(k) *. f1.(k) done done; Some result @@ -230,11 +227,21 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let p1 = match v1, v2 with | None, None -> None - | Some v1, Some v2 -> - Some (Array.init np (fun l -> Array.init nq (fun k -> - v1.(l).(k) +. v2.(l).(k) ) ) ) | None, Some v2 -> Some v2 | Some v1, None -> Some v1 + | Some v1, Some v2 -> + begin + for l=0 to np-1 do + for k=0 to nq-1 do + v2.(l).(k) <- v2.(l).(k) +. v1.(l).(k) + done + done; + Some v2 + end + (* + Some (Array.init np (fun l -> Array.init nq (fun k -> + v1.(l).(k) +. v2.(l).(k) ) ) ) + *) in let p2 = @@ -257,7 +264,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) for l=0 to np-1 do for k=0 to nq-1 do result.(l).(k) <- v1.(l).(k) *. f1.(k) - done + done; done; Some result end @@ -289,21 +296,45 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in match p1, v1, v2 with | None, None, None -> None - | Some p1, Some v1, Some v2 -> - Some (Array.init np (fun l -> Array.init nq (fun k -> - p1.(l).(k) +. v1.(l).(k) +. v2.(l).(k)) ) ) - | Some p1, Some v1, None -> - Some (Array.init np (fun l -> Array.init nq (fun k -> - p1.(l).(k) +. v1.(l).(k) ) ) ) - | Some p1, None, Some v2 -> - Some (Array.init np (fun l -> Array.init nq (fun k -> - p1.(l).(k) +. v2.(l).(k)) ) ) - | None , Some v1, Some v2 -> - Some (Array.init np (fun l -> Array.init nq (fun k -> - v1.(l).(k) +. v2.(l).(k)) ) ) | Some p1, None, None -> Some p1 | None, Some v1, None -> Some v1 | None, None, Some v2 -> Some v2 + | Some p1, Some v1, Some v2 -> + begin + for l=0 to np-1 do + for k=0 to nq-1 do + v2.(l).(k) <- p1.(l).(k) +. v1.(l).(k) +. v2.(l).(k) + done + done; + Some v2 + end + | Some p1, Some v1, None -> + begin + for l=0 to np-1 do + for k=0 to nq-1 do + p1.(l).(k) <- v1.(l).(k) +. p1.(l).(k) + done + done; + Some p1 + end + | Some p1, None, Some v2 -> + begin + for l=0 to np-1 do + for k=0 to nq-1 do + p1.(l).(k) <- p1.(l).(k) +. v2.(l).(k) + done + done; + Some p1 + end + | None , Some v1, Some v2 -> + begin + for l=0 to np-1 do + for k=0 to nq-1 do + v2.(l).(k) <- v1.(l).(k) +. v2.(l).(k) + done + done; + Some v2 + end in if (axyz < 1) || (cxyz < 1) then p2 else let v =