From 768e1ac5d896afa2192d10deb20039698f3c3c84 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 20 Nov 2022 22:29:58 +0100 Subject: [PATCH] Vectorizing integrals --- .../prim_int_gauss_gauss.irp.f | 2 +- src/utils/integration.irp.f | 135 +++++++++++++----- src/utils/one_e_integration.irp.f | 2 +- 3 files changed, 104 insertions(+), 35 deletions(-) diff --git a/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f b/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f index ae62a86c..4cba0965 100644 --- a/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f +++ b/src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f @@ -90,7 +90,7 @@ subroutine overlap_gauss_r12_v(D_center_,delta,A_center,B_center,power_A,power_B thr = 1.d-10 d(:) = 0 - maxab = maxval(d(1:3)) + maxab = maxval(power_A(1:3)) double precision, allocatable :: D_center(:,:) allocate(D_center(3,n_points)) diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index 7bf73c15..29c6d1ff 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -161,52 +161,82 @@ subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center,p,fact_k,iorde call gaussian_product_v(alpha,A_center,beta,B_center,fact_k,p,P_center,n_points) +!TODO TRANSP +double precision, allocatable :: P_a_(:,:,:), P_b_(:,:,:), A_center_(:,:), P_center_(:,:), P_new_(:,:,:) +allocate(A_center_(n_points,3), P_center_(n_points,3), P_new_(n_points,0:ldp,3)) +A_center_(1:n_points,1) = A_center(1,1:n_points) +A_center_(1:n_points,2) = A_center(2,1:n_points) +A_center_(1:n_points,3) = A_center(3,1:n_points) +P_center_(1:n_points,1) = P_center(1,1:n_points) +P_center_(1:n_points,2) = P_center(2,1:n_points) +P_center_(1:n_points,3) = P_center(3,1:n_points) + if ( ior(ior(b(1),b(2)),b(3)) == 0 ) then ! b == (0,0,0) + lda = maxval(a) ldb = 0 - allocate(P_a(0:lda,3,n_points),P_b(0:0,3,n_points)) + allocate(P_a_(n_points,0:lda,3), P_b_(n_points,0:0,3)) - call recentered_poly2_v0(P_a,lda,A_center,P_center,a,P_b,B_center,P_center,n_points) + + call recentered_poly2_v0(P_a_,lda,A_center_,P_center_,a,P_b_,B_center,P_center_,n_points) iorder(1:3) = a(1:3) do ipoint=1,n_points do xyz=1,3 - P_new(0,xyz,ipoint) = P_b(0,xyz,ipoint) * P_a(0,xyz,ipoint) + P_new_(ipoint,0,xyz) = P_a_(ipoint,0,xyz) * P_b_(ipoint,0,xyz) do i=1,a(xyz) - P_new(i,xyz,ipoint) = P_new(i,xyz,ipoint) + P_b(0,xyz,ipoint) * P_a(i,xyz,ipoint) + P_new_(ipoint,i,xyz) = P_new_(ipoint,i,xyz) + P_b_(ipoint,0,xyz) * P_a_(ipoint,i,xyz) enddo enddo enddo + +do ipoint=1,n_points +do i=0,ldp +P_new(i,1,ipoint) = P_new_(ipoint,i,1) +P_new(i,2,ipoint) = P_new_(ipoint,i,2) +P_new(i,3,ipoint) = P_new_(ipoint,i,3) +enddo +enddo return endif lda = maxval(a) ldb = maxval(b) - allocate(P_a(0:lda,3,n_points),P_b(0:ldb,3,n_points)) - call recentered_poly2_v(P_a,lda,A_center,P_center,a,P_b,ldb,B_center,P_center,b,n_points) + allocate(P_a_(n_points,0:lda,3), P_b_(n_points,0:ldb,3)) + + call recentered_poly2_v(P_a_,lda,A_center_,P_center_,a,P_b_,ldb,B_center,P_center_,b,n_points) + iorder(1:3) = a(1:3) + b(1:3) + do xyz=1,3 if (b(xyz) == 0) then do ipoint=1,n_points - P_new(0,xyz,ipoint) = P_b(0,xyz,ipoint) * P_a(0,xyz,ipoint) + P_new_(ipoint,0,xyz) = P_a_(ipoint,0,xyz) * P_b_(ipoint,0,xyz) do i=1,a(xyz) - P_new(i,xyz,ipoint) = P_new(i,xyz,ipoint) + P_b(0,xyz,ipoint) * P_a(i,xyz,ipoint) + P_new_(ipoint,i,xyz) = P_new_(ipoint,i,xyz) + P_b_(ipoint,0,xyz) * P_a_(ipoint,i,xyz) enddo enddo else - do ipoint=1,n_points - do i=0,iorder(xyz) - P_new(i,xyz,ipoint) = 0.d0 + do i=0,iorder(xyz) + do ipoint=1,n_points + P_new_(ipoint,i,xyz) = 0.d0 enddo - n_new=0 - call multiply_poly(P_a(0,xyz,ipoint),a(xyz),P_b(0,xyz,ipoint),b(xyz),P_new(0,xyz,ipoint),n_new) enddo + call multiply_poly_v(P_a_(1,0,xyz), a(xyz),P_b_(1,0,xyz),b(xyz),P_new_(1,0,xyz),ldp,n_points) endif enddo + +do ipoint=1,n_points +do i=0,ldp +P_new(i,1,ipoint) = P_new_(ipoint,i,1) +P_new(i,2,ipoint) = P_new_(ipoint,i,2) +P_new(i,3,ipoint) = P_new_(ipoint,i,3) +enddo +enddo end !- @@ -487,6 +517,45 @@ subroutine multiply_poly(b,nb,c,nc,d,nd) end +subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points) + implicit none + BEGIN_DOC + ! Multiply pairs of polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nb, nc, n_points + integer, intent(in) :: nd + double precision, intent(in) :: b(n_points,0:nb), c(n_points,0:nc) + double precision, intent(inout) :: d(n_points,0:nd) + + integer :: ib, ic, id, k, ipoint + if (nd < nb+nc) then + print *, nd, nb, nc + print *, irp_here, ': nd < nb+nc' + stop 1 + endif + + do ic = 0,nc + do ipoint=1, n_points + d(ipoint,ic) = d(ipoint,ic) + c(ipoint,ic) * b(ipoint,0) + enddo + enddo + + do ib=1,nb + do ipoint=1, n_points + d(ipoint, ib) = d(ipoint, ib) + c(ipoint,0) * b(ipoint, ib) + enddo + do ic = 1,nc + do ipoint=1, n_points + d(ipoint, ib+ic) = d(ipoint, ib+ic) + c(ipoint,ic) * b(ipoint, ib) + enddo + enddo + enddo + +end + + subroutine add_poly(b,nb,c,nc,d,nd) implicit none BEGIN_DOC @@ -593,8 +662,8 @@ subroutine recentered_poly2_v(P_new,lda,x_A,x_P,a,P_new2,ldb,x_B,x_Q,b,n_points) ! Recenter two polynomials END_DOC integer, intent(in) :: a(3),b(3), n_points, lda, ldb - double precision, intent(in) :: x_A(3,n_points),x_P(3,n_points),x_B(3),x_Q(3,n_points) - double precision, intent(out) :: P_new(0:lda,3,n_points),P_new2(0:ldb,3,n_points) + double precision, intent(in) :: x_A(n_points,3),x_P(n_points,3),x_B(3),x_Q(n_points,3) + double precision, intent(out) :: P_new(n_points,0:lda,3),P_new2(n_points,0:ldb,3) double precision :: binom_func integer :: i,j,k,l, minab(3), maxab(3),ipoint, xyz double precision, allocatable :: pows_a(:,:), pows_b(:,:) @@ -610,9 +679,9 @@ subroutine recentered_poly2_v(P_new,lda,x_A,x_P,a,P_new2,ldb,x_B,x_Q,b,n_points) if ((a(xyz)<0).or.(b(xyz)<0) ) cycle do ipoint=1,n_points pows_a(ipoint,0) = 1.d0 - pows_a(ipoint,1) = (x_P(xyz,ipoint) - x_A(xyz,ipoint)) + pows_a(ipoint,1) = (x_P(ipoint,xyz) - x_A(ipoint,xyz)) pows_b(ipoint,0) = 1.d0 - pows_b(ipoint,1) = (x_Q(xyz,ipoint) - x_B(xyz)) + pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz)) enddo do i = 2,maxab(xyz) do ipoint=1,n_points @@ -621,39 +690,39 @@ subroutine recentered_poly2_v(P_new,lda,x_A,x_P,a,P_new2,ldb,x_B,x_Q,b,n_points) enddo enddo do ipoint=1,n_points - P_new (0,xyz,ipoint) = pows_a(ipoint,a(xyz)) - P_new2(0,xyz,ipoint) = pows_b(ipoint,b(xyz)) + P_new (ipoint,0,xyz) = pows_a(ipoint,a(xyz)) + P_new2(ipoint,0,xyz) = pows_b(ipoint,b(xyz)) enddo do i = 1,min(minab(xyz),20) fa = binom_transp(a(xyz)-i,a(xyz)) fb = binom_transp(b(xyz)-i,b(xyz)) do ipoint=1,n_points - P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i) - P_new2(i,xyz,ipoint) = fb * pows_b(ipoint,b(xyz)-i) + P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i) + P_new2(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i) enddo enddo do i = minab(xyz)+1,min(a(xyz),20) fa = binom_transp(a(xyz)-i,a(xyz)) do ipoint=1,n_points - P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i) + P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i) enddo enddo do i = minab(xyz)+1,min(b(xyz),20) fb = binom_transp(b(xyz)-i,b(xyz)) do ipoint=1,n_points - P_new2(i,xyz,ipoint) = fb * pows_b(ipoint,b(xyz)-i) + P_new2(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i) enddo enddo do i = 21,a(xyz) fa = binom_func(a(xyz),a(xyz)-i) do ipoint=1,n_points - P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i) + P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i) enddo enddo do i = 21,b(xyz) fb = binom_func(b(xyz),b(xyz)-i) do ipoint=1,n_points - P_new2(i,xyz,ipoint) = fb * pows_b(ipoint,b(xyz)-i) + P_new2(ipoint,i,xyz) = fb * pows_b(ipoint,b(xyz)-i) enddo enddo enddo @@ -666,8 +735,8 @@ subroutine recentered_poly2_v0(P_new,lda,x_A,x_P,a,P_new2,x_B,x_Q,n_points) ! Recenter two polynomials. Special case for b=(0,0,0) END_DOC integer, intent(in) :: a(3), n_points, lda - double precision, intent(in) :: x_A(3,n_points),x_P(3,n_points),x_B(3),x_Q(3,n_points) - double precision, intent(out) :: P_new(0:lda,3,n_points),P_new2(3,n_points) + double precision, intent(in) :: x_A(n_points,3),x_P(n_points,3),x_B(3),x_Q(n_points,3) + double precision, intent(out) :: P_new(n_points,0:lda,3),P_new2(n_points,3) double precision :: binom_func integer :: i,j,k,l, xyz, ipoint, maxab(3) double precision, allocatable :: pows_a(:,:), pows_b(:,:) @@ -681,9 +750,9 @@ subroutine recentered_poly2_v0(P_new,lda,x_A,x_P,a,P_new2,x_B,x_Q,n_points) if (a(xyz)<0) cycle do ipoint=1,n_points pows_a(ipoint,0) = 1.d0 - pows_a(ipoint,1) = (x_P(xyz,ipoint) - x_A(xyz,ipoint)) + pows_a(ipoint,1) = (x_P(ipoint,xyz) - x_A(ipoint,xyz)) pows_b(ipoint,0) = 1.d0 - pows_b(ipoint,1) = (x_Q(xyz,ipoint) - x_B(xyz)) + pows_b(ipoint,1) = (x_Q(ipoint,xyz) - x_B(xyz)) enddo do i = 2,maxab(xyz) do ipoint=1,n_points @@ -692,19 +761,19 @@ subroutine recentered_poly2_v0(P_new,lda,x_A,x_P,a,P_new2,x_B,x_Q,n_points) enddo enddo do ipoint=1,n_points - P_new (0,xyz,ipoint) = pows_a(ipoint,a(xyz)) - P_new2(xyz,ipoint) = pows_b(ipoint,0) + P_new (ipoint,0,xyz) = pows_a(ipoint,a(xyz)) + P_new2(ipoint,xyz) = pows_b(ipoint,0) enddo do i = 1,min(a(xyz),20) fa = binom_transp(a(xyz)-i,a(xyz)) do ipoint=1,n_points - P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i) + P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i) enddo enddo do i = 21,a(xyz) fa = binom_func(a(xyz),a(xyz)-i) do ipoint=1,n_points - P_new (i,xyz,ipoint) = fa * pows_a(ipoint,a(xyz)-i) + P_new (ipoint,i,xyz) = fa * pows_a(ipoint,a(xyz)-i) enddo enddo diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f index ef0c4ffc..2a32d1b9 100644 --- a/src/utils/one_e_integration.irp.f +++ b/src/utils/one_e_integration.irp.f @@ -176,7 +176,7 @@ subroutine overlap_gaussian_xyz_v(A_center,B_center,alpha,beta,power_A,& integer :: nmax double precision :: F_integral - ldp = max_dim + ldp = maxval( power_A(1:3) + power_B(1:3) ) allocate(P_new(0:ldp,3,n_points), P_center(3,n_points), fact_p(n_points), & fact_pp(n_points), pp(n_points))