10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-22 04:13:33 +01:00

Better vectorization

This commit is contained in:
Anthony Scemama 2018-02-11 01:05:30 +01:00
parent 58d0e65b4e
commit 03380374a1

View File

@ -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 =