From e9c25e7fe43cdbd81a1a8a3d2fe0322b34255d1b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 2 Feb 2018 11:59:30 +0100 Subject: [PATCH] Introducing tuples in Vectorized --- Basis/TwoElectronRRVectorized.ml | 333 ++++++++++++++++--------------- 1 file changed, 174 insertions(+), 159 deletions(-) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 307370e..d971aec 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -34,156 +34,155 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0_v m angMom_a = function - | 1 -> let i = if angMom_a.(0) = 1 then 0 else if angMom_a.(1) = 1 then 1 else 2 + | 1 -> + let xyz = + match angMom_a with + | (1,_,_) -> 0 + | (_,1,_) -> 1 + | _ -> 2 in - let f = expo_b *. (Coordinate.coord center_ab i) in + let f = expo_b *. (Coordinate.coord center_ab xyz) in Array.mapi (fun k c -> c *. expo_inv_p *. - ( (Coordinate.coord center_pq.(k) i) *. zero_m_array.(k).(m+1) + ( (Coordinate.coord center_pq.(k) xyz) *. zero_m_array.(k).(m+1) -. f *. zero_m_array.(k).(m) ) ) coef_prod | 0 -> Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod | totAngMom_a -> - let key = Zkey.of_int_tuple (Zkey.Three - (angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) ) + let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in - let (found, result) = - try (true, Zmap.find map_1d.(m) key) with - | Not_found -> (false, - let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] - and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] - and xyz = - match (angMom_a.(0),angMom_a.(1),angMom_a.(2)) with - | (0,0,_) -> 2 - | (0,_,_) -> 1 - | _ -> 0 - in - am.(xyz) <- am.(xyz) - 1; - amm.(xyz) <- amm.(xyz) - 2; - if am.(xyz) < 0 then empty else - let v1 = - let f = - -. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) - in - if (abs_float f < cutoff) then empty else - Array.mapi (fun k v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) ) - in - let p1 = - - Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p *. (Coordinate.coord center_pq.(k) xyz) *. v2k) (vrr0_v (m+1) am (totAngMom_a-1)) - in - if amm.(xyz) < 0 then p1 else - let f = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5 - in - if (abs_float f < cutoff) then empty else - let v1 = vrr0_v m amm (totAngMom_a-2) - in - let v2 = - if (abs_float (f *. expo_inv_p)) < cutoff then empty else - vrr0_v (m+1) amm (totAngMom_a-2) - in - Array.mapi (fun k _ -> p1.(k) +. - f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) coef_prod - ) - in - if not found then - Zmap.add map_1d.(m) key result; - result + try Zmap.find map_1d.(m) key with + | Not_found -> + let result = + let am, amm, amxyz, xyz = + match angMom_a with + | (x,0,0) -> (* 28_336 *) (x-1,0,0),(x-2,0,0), x-1, 0 + | (x,y,0) -> (* 52_221 *) (x,y-1,0),(x,y-2,0), y-1, 1 + | (x,y,z) -> (* 87_215 *) (x,y,z-1),(x,y,z-2), z-1, 2 + in + if amxyz < 0 then empty else + let v1 = + let f = + -. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) + in + if (abs_float f < cutoff) then empty else + Array.mapi (fun k v1k -> f *. v1k) (vrr0_v m am (totAngMom_a-1) ) + in + let p1 = + + Array.mapi (fun k v2k -> v1.(k) +. expo_inv_p *. (Coordinate.coord center_pq.(k) xyz) *. v2k) (vrr0_v (m+1) am (totAngMom_a-1)) + in + if amxyz < 1 then p1 else + let f = (float_of_int amxyz) *. expo_inv_p *. 0.5 + in + if (abs_float f < cutoff) then empty else + let v1 = vrr0_v m amm (totAngMom_a-2) + in + let v2 = + if (abs_float (f *. expo_inv_p)) < cutoff then empty else + vrr0_v (m+1) amm (totAngMom_a-2) + in + Array.mapi (fun k _ -> p1.(k) +. + f *. (v1.(k) +. v2.(k) *. expo_inv_p ) ) coef_prod + in Zmap.add map_1d.(m) key result; + result and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = match (totAngMom_a, totAngMom_c) with - | (i,0) -> if (i=0) then - Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod - else + | (i,0) -> if (i>0) then vrr0_v m angMom_a totAngMom_a + else + (* TODO : TRANSPOSE zero_m_array *) + Array.mapi (fun k c -> c *. zero_m_array.(k).(m)) coef_prod | (_,_) -> - let key = Zkey.of_int_tuple (Zkey.Six - ((angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1), - (angMom_c.(0)+1, angMom_c.(1)+1, angMom_c.(2)+1)) ) - + let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c)) in - let (found, result) = - try (true, Zmap.find map_2d.(m) key) with - | Not_found -> (false, - let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] - and cm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] - and cmm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] - and xyz = - match (angMom_c.(0),angMom_c.(1),angMom_c.(2)) with - | (0,0,_) -> 2 - | (0,_,_) -> 1 - | _ -> 0 - in - am.(xyz) <- am.(xyz) - 1; - cm.(xyz) <- cm.(xyz) - 1; - cmm.(xyz) <- cmm.(xyz) - 2; - if cm.(xyz) < 0 then - empty - else - let f1 = - Array.mapi (fun k _ -> - expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) ) expo_inv_q - in - let f2 = - Array.mapi (fun k _ -> - expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) ) expo_inv_q - in - let v1 = - if (at_least_one_valid f1) then - vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) - else empty - and v2 = - if (at_least_one_valid f2) then - vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) - else empty - in - let p1 = - Array.mapi (fun k _ -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k)) coef_prod - in - let p2 = - if cmm.(xyz) < 0 then p1 else - let fcm = - (float_of_int cm.(xyz)) *. 0.5 - in - let f1 = - Array.mapi (fun k _ -> fcm *. expo_inv_q.(k) ) coef_prod - in - let f2 = - Array.mapi (fun k _ -> f1.(k) *. expo_inv_q.(k) ) coef_prod - in - let v1 = - if (at_least_one_valid f1) then - vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) - else empty - in - let v2 = - if (at_least_one_valid f2) then - vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) - else empty - in - Array.mapi (fun k _ -> p1.(k) +. f1.(k) *. v1.(k) +. f2.(k) *. v2.(k)) coef_prod - in - if (am.(xyz) < 0) || (cm.(xyz) < 0) then p2 else - let fa = - (float_of_int angMom_a.(xyz)) *. expo_inv_p *. 0.5 - in - let f1 = - Array.mapi (fun k _ -> fa *. expo_inv_q.(k) ) coef_prod - in - if (at_least_one_valid f1) then - let v = - vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) - in - Array.mapi (fun k _ -> p2.(k) -. f1.(k) *. v.(k)) coef_prod - else p2 - - ) - in - if not found then - Zmap.add map_2d.(m) key result; + try Zmap.find map_2d.(m) key with + | Not_found -> + let result = + begin + let am, cm, cmm, axyz, cxyz, xyz = + let angMom_ax, angMom_ay, angMom_az = angMom_a + and angMom_cx, angMom_cy, angMom_cz = angMom_c in + match angMom_c with + | (_,0,0) -> (* 321_984 *) + (angMom_ax-1, angMom_ay, angMom_az), + (angMom_cx-1, angMom_cy, angMom_cz), + (angMom_cx-2, angMom_cy, angMom_cz), + angMom_ax,angMom_cx, 0 + | (_,_,0) -> (* 612_002 *) + (angMom_ax, angMom_ay-1, angMom_az), + (angMom_cx, angMom_cy-1, angMom_cz), + (angMom_cx, angMom_cy-2, angMom_cz), + angMom_ay,angMom_cy, 1 + | _ -> (* 1_067_324 *) + (angMom_ax, angMom_ay, angMom_az-1), + (angMom_cx, angMom_cy, angMom_cz-1), + (angMom_cx, angMom_cy, angMom_cz-2), + angMom_az,angMom_cz, 2 + in + if cxyz < 1 then empty else + let f1 = + Array.mapi (fun k _ -> + expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) ) expo_inv_q + in + let f2 = + Array.mapi (fun k _ -> + expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz) ) expo_inv_q + in + let v1 = + if (at_least_one_valid f1) then + vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) + else empty + and v2 = + if (at_least_one_valid f2) then + vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + else empty + in + let p1 = + Array.mapi (fun k _ -> -. v1.(k) *. f1.(k) -. v2.(k) *. f2.(k)) coef_prod + in + let p2 = + if cxyz < 2 then p1 else + let fcm = + (float_of_int (cxyz-1)) *. 0.5 + in + let f1 = + Array.mapi (fun k _ -> fcm *. expo_inv_q.(k) ) coef_prod + in + let f2 = + Array.mapi (fun k _ -> f1.(k) *. expo_inv_q.(k) ) coef_prod + in + let v1 = + if (at_least_one_valid f1) then + vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) + else empty + in + let v2 = + if (at_least_one_valid f2) then + vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) + else empty + in + Array.mapi (fun k _ -> p1.(k) +. f1.(k) *. v1.(k) +. f2.(k) *. v2.(k)) coef_prod + in + if (axyz < 1) || (cxyz < 1) then p2 else + let fa = + (float_of_int axyz) *. expo_inv_p *. 0.5 + in + let f1 = + Array.mapi (fun k _ -> fa *. expo_inv_q.(k) ) coef_prod + in + if (at_least_one_valid f1) then + let v = + vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) + in + Array.mapi (fun k _ -> p2.(k) -. f1.(k) *. v.(k)) coef_prod + else p2 + end + in Zmap.add map_2d.(m) key result; result @@ -201,9 +200,15 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | (_,0) -> vrr0_v 0 angMom_a totAngMom_a | (_,_) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c end - | 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; + | 1 -> + let angMom_ax, angMom_ay, angMom_az = angMom_a in + (* 5_045_008 *) + let ap, xyz = + match angMom_b with + | (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), 0 + | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), 1 + | (_,_,_) -> (angMom_ax,angMom_ay,angMom_az+1), 2 + in let f = Coordinate.coord center_ab xyz in let v1 = vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c @@ -214,17 +219,23 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2 | _ -> - let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] - and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |] - and xyz = - match (angMom_b.(0),angMom_b.(1),angMom_b.(2)) with - | (0,0,_) -> 2 - | (0,_,_) -> 1 - | _ -> 0 + let angMom_ax, angMom_ay, angMom_az = angMom_a + and angMom_bx, angMom_by, angMom_bz = angMom_b in + (* 3_403_478 *) + let bxyz, xyz = + match angMom_b with + | (1,_,_) -> angMom_bx, 0 + | (_,1,_) -> angMom_by, 1 + | (_,_,_) -> angMom_bz, 2 in - ap.(xyz) <- ap.(xyz) + 1; - bm.(xyz) <- bm.(xyz) - 1; - if (bm.(xyz) < 0) then empty else + if (bxyz < 1) then empty else + let ap, bm = + match xyz with + | 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) + | 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) + | _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) + in + let h1 = hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in @@ -238,19 +249,19 @@ let hvrr_two_e_vector (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_v 0 angMom_a angMom_c totAngMom_a totAngMom_c - | (_,0) -> hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c + | (_,0) -> if (totAngMom_b = 0) then + vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c + else + hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> - let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] - and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |] - and xyz = - match (angMom_d.(0),angMom_d.(1),angMom_d.(2)) with - | (0,0,_) -> 2 - | (0,_,_) -> 1 - | _ -> 0 + let (angMom_cx, angMom_cy, angMom_cz) = angMom_c + and (angMom_dx, angMom_dy, angMom_dz) = angMom_d in + let cp, dm, xyz = + match angMom_d with + | (_,0,0) -> (* 1_524_451 *) (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), 0 + | (_,_,0) -> (* 1_302_937 *) (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), 1 + | _ -> (* 1_302_937 *) (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), 2 in - cp.(xyz) <- cp.(xyz) + 1; - dm.(xyz) <- dm.(xyz) - 1; let h1 = hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) and h2 = @@ -258,8 +269,12 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in Array.mapi (fun k center_cd -> h1.(k) +. h2.(k) *. (Coordinate.coord center_cd xyz)) center_cd in - hrr_v angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b - totAngMom_c totAngMom_d + hrr_v + (angMom_a.(0),angMom_a.(1),angMom_a.(2)) + (angMom_b.(0),angMom_b.(1),angMom_b.(2)) + (angMom_c.(0),angMom_c.(1),angMom_c.(2)) + (angMom_d.(0),angMom_d.(1),angMom_d.(2)) + totAngMom_a totAngMom_b totAngMom_c totAngMom_d