From 9e9fe01f8de1c4a11fb7ab7283426b7a9d678e06 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 10 Feb 2018 01:07:34 +0100 Subject: [PATCH] Removed Mat and Vec --- Basis/ERI.ml | 4 +- Basis/TwoElectronRRVectorized.ml | 270 ++++++++++++++++--------------- 2 files changed, 143 insertions(+), 131 deletions(-) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index d01c185..58ee371 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -4,7 +4,7 @@ open Util open Constants open Bigarray -type t = (float, float64_elt, fortran_layout) Bigarray.Genarray.t +type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t (** (00|00)^m : Fundamental electron repulsion integral $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $ @@ -114,7 +114,7 @@ let of_basis basis = (* 4D data initialization *) let eri_array = - Genarray.create Float64 fortran_layout [| n ; n ; n ; n|] + Genarray.create Float32 fortran_layout [| n ; n ; n ; n|] in Genarray.fill eri_array 0.; diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 7cefe8e..01d45c7 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -8,6 +8,11 @@ let cutoff2 = cutoff *. cutoff exception NullQuartet exception Found +let at_least_one_valid arr = + try + Array.iter (fun x -> if (abs_float x > cutoff) then raise Found) arr ; false + with Found -> true + (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) @@ -16,11 +21,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) - coef_prod map_1d map_2d + coef_prod map_1d map_2d np nq = - let nq = Mat.dim1 coef_prod in - let np = Mat.dim2 coef_prod in let empty = Array.make nq 0. @@ -35,23 +38,7 @@ 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 xyz = - match angMom_a with - | (1,_,_) -> 0 - | (_,1,_) -> 1 - | _ -> 2 - in - let a = Mat.mul center_pq.(xyz) zero_m_array.(m+1) in - let b = copy expo_b in - scal (Coordinate.coord center_ab xyz) b; - let c = Mat.map (fun x -> x) zero_m_array.(m) in - Mat.scal_cols c b; - let d = Mat.sub a c in - Mat.scal_cols d expo_inv_p; - Some (Mat.mul coef_prod d) - *) + (* | 1 -> let xyz = match angMom_a with @@ -60,14 +47,18 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | _ -> 2 in Some ( - Array.init np (fun ab -> let l=ab+1 in + Array.init np (fun l -> Array.init nq (fun k -> - let f = expo_b.{l} *. (Coordinate.coord center_ab xyz) in - coef_prod.{k+1,l} *. expo_inv_p.{l} *. - (center_pq.(xyz).{k+1,l} *. zero_m_array.(m+1).{k+1,l} - -. f *. zero_m_array.(m).{k+1,l} ) )) + let f = expo_b.(l) *. (Coordinate.coord center_ab xyz) in + coef_prod.(l).(k) *. expo_inv_p.(l) *. + (center_pq.(xyz).(l).(k) *. zero_m_array.(m+1).(l).(k) + -. f *. zero_m_array.(m).(l).(k) ) )) ) - | 0 -> Some (Mat.mul zero_m_array.(m) coef_prod |> Mat.transpose_copy |> Mat.to_array) + *) + | 0 -> Some ( + Array.init np (fun l -> + Array.init nq (fun k -> + zero_m_array.(m).(l).(k) *. coef_prod.(l).(k) ) ) ) | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in @@ -99,42 +90,42 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in Some ( - Array.init np (fun ab -> let l = ab+1 in + Array.init np (fun l -> let v1 = let f = - -. expo_b.{l} *. expo_inv_p.{l} *. (Coordinate.coord center_ab xyz) + -. expo_b.(l) *. expo_inv_p.(l) *. (Coordinate.coord center_ab xyz) in match v1_top with | Some v1_top -> - v1_top.(l-1) + v1_top.(l) |> Array.map (fun x -> f *. x) | None -> empty in let p1 = match p1_top with | Some p1_top -> - p1_top.(l-1) + p1_top.(l) | _ -> assert false in let p1 = Array.init nq (fun k -> - v1.(k) +. expo_inv_p.{l} *. center_pq.(xyz).{k+1,l} *. p1.(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 + 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-1) + | Some v1_top2 -> v1_top2.(l) | None -> assert false in let v2 = match p1_top2 with - | Some p1_top2 -> p1_top2.(l-1) + | 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} ) ) + p1.(k) +. f *. (v1.(k) +. v2.(k) *. expo_inv_p.(l) ) ) ) ) in Zmap.add map_1d.(m) key result; @@ -145,14 +136,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) match (totAngMom_a, totAngMom_c) with | (i,0) -> - if (i>0) then - begin - match vrr0_v m angMom_a totAngMom_a with - | Some x -> Some (Mat.of_array x |> Mat.transpose_copy) - | None -> None - end - else - Some (Mat.mul zero_m_array.(m) coef_prod ) + vrr0_v m angMom_a totAngMom_a | (_,_) -> let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c)) @@ -183,34 +167,37 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in let f1 = let f = (Coordinate.coord center_cd xyz) in - Vec.init nq (fun k -> - expo_d.{k} *. expo_inv_q.{k} *. f) + Array.init nq (fun k -> expo_d.(k) *. expo_inv_q.(k) *. f) in let v1 = - if (abs_float @@ amax f1 > cutoff) then + if (at_least_one_valid f1) then vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) else None in let f2 = - Mat.init_cols nq np (fun k l -> - expo_inv_q.{k} *. center_pq.(xyz).{k,l} ) + Array.init np (fun l -> + Array.init nq (fun k -> + expo_inv_q.(k) *. center_pq.(xyz).(l).(k) ) ) in let v2 = - if (Mat.as_vec f2 |> amax |> abs_float) < cutoff then - None - else + if (at_least_one_valid (Array.to_list f2 |> Array.concat)) then vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + else + None in let p1 = match v1, v2 with | Some v1, Some v2 -> - Some (Mat.init_cols nq np (fun k l -> - -. v1.{k,l} *. f1.{k} -. v2.{k,l} *. f2.{k,l}) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> + -. v1.(l).(k) *. f1.(k) -. v2.(l).(k) *. f2.(l).(k)) ) ) | None, Some v2 -> - Some (Mat.init_cols nq np (fun k l -> -. v2.{k,l} *. f2.{k,l}) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> -. v2.(l).(k) *. f2.(l).(k)) ) ) | Some v1, None -> - Some (Mat.init_cols nq np (fun k l -> -. v1.{k,l} *. f1.{k} ) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> -. v1.(l).(k) *. f1.(k) ) ) ) | None, None -> None in @@ -220,35 +207,43 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (float_of_int (cxyz-1)) *. 0.5 in let f1 = - Vec.map (fun e -> fcm *. e) expo_inv_q + Array.init nq (fun k -> fcm *. expo_inv_q.(k)) in let f2 = - Vec.mul f1 expo_inv_q + Array.mapi (fun k x -> x *. expo_inv_q.(k)) f1 in let v1 = - if (abs_float @@ amax f1 > cutoff) then + if (at_least_one_valid f1) then vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) else None in let v2 = - if (abs_float @@ amax f2 > cutoff) then + if (at_least_one_valid f2) then vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) else None in match p1, v1, v2 with | Some p1, Some v1, Some v2 -> - Some (Mat.init_cols nq np (fun k l -> p1.{k,l} +. f1.{k} *. v1.{k,l} +. f2.{k} *. v2.{k,l}) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> + p1.(l).(k) +. f1.(k) *. v1.(l).(k) +. f2.(k) *. v2.(l).(k)) ) ) | Some p1, Some v1, None -> - Some (Mat.init_cols nq np (fun k l -> p1.{k,l} +. f1.{k} *. v1.{k,l} )) + Some (Array.init np (fun l -> + Array.init nq (fun k -> p1.(l).(k) +. f1.(k) *. v1.(l).(k) ) ) ) | Some p1, None, Some v2 -> - Some (Mat.init_cols nq np (fun k l -> p1.{k,l} +. f2.{k} *. v2.{k,l}) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> p1.(l).(k) +. f2.(k) *. v2.(l).(k)) ) ) | None , Some v1, Some v2 -> - Some (Mat.init_cols nq np (fun k l -> f1.{k} *. v1.{k,l} +. f2.{k} *. v2.{k,l}) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> + f1.(k) *. v1.(l).(k) +. f2.(k) *. v2.(l).(k)) ) ) | Some p1, None, None -> Some p1 | None , Some v1, None -> - Some (Mat.init_cols nq np (fun k l -> f1.{k} *. v1.{k,l})) + Some (Array.init np (fun l -> + Array.init nq (fun k -> f1.(k) *. v1.(l).(k) ) ) ) | None, None, Some v2 -> - Some (Mat.init_cols nq np (fun k l -> f2.{k} *. v2.{k,l}) ) + Some (Array.init np (fun l -> + Array.init nq (fun k -> f2.(k) *. v2.(l).(k)) ) ) | None, None, None -> None in if (axyz < 1) || (cxyz < 1) then p2 else @@ -258,28 +253,20 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) begin match (p2, v) with | Some p2, Some v -> Some ( - Array.init np (fun ab -> let l = ab+1 in + Array.init np (fun l -> let fa = - (float_of_int axyz) *. expo_inv_p.{l} *. 0.5 + (float_of_int axyz) *. expo_inv_p.(l) *. 0.5 in - let f1 = - Vec.map (fun e -> fa *. e ) expo_inv_q - in - Vec.init nq (fun k -> p2.{k,l} -. f1.{k} *. v.{k,l}) - ) - |> Mat.of_col_vecs ) + Array.init nq (fun k -> p2.(l).(k) -. expo_inv_q.(k) *. fa *. v.(l).(k)) + ) ) | Some p2, None -> Some p2 | None, Some v -> Some ( - Array.init np (fun ab -> let l = ab+1 in + Array.init np (fun l -> let fa = - (float_of_int axyz) *. expo_inv_p.{l} *. 0.5 + (float_of_int axyz) *. expo_inv_p.(l) *. 0.5 in - let f1 = - Vec.map (fun e -> fa *. e ) expo_inv_q - in - Vec.init nq (fun k -> -. f1.{k} *. v.{k,l}) - ) - |> Mat.of_col_vecs ) + Array.init nq (fun k -> -. fa *. expo_inv_q.(k) *. v.(l).(k)) + ) ) | None, None -> None end end @@ -297,16 +284,25 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | 0 -> begin match (totAngMom_a, totAngMom_c) with - | (0,0) -> Mat.gemm_trace zero_m_array.(0) coef_prod + | (0,0) -> + begin + let result = ref 0. in + for l=0 to np-1 do + for k=0 to nq-1 do + result := !result +. zero_m_array.(0).(l).(k) *. coef_prod.(l).(k) + done + done; + !result + end | (_,0) -> begin match vrr0_v 0 angMom_a totAngMom_a with - | Some matrix -> Mat.sum (Mat.of_array matrix) + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end | (_,_) -> match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with - | Some matrix -> Mat.sum matrix + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end | 1 -> @@ -320,13 +316,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let f = Coordinate.coord center_ab xyz in let v1 = match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with - | Some matrix -> Mat.sum matrix + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. in if (abs_float f < cutoff) then v1 else let v2 = match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with - | Some matrix -> Mat.sum matrix + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. in v1 +. v2 *. f @@ -364,7 +360,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | (_,0) -> if (totAngMom_b = 0) then begin match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with - | Some matrix -> Mat.sum matrix + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end else @@ -428,30 +424,30 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (* Compute all integrals in the shell for each pair of significant shell pairs *) - let expo_inv_p = - Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp - |> Vec.of_array - and expo_inv_q = - Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq - |> Vec.of_array - in - let np, nq = - Vec.dim expo_inv_p, Vec.dim expo_inv_q - in - let coef = - let result = Mat.make0 nq np in - Lacaml.D.ger - (Vec.of_array shell_q.ContractedShellPair.coef) - (Vec.of_array shell_p.ContractedShellPair.coef) - result; - result - in - begin match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, totAngMom shell_c, totAngMom shell_d) with | Angular_momentum.(S,S,S,S) -> contracted_class.(0) <- + let expo_inv_p = + Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp + |> Vec.of_array + and expo_inv_q = + Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq + |> Vec.of_array + in + let np, nq = + Vec.dim expo_inv_p, Vec.dim expo_inv_q + in + + let coef = + let result = Mat.make0 nq np in + Lacaml.D.ger + (Vec.of_array shell_q.ContractedShellPair.coef) + (Vec.of_array shell_p.ContractedShellPair.coef) + result; + result + in let zm_array = Mat.init_cols np nq (fun i j -> (** Screening on the product of coefficients *) try @@ -476,47 +472,63 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q ) in Mat.gemm_trace zm_array coef | _ -> + let expo_inv_p = + Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp + and expo_inv_q = + Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq + in + let np, nq = + Array.length expo_inv_p, Array.length expo_inv_q + in let expo_b = Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp - |> Vec.of_array and expo_d = Array.map (fun shell_cd -> Contracted_shell.expo shell_d shell_cd.ShellPair.j) sq - |> Vec.of_array in let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in + let coef = + Array.init np (fun l -> + Array.init nq (fun k -> + shell_q.ContractedShellPair.coef.(k) *. + shell_p.ContractedShellPair.coef.(l) + ) ) + in let center_pq = Array.init 3 (fun xyz -> - Mat.init_cols nq np (fun cd ab -> - let shell_ab = sp.(ab-1) - and shell_cd = sq.(cd-1) - in - let cpq = - Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) - in - match xyz with - | 0 -> Coordinate.x cpq; - | 1 -> Coordinate.y cpq; - | 2 -> Coordinate.z cpq; - | _ -> assert false - ) + Array.init np (fun ab -> + let shell_ab = sp.(ab) in + Array.init nq (fun cd -> + let shell_cd = sq.(cd) + in + let cpq = + Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) + in + match xyz with + | 0 -> Coordinate.x cpq; + | 1 -> Coordinate.y cpq; + | 2 -> Coordinate.z cpq; + | _ -> assert false + ) ) + ) in let zero_m_array = let result = - Array.init (maxm+1) (fun _ -> Mat.make0 nq np) + Array.init (maxm+1) (fun _ -> + Array.init np (fun _ -> Array.make nq 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+1} +. expo_inv_q.{cd+1} + expo_inv_p.(ab) +. expo_inv_q.(cd) in let norm_pq_sq = - center_pq.(0).{cd+1,ab+1} *. center_pq.(0).{cd+1,ab+1} +. - center_pq.(1).{cd+1,ab+1} *. center_pq.(1).{cd+1,ab+1} +. - center_pq.(2).{cd+1,ab+1} *. center_pq.(2).{cd+1,ab+1} + 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 @@ -530,8 +542,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in (* Transpose result *) for m=0 to maxm do - for cd=1 to nq do - result.(m).{cd,ab+1} <- zero_m_array_tmp.(cd-1).(m) + for cd=0 to nq-1 do + result.(m).(ab).(cd) <- zero_m_array_tmp.(cd).(m) done done ) sp; @@ -568,7 +580,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (expo_inv_p, expo_inv_q) (shell_p.ContractedShellPair.center_ab, shell_q.ContractedShellPair.center_ab, center_pq) - coef map_1d map_2d + coef map_1d map_2d np nq in contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) ) class_indices