From c70e4591a9874a41f93795ebe20e0fb040bb1727 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 7 Apr 2014 20:01:30 +0200 Subject: [PATCH] Tests for integration --- scripts/run_tests.py | 8 +- scripts/unarchive_ezfio.sh | 3 + src/Ezfio_files/ezfio.irp.f | 2 +- src/Utils/LinearAlgebra.irp.f | 41 +- src/Utils/README.rst | 390 ++++++++- src/Utils/integration.irp.f | 1026 ++++++++++++------------ src/Utils/one_e_integration.irp.f | 252 +++--- src/Utils/sort.irp.f | 485 ++++++----- src/Utils/tests/Makefile | 33 + src/Utils/tests/test_integration.irp.f | 130 +++ src/Utils/tests/test_integration.ref | 514 ++++++++++++ src/Utils/util.irp.f | 207 ++--- 12 files changed, 2091 insertions(+), 1000 deletions(-) create mode 100644 src/Utils/tests/Makefile create mode 100644 src/Utils/tests/test_integration.irp.f create mode 100644 src/Utils/tests/test_integration.ref diff --git a/scripts/run_tests.py b/scripts/run_tests.py index 8043ce83..f8ec160f 100755 --- a/scripts/run_tests.py +++ b/scripts/run_tests.py @@ -25,6 +25,8 @@ def run_test(test_name,inp): template = """ class $test(unittest.TestCase): + default_precision = 1.e-10 + execfile('$test.ref') def setUp(self): @@ -38,9 +40,13 @@ class $test(unittest.TestCase): continue l,r = buffer l,r = l.strip(), eval(r) + if 'precision' in self.__dict__: + precision = self.precision[l] + else: + precision = self.default_precision if type(r) == float: self.assertAlmostEqual(self.data[inp][l], r, - places=abs(int(log10(self.precision[l]*max(abs(self.data[inp][l]),1.e-12)))), msg=None) + places=abs(int(log10(precision*max(abs(self.data[inp][l]),1.e-12)))), msg=None) else: self.assertEqual(self.data[inp][l], r, msg=None) diff --git a/scripts/unarchive_ezfio.sh b/scripts/unarchive_ezfio.sh index debc8527..7f70fa6e 100755 --- a/scripts/unarchive_ezfio.sh +++ b/scripts/unarchive_ezfio.sh @@ -35,6 +35,7 @@ then exit 1 fi +VERSION=$( cut -d '=' -f 2 < ${QPACKAGE_ROOT}/EZFIO/version) for i in ${!key[@]} do MD5=${key[$i]} @@ -42,6 +43,7 @@ do if [[ ! -d $file ]] then mkdir -p $(dirname $file) + echo ${VERSION} > $(dirname $file)/.version fi if [[ ! -f ${QPACKAGE_ROOT}/data/cache/${MD5} ]] then @@ -49,3 +51,4 @@ do fi cp ${QPACKAGE_ROOT}/data/cache/${MD5} ${file} done +echo ${VERSION} > ${EZFIO_FILE}.ezfio/.version diff --git a/src/Ezfio_files/ezfio.irp.f b/src/Ezfio_files/ezfio.irp.f index 79e71590..9432020f 100644 --- a/src/Ezfio_files/ezfio.irp.f +++ b/src/Ezfio_files/ezfio.irp.f @@ -21,7 +21,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ] ! Check that file exists logical :: exists - inquire(file=trim(ezfio_filename)//'/ezfio/.version',exist=exists) + inquire(file=trim(ezfio_filename)//'/ezfio/creation',exist=exists) if (.not.exists) then print *, 'Error: file '//trim(ezfio_filename)//' does not exist' stop 1 diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index b6023179..5912f0ac 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -1,7 +1,8 @@ subroutine ortho_lowdin(overlap,lda,n,C,ldc,m) implicit none - + BEGIN_DOC ! Compute U.S^-1/2 canonical orthogonalization + END_DOC integer, intent(in) :: lda, ldc, n, m double precision, intent(in) :: overlap(lda,n) @@ -70,8 +71,10 @@ end subroutine get_pseudo_inverse(A,m,n,C,LDA) - ! Find C = A^-1 implicit none + BEGIN_DOC + ! Find C = A^-1 + END_DOC integer, intent(in) :: m,n, LDA double precision, intent(in) :: A(LDA,n) double precision, intent(out) :: C(n,m) @@ -97,7 +100,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA) call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info) if (info /= 0) then print *, info, ': SVD failed' - stop + stop 1 endif do i=1,n @@ -122,8 +125,10 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA) end subroutine find_rotation(A,LDA,B,m,C,n) - ! Find A.C = B implicit none + BEGIN_DOC + ! Find A.C = B + END_DOC integer, intent(in) :: m,n, LDA double precision, intent(in) :: A(LDA,n), B(LDA,n) double precision, intent(out) :: C(n,n) @@ -138,10 +143,11 @@ subroutine find_rotation(A,LDA,B,m,C,n) end - - subroutine apply_rotation(A,LDA,R,LDR,B,LDB,m,n) implicit none + BEGIN_DOC + ! Apply the rotation found by find_rotation + END_DOC double precision, intent(in) :: R(LDR,n) double precision, intent(in) :: A(LDA,n) double precision, intent(out) :: B(LDB,n) @@ -149,8 +155,11 @@ subroutine apply_rotation(A,LDA,R,LDR,B,LDB,m,n) call dgemm('N','N',m,n,n,1.d0,A,LDA,R,LDR,0.d0,B,LDB) end -subroutine jacobi_lapack(eigvalues,eigvectors,H,nmax,n) +subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) implicit none + BEGIN_DOC + ! Diagonalize matrix H + END_DOC integer, intent(in) :: n,nmax double precision, intent(out) :: eigvectors(nmax,n) double precision, intent(out) :: eigvalues(n) @@ -159,31 +168,19 @@ subroutine jacobi_lapack(eigvalues,eigvectors,H,nmax,n) double precision,allocatable :: work(:) double precision,allocatable :: A(:,:) !eigvectors(i,j) = where d_i is the basis function and psi_j is the j th eigenvector - print*,nmax,n allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax)) integer :: LWORK, info, i,j,l,k - double precision :: cpu_2, cpu_1 A=H - call cpu_time (cpu_1) LWORK = 4*nmax - call dsyev( 'V', & - 'U', & - n, & - A, & - nmax, & - eigenvalues, & - work, & - LWORK, & - info ) + call dsyev( 'V', 'U', n, A, nmax, eigenvalues, work, LWORK, info ) if (info < 0) then print *, irp_here, ': the ',-info,'-th argument had an illegal value' - stop + stop 1 else if (info > 0) then print *, irp_here, ': the algorithm failed to converge; ',info,' off-diagonal' print *, 'elements of an intermediate tridiagonal form did not converge to zero.' - stop + stop 1 endif - call cpu_time (cpu_2) eigvectors = 0.d0 eigvalues = 0.d0 do j = 1, n diff --git a/src/Utils/README.rst b/src/Utils/README.rst index 43d12161..c7dec5b5 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -2,19 +2,387 @@ Utils Module ============ +Contains general purpose utilities. -Assumptions ------------ - -.. include:: ./ASSUMPTIONS.rst - - -Needed Modules --------------- - -.. include:: ./NEEDED_MODULES - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f index c3ebaf5b..2f06687d 100644 --- a/src/Utils/integration.irp.f +++ b/src/Utils/integration.irp.f @@ -1,579 +1,547 @@ - subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) -! subroutine that transform the product of -! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) -! into -! fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2) - implicit none - integer, intent(in) :: dim - integer, intent(in) :: a,b ! powers : (x-xa)**a_x = (x-A(1))**a(1) - double precision, intent(in) :: alpha, beta ! exponents - double precision, intent(in) :: A_center ! A center - double precision, intent(in) :: B_center ! B center - double precision, intent(out) :: P_center ! new center - double precision, intent(out) :: p ! new exponent - double precision, intent(out) :: fact_k ! constant factor - include 'constants.F' - double precision, intent(out) :: P_new(0:max_dim) ! polynom - integer, intent(out) :: iorder ! order of the polynoms - double precision :: P_a(0:max_dim), P_b(0:max_dim) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b +subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) + BEGIN_DOC + ! Transform the product of + ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) + ! into + ! fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2) + END_DOC + implicit none + include 'constants.F' + integer, intent(in) :: dim + integer, intent(in) :: a,b ! powers : (x-xa)**a_x = (x-A(1))**a(1) + double precision, intent(in) :: alpha, beta ! exponents + double precision, intent(in) :: A_center ! A center + double precision, intent(in) :: B_center ! B center + double precision, intent(out) :: P_center ! new center + double precision, intent(out) :: p ! new exponent + double precision, intent(out) :: fact_k ! constant factor + double precision, intent(out) :: P_new(0:max_dim) ! polynomial + integer, intent(out) :: iorder ! order of the polynomials - integer :: n_new,i,j - P_new = 0.d0 -! you do the gaussiant product to get the new center and the new exponent - double precision :: p_inv,ab,d_AB - p = alpha+beta - p_inv = 1.d0/p - ab = alpha * beta - d_AB = (A_center - B_center) * (A_center - B_center) - P_center = (alpha * A_center + beta * B_center) * p_inv - fact_k = exp(-ab*p_inv * d_AB) -! you recenter the polynomw P_a and P_b on x -!call recentered_poly(P_a(0),A_center,P_center,a,dim) -!call recentered_poly(P_b(0),B_center,P_center,b,dim) - !DIR$ FORCEINLINE - call recentered_poly2(P_a(0),A_center,P_center,a,P_b(0),B_center,P_center,b) - n_new = 0 - !DEC$ FORCEINLINE - call multiply_poly(P_a(0),a,P_b(0),b,P_new(0),n_new) - iorder = a + b - end + double precision :: P_a(0:max_dim), P_b(0:max_dim) + integer :: n_new,i,j + double precision :: p_inv,ab,d_AB + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b + + ! Do the gaussian product to get the new center and the new exponent + P_new = 0.d0 + p = alpha+beta + p_inv = 1.d0/p + ab = alpha * beta + d_AB = (A_center - B_center) * (A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + fact_k = exp(-ab*p_inv * d_AB) + + ! Recenter the polynomials P_a and P_b on x + !DIR$ FORCEINLINE + call recentered_poly2(P_a(0),A_center,P_center,a,P_b(0),B_center,P_center,b) + n_new = 0 + + !DEC$ FORCEINLINE + call multiply_poly(P_a(0),a,P_b(0),b,P_new(0),n_new) + iorder = a + b +end - subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) -! subroutine that transform the product of -! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) -! into -! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) -! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) -! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) - implicit none - include 'constants.F' - integer, intent(in) :: dim - integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) - double precision, intent(in) :: alpha, beta ! exponents - double precision, intent(in) :: A_center(3) ! A center - double precision, intent(in) :: B_center (3) ! B center - double precision, intent(out) :: P_center(3) ! new center - double precision, intent(out) :: p ! new exponent - double precision, intent(out) :: fact_k ! constant factor - double precision, intent(out) :: P_new(0:max_dim,3) ! polynom - integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynoms - double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b - integer :: n_new,i,j - - iorder(1) = 0 - iorder(2) = 0 - iorder(3) = 0 - P_new(0,1) = 0.d0 - P_new(0,2) = 0.d0 - P_new(0,3) = 0.d0 - - !DEC$ FORCEINLINE - call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) - if (fact_k < 1.d-8) then - fact_k = 0.d0 - return - endif - - !DEC$ FORCEINLINE - call recentered_poly2(P_a(0,1),A_center(1),P_center(1),a(1),P_b(0,1),B_center(1),P_center(1),b(1)) - iorder(1) = a(1) + b(1) - !DIR$ VECTOR ALIGNED - do i=0,iorder(1) - P_new(i,1) = 0.d0 - enddo - n_new=0 - !DEC$ FORCEINLINE - call multiply_poly(P_a(0,1),a(1),P_b(0,1),b(1),P_new(0,1),n_new) - - !DEC$ FORCEINLINE - call recentered_poly2(P_a(0,2),A_center(2),P_center(2),a(2),P_b(0,2),B_center(2),P_center(2),b(2)) - iorder(2) = a(2) + b(2) - !DIR$ VECTOR ALIGNED - do i=0,iorder(2) - P_new(i,2) = 0.d0 - enddo - n_new=0 - !DEC$ FORCEINLINE - call multiply_poly(P_a(0,2),a(2),P_b(0,2),b(2),P_new(0,2),n_new) - - !DEC$ FORCEINLINE - call recentered_poly2(P_a(0,3),A_center(3),P_center(3),a(3),P_b(0,3),B_center(3),P_center(3),b(3)) - iorder(3) = a(3) + b(3) - !DIR$ VECTOR ALIGNED - do i=0,iorder(3) - P_new(i,3) = 0.d0 - enddo - n_new=0 - !DEC$ FORCEINLINE - call multiply_poly(P_a(0,3),a(3),P_b(0,3),b(3),P_new(0,3),n_new) - - end - +subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,beta,a,b,A_center,B_center,dim) + BEGIN_DOC + ! Transforms the product of + ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) + ! into + ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) + ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) + ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) + END_DOC + implicit none + include 'constants.F' + integer, intent(in) :: dim + integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1) + double precision, intent(in) :: alpha, beta ! exponents + double precision, intent(in) :: A_center(3) ! A center + double precision, intent(in) :: B_center (3) ! B center + double precision, intent(out) :: P_center(3) ! new center + double precision, intent(out) :: p ! new exponent + double precision, intent(out) :: fact_k ! constant factor + double precision, intent(out) :: P_new(0:max_dim,3)! polynomial + integer, intent(out) :: iorder(3) ! i_order(i) = order of the polynomials + double precision :: P_a(0:max_dim,3), P_b(0:max_dim,3) + integer :: n_new,i,j + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b + + iorder(1) = 0 + iorder(2) = 0 + iorder(3) = 0 + P_new(0,1) = 0.d0 + P_new(0,2) = 0.d0 + P_new(0,3) = 0.d0 + + !DEC$ FORCEINLINE + call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) + if (fact_k < 1.d-8) then + fact_k = 0.d0 + return + endif + + !DEC$ FORCEINLINE + call recentered_poly2(P_a(0,1),A_center(1),P_center(1),a(1),P_b(0,1),B_center(1),P_center(1),b(1)) + iorder(1) = a(1) + b(1) + !DIR$ VECTOR ALIGNED + do i=0,iorder(1) + P_new(i,1) = 0.d0 + enddo + n_new=0 + !DEC$ FORCEINLINE + call multiply_poly(P_a(0,1),a(1),P_b(0,1),b(1),P_new(0,1),n_new) + + !DEC$ FORCEINLINE + call recentered_poly2(P_a(0,2),A_center(2),P_center(2),a(2),P_b(0,2),B_center(2),P_center(2),b(2)) + iorder(2) = a(2) + b(2) + !DIR$ VECTOR ALIGNED + do i=0,iorder(2) + P_new(i,2) = 0.d0 + enddo + n_new=0 + !DEC$ FORCEINLINE + call multiply_poly(P_a(0,2),a(2),P_b(0,2),b(2),P_new(0,2),n_new) + + !DEC$ FORCEINLINE + call recentered_poly2(P_a(0,3),A_center(3),P_center(3),a(3),P_b(0,3),B_center(3),P_center(3),b(3)) + iorder(3) = a(3) + b(3) + !DIR$ VECTOR ALIGNED + do i=0,iorder(3) + P_new(i,3) = 0.d0 + enddo + n_new=0 + !DEC$ FORCEINLINE + call multiply_poly(P_a(0,3),a(3),P_b(0,3),b(3),P_new(0,3),n_new) + +end subroutine gaussian_product(a,xa,b,xb,k,p,xp) - implicit none -! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} - - double precision , intent(in) :: a,b ! Exponents - double precision , intent(in) :: xa(3),xb(3) ! Centers - double precision , intent(out) :: p ! New exponent - double precision , intent(out) :: xp(3) ! New center - double precision , intent(out) :: k ! Constant - - double precision :: p_inv - - ASSERT (a>0.) - ASSERT (b>0.) - - double precision :: xab(3), ab - !DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab - - p = a+b - p_inv = 1.d0/(a+b) - ab = a*b - xab(1) = xa(1)-xb(1) - xab(2) = xa(2)-xb(2) - xab(3) = xa(3)-xb(3) - ab = ab*p_inv - k = ab*(xab(1)*xab(1)+xab(2)*xab(2)+xab(3)*xab(3)) - if (k > 20.d0) then - k=0.d0 - return - endif - k = exp(-k) - xp(1) = (a*xa(1)+b*xb(1))*p_inv - xp(2) = (a*xa(2)+b*xb(2))*p_inv - xp(3) = (a*xa(3)+b*xb(3))*p_inv + implicit none + BEGIN_DOC + ! Gaussian product in 1D. + ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} + END_DOC + + double precision, intent(in) :: a,b ! Exponents + double precision, intent(in) :: xa(3),xb(3) ! Centers + double precision, intent(out) :: p ! New exponent + double precision, intent(out) :: xp(3) ! New center + double precision, intent(out) :: k ! Constant + + double precision :: p_inv + + ASSERT (a>0.) + ASSERT (b>0.) + + double precision :: xab(3), ab + !DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab + + p = a+b + p_inv = 1.d0/(a+b) + ab = a*b + xab(1) = xa(1)-xb(1) + xab(2) = xa(2)-xb(2) + xab(3) = xa(3)-xb(3) + ab = ab*p_inv + k = ab*(xab(1)*xab(1)+xab(2)*xab(2)+xab(3)*xab(3)) + if (k > 20.d0) then + k=0.d0 + return + endif + k = exp(-k) + xp(1) = (a*xa(1)+b*xb(1))*p_inv + xp(2) = (a*xa(2)+b*xb(2))*p_inv + xp(3) = (a*xa(3)+b*xb(3))*p_inv end subroutine subroutine gaussian_product_x(a,xa,b,xb,k,p,xp) - implicit none -! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} - - double precision , intent(in) :: a,b ! Exponents - double precision , intent(in) :: xa,xb ! Centers - double precision , intent(out) :: p ! New exponent - double precision , intent(out) :: xp ! New center - double precision , intent(out) :: k ! Constant - - double precision :: p_inv - - ASSERT (a>0.) - ASSERT (b>0.) - - double precision :: xab, ab - - p = a+b - p_inv = 1.d0/(a+b) - ab = a*b - xab = xa-xb - ab = ab*p_inv - k = ab*xab*xab - if (k > 20.d0) then - k=0.d0 - return - endif - k = exp(-k) - xp = (a*xa+b*xb)*p_inv + implicit none + BEGIN_DOC + ! Gaussian product in 1D. + ! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} + END_DOC + + double precision , intent(in) :: a,b ! Exponents + double precision , intent(in) :: xa,xb ! Centers + double precision , intent(out) :: p ! New exponent + double precision , intent(out) :: xp ! New center + double precision , intent(out) :: k ! Constant + + double precision :: p_inv + + ASSERT (a>0.) + ASSERT (b>0.) + + double precision :: xab, ab + + p = a+b + p_inv = 1.d0/(a+b) + ab = a*b + xab = xa-xb + ab = ab*p_inv + k = ab*xab*xab + if (k > 20.d0) then + k=0.d0 + return + endif + k = exp(-k) + xp = (a*xa+b*xb)*p_inv end subroutine - - - - + + + + subroutine multiply_poly(b,nb,c,nc,d,nd) - implicit none - ! D(t) += B(t)*C(t) - ! 4251722 + 3928617 - ! 4481076 - ! 185461 - ! 418740 - - integer, intent(in) :: nb, nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:nb), c(0:nc) - double precision, intent(inout) :: d(0:nb+nc) - - integer :: ndtmp - integer :: ib, ic, id, k - if(ior(nc,nb) >= 0) then ! True if nc>=0 and nb>=0 - continue - else - return - endif - ndtmp = nb+nc - - !DIR$ VECTOR ALIGNED - do ic = 0,nc - d(ic) = d(ic) + c(ic) * b(0) - enddo - - do ib=1,nb - d(ib) = d(ib) + c(0) * b(ib) - do ic = 1,nc - d(ib+ic) = d(ib+ic) + c(ic) * b(ib) - enddo - enddo - - do nd = ndtmp,0,-1 - if (d(nd) == 0.d0) then - cycle + implicit none + BEGIN_DOC + ! Multiply two polynomials + ! D(t) += B(t)*C(t) + END_DOC + + integer, intent(in) :: nb, nc + integer, intent(out) :: nd + double precision, intent(in) :: b(0:nb), c(0:nc) + double precision, intent(inout) :: d(0:nb+nc) + + integer :: ndtmp + integer :: ib, ic, id, k + if(ior(nc,nb) >= 0) then ! True if nc>=0 and nb>=0 + continue + else + return endif - exit - enddo - + ndtmp = nb+nc + + !DIR$ VECTOR ALIGNED + do ic = 0,nc + d(ic) = d(ic) + c(ic) * b(0) + enddo + + do ib=1,nb + d(ib) = d(ib) + c(0) * b(ib) + do ic = 1,nc + d(ib+ic) = d(ib+ic) + c(ic) * b(ib) + enddo + enddo + + do nd = ndtmp,0,-1 + if (d(nd) == 0.d0) then + cycle + endif + exit + enddo + end subroutine add_poly(b,nb,c,nc,d,nd) - implicit none - ! D(t) += B(t)+C(t) - integer, intent(inout) :: nb, nc - integer, intent(out) :: nd - double precision, intent(in) :: b(0:nb), c(0:nc) - double precision, intent(out) :: d(0:nb+nc) - - nd = nb+nc - integer :: ib, ic, id - do ib=0,max(nb,nc) - d(ib) = d(ib) + c(ib) + b(ib) - enddo - do while ( (d(nd) == 0.d0).and.(nd>=0) ) + implicit none + BEGIN_DOC + ! Add two polynomials + ! D(t) += B(t)+C(t) + END_DOC + integer, intent(inout) :: nb, nc + integer, intent(out) :: nd + double precision, intent(in) :: b(0:nb), c(0:nc) + double precision, intent(out) :: d(0:nb+nc) + + nd = nb+nc + integer :: ib, ic, id + do ib=0,max(nb,nc) + d(ib) = d(ib) + c(ib) + b(ib) + enddo + do while ( (d(nd) == 0.d0).and.(nd>=0) ) nd -= 1 if (nd < 0) then exit endif - enddo - + enddo + end subroutine add_poly_multiply(b,nb,cst,d,nd) - implicit none - ! D(t) += cst * B(t) - integer, intent(in) :: nb - integer, intent(inout) :: nd - double precision, intent(in) :: b(0:nb),cst - double precision, intent(inout) :: d(0:max(nb,nd)) - - nd = max(nd,nb) - if (nd /= -1) then - integer :: ib, ic, id - do ib=0,nb - d(ib) = d(ib) + cst*b(ib) - enddo - do while ( d(nd) == 0.d0 ) + implicit none + BEGIN_DOC + ! Add a polynomial multiplied by a constant + ! D(t) += cst * B(t) + END_DOC + integer, intent(in) :: nb + integer, intent(inout) :: nd + double precision, intent(in) :: b(0:nb),cst + double precision, intent(inout) :: d(0:max(nb,nd)) + + nd = max(nd,nb) + if (nd /= -1) then + integer :: ib, ic, id + do ib=0,nb + d(ib) = d(ib) + cst*b(ib) + enddo + do while ( d(nd) == 0.d0 ) nd -= 1 if (nd < 0) then exit endif - enddo - endif - + enddo + endif + end subroutine recentered_poly2(P_new,x_A,x_P,a,P_new2,x_B,x_Q,b) -! you enter with (x-x_A)^a -! you leave with sum_i=0,a c_i * (x-x_P)^i ==== P_new(i) * (x-x_P)^i - implicit none - integer, intent(in) :: a,b - double precision, intent(in) :: x_A,x_P,x_B,x_Q - double precision, intent(out) :: P_new(0:a),P_new2(0:b) - double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) - double precision :: binom_func - integer :: i,j,k,l, minab, maxab - if ((a<0).or.(b<0) ) return - maxab = max(a,b) - minab = max(min(a,b),0) - pows_a(0) = 1.d0 - pows_a(1) = (x_P - x_A) - pows_b(0) = 1.d0 - pows_b(1) = (x_Q - x_B) - do i = 2,maxab - pows_a(i) = pows_a(i-1)*pows_a(1) - pows_b(i) = pows_b(i-1)*pows_b(1) - enddo - P_new (0) = pows_a(a) - P_new2(0) = pows_b(b) - do i = 1,min(minab,20) - P_new (i) = binom_transp(a-i,a) * pows_a(a-i) - P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) - enddo - do i = minab+1,min(a,20) - P_new (i) = binom_transp(a-i,a) * pows_a(a-i) - enddo - do i = minab+1,min(b,20) - P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) - enddo - do i = 101,a - P_new(i) = binom_func(a,a-i) * pows_a(a-i) - enddo - do i = 101,b - P_new2(i) = binom_func(b,b-i) * pows_b(b-i) - enddo + implicit none + BEGIN_DOC + ! Recenter two polynomials + END_DOC + integer, intent(in) :: a,b + double precision, intent(in) :: x_A,x_P,x_B,x_Q + double precision, intent(out) :: P_new(0:a),P_new2(0:b) + double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4) + double precision :: binom_func + integer :: i,j,k,l, minab, maxab + if ((a<0).or.(b<0) ) return + maxab = max(a,b) + minab = max(min(a,b),0) + pows_a(0) = 1.d0 + pows_a(1) = (x_P - x_A) + pows_b(0) = 1.d0 + pows_b(1) = (x_Q - x_B) + do i = 2,maxab + pows_a(i) = pows_a(i-1)*pows_a(1) + pows_b(i) = pows_b(i-1)*pows_b(1) + enddo + P_new (0) = pows_a(a) + P_new2(0) = pows_b(b) + do i = 1,min(minab,20) + P_new (i) = binom_transp(a-i,a) * pows_a(a-i) + P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) + enddo + do i = minab+1,min(a,20) + P_new (i) = binom_transp(a-i,a) * pows_a(a-i) + enddo + do i = minab+1,min(b,20) + P_new2(i) = binom_transp(b-i,b) * pows_b(b-i) + enddo + do i = 101,a + P_new(i) = binom_func(a,a-i) * pows_a(a-i) + enddo + do i = 101,b + P_new2(i) = binom_func(b,b-i) * pows_b(b-i) + enddo end - double precision function F_integral(n,p) - ! function that calculates the following integral - ! sum (x) between [-infty;+infty] of x^n exp(-p*x^2) - implicit none - integer :: n - double precision :: p - integer :: i,j - double precision :: accu,sqrt_p,fact_ratio,tmp,fact - include 'include/constants.F' - if(n < 0)then - F_integral = 0.d0 - endif - if(iand(n,1).ne.0)then - F_integral = 0.d0 - return - endif - sqrt_p = 1.d0/dsqrt(p) - if(n==0)then - F_integral = sqpi * sqrt_p - return - endif - F_integral = sqpi * 0.5d0**n * sqrt_p**(n+1) * fact(n)/fact(ishft(n,-1)) - end +double precision function F_integral(n,p) + BEGIN_DOC + ! function that calculates the following integral + ! \int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx + END_DOC + implicit none + integer :: n + double precision :: p + integer :: i,j + double precision :: accu,sqrt_p,fact_ratio,tmp,fact + include 'include/constants.F' + if(n < 0)then + F_integral = 0.d0 + endif + if(iand(n,1).ne.0)then + F_integral = 0.d0 + return + endif + sqrt_p = 1.d0/dsqrt(p) + if(n==0)then + F_integral = sqpi * sqrt_p + return + endif + F_integral = sqpi * 0.5d0**n * sqrt_p**(n+1) * fact(n)/fact(ishft(n,-1)) +end - double precision function rint(n,rho) - implicit none - double precision :: rho,u,rint1,v,val0,rint_large_n,u_inv - integer :: n,k - double precision, parameter :: pi=3.14159265359d0 - double precision, parameter :: dsqpi=1.77245385091d0 - double precision :: two_rho_inv - - -! double precision :: rint_brut -! rint = rint_brut(n,rho,10000) -! return - - if(n.eq.0)then - if(rho == 0.d0)then - rint=1.d0 - else - u_inv=1.d0/dsqrt(rho) - u=rho*u_inv - rint=0.5d0*u_inv*dsqpi*erf(u) - endif - return - endif - if(rho.lt.1.d0)then - rint=rint1(n,rho) - else - if(n.le.20)then - u_inv=1.d0/dsqrt(rho) - v=dexp(-rho) - u=rho*u_inv - two_rho_inv = 0.5d0*u_inv*u_inv - val0=0.5d0*u_inv*dsqpi*erf(u) - rint=(val0-v)*two_rho_inv - do k=2,n - rint=(rint*dfloat(k+k-1)-v)*two_rho_inv - enddo - else - rint=rint_large_n(n,rho) - endif - endif - end - - - - double precision function rint_sum(n_pt_out,rho,d1) - implicit none - integer, intent(in) :: n_pt_out - double precision, intent(in) :: rho,d1(0:n_pt_out) - double precision :: u,rint1,v,val0,rint_large_n,u_inv - integer :: n,k,i - double precision, parameter :: pi=3.14159265359d0 - double precision, parameter :: dsqpi=1.77245385091d0 - double precision :: two_rho_inv, rint_tmp, di - - - if(rho < 1.d0)then - - if(rho == 0.d0)then - rint_sum=d1(0) - else - u_inv=1.d0/dsqrt(rho) - u=rho*u_inv - rint_sum=0.5d0*u_inv*dsqpi*erf(u) *d1(0) - endif - - do i=2,n_pt_out,2 - n = ishft(i,-1) - rint_sum = rint_sum + d1(i)*rint1(n,rho) - enddo - - else - - v=dexp(-rho) - u_inv=1.d0/dsqrt(rho) - u=rho*u_inv - two_rho_inv = 0.5d0*u_inv*u_inv - val0=0.5d0*u_inv*dsqpi*erf(u) - rint_sum=val0*d1(0) - rint_tmp=(val0-v)*two_rho_inv - di = 3.d0 - do i=2,min(n_pt_out,40),2 - rint_sum = rint_sum + d1(i)*rint_tmp - rint_tmp = (rint_tmp*di-v)*two_rho_inv - di = di+2.d0 - enddo - do i=42,n_pt_out,2 - n = ishft(i,-1) - rint_sum = rint_sum + d1(i)*rint_large_n(n,rho) - enddo - - endif - end - - double precision function hermite(n,x) - implicit none - integer :: n,k - double precision :: h0,x,h1,h2 - h0=1.d0 - if(n.eq.0)then - hermite=h0 - return - endif - h1=x+x - if(n.eq.1)then - hermite=h1 - return - endif - do k=1,n-1 - h2=(x+x)*h1-dfloat(k+k)*h0 - h0=h1 - h1=h2 +double precision function rint(n,rho) + implicit none + BEGIN_DOC +!.. math:: +! +! \int_0^1 dx \exp(-p x^2) x^n +! + END_DOC + include 'include/constants.F' + double precision :: rho,u,rint1,v,val0,rint_large_n,u_inv + integer :: n,k + double precision :: two_rho_inv + + if(n.eq.0)then + if(rho == 0.d0)then + rint=1.d0 + else + u_inv=1.d0/dsqrt(rho) + u=rho*u_inv + rint=0.5d0*u_inv*sqpi*erf(u) + endif + return + endif + if(rho.lt.1.d0)then + rint=rint1(n,rho) + else + if(n.le.20)then + u_inv=1.d0/dsqrt(rho) + v=dexp(-rho) + u=rho*u_inv + two_rho_inv = 0.5d0*u_inv*u_inv + val0=0.5d0*u_inv*sqpi*erf(u) + rint=(val0-v)*two_rho_inv + do k=2,n + rint=(rint*dfloat(k+k-1)-v)*two_rho_inv enddo - hermite=h2 - end - - double precision function rint_large_n(n,rho) - implicit none - integer :: n,k,l - double precision :: rho,u,accu,eps,t1,t2,fact,alpha_k,rajout,hermite - u=dsqrt(rho) - accu=0.d0 - k=0 - eps=1.d0 - do while (eps.gt.1.d-15) - t1=1.d0 - do l=0,k - t1=t1*(n+n+l+1.d0) - enddo - t2=0.d0 - do l=0,k - t2=t2+(-1.d0)**l/fact(l+1)/fact(k-l) - enddo - alpha_k=t2*fact(k+1)*fact(k)*(-1.d0)**k - alpha_k= alpha_k/t1 - rajout=(-1.d0)**k*u**k*hermite(k,u)*alpha_k/fact(k) - accu=accu+rajout - eps=dabs(rajout)/accu - k=k+1 - enddo - rint_large_n=dexp(-rho)*accu - end + else + rint=rint_large_n(n,rho) + endif + endif +end - double precision function rint_brut(n,rho,npts) - implicit double precision(a-h,o-z) - double precision :: fi(4), t2(4), accu(4), dt(4), rho4(4), four(4) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fi, t2, accu, dt, rho4, four - accu(1:4)=0.d0 - dt(1)=1.d0/dfloat(npts) - dt(2)=dt(1) - dt(3)=dt(1) - dt(4)=dt(1) - rho4(1:4) = (/ rho, rho, rho, rho /) - fi(1:4)= (/ 0.5d0, 1.5d0, 2.5d0, 3.5d0 /) - four(1:4) = (/ 4.d0, 4.d0, 4.d0, 4.d0 /) - select case(n) - case (1) - do i=1,npts,4 - !DIR$ VECTOR ALIGNED - t2(1:4)=fi(1:4)*dt(1:4) - !DIR$ VECTOR ALIGNED - t2(1:4) = t2(1:4)*t2(1:4) - !DIR$ VECTOR ALIGNED - accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*t2(1:4) - !DIR$ VECTOR ALIGNED - fi(1:4) = fi(1:4)+four(1:4) - enddo - case (2) - do i=1,npts,4 - !DIR$ VECTOR ALIGNED - t2(1:4)=fi(1:4)*dt(1:4) - !DIR$ VECTOR ALIGNED - t2(1:4) = t2(1:4)*t2(1:4) - !DIR$ VECTOR ALIGNED - accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*t2(1:4)*t2(1:4) - !DIR$ VECTOR ALIGNED - fi(1:4) = fi(1:4)+four(1:4) - enddo - case (3) - do i=1,npts,4 - !DIR$ VECTOR ALIGNED - t2(1:4)=fi(1:4)*dt(1:4) - !DIR$ VECTOR ALIGNED - t2(1:4) = t2(1:4)*t2(1:4) - !DIR$ VECTOR ALIGNED - accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*t2(1:4)*t2(1:4)*t2(1:4) - !DIR$ VECTOR ALIGNED - fi(1:4) = fi(1:4)+four(1:4) - enddo - case default - do i=1,npts,4 - !DIR$ VECTOR ALIGNED - t2(1:4)=fi(1:4)*dt(1:4) - !DIR$ VECTOR ALIGNED - t2(1:4) = t2(1:4)*t2(1:4) - !DIR$ VECTOR ALIGNED - accu(1:4)=accu(1:4)+dexp(-rho4(1:4)*t2(1:4))*(t2(1:4)**(n)) - !DIR$ VECTOR ALIGNED - fi(1:4) = fi(1:4)+four(1:4) - enddo - end select - rint_brut=sum(accu)*dt(1) - end - double precision function rint1(n,rho) - implicit none - integer, intent(in) :: n - double precision, intent(in) :: rho - double precision, parameter :: eps=1.d-15 - double precision :: rho_tmp, diff - integer :: k - rint1=inv_int(n+n+1) - rho_tmp = 1.d0 - do k=1,20 - rho_tmp = -rho_tmp*rho - diff=rho_tmp*fact_inv(k)*inv_int(ishft(k+n,1)+1) - rint1=rint1+diff - if (dabs(diff) > eps) then - cycle - endif - return - enddo - write(*,*)'pb in rint1 k too large!' - stop - end +double precision function rint_sum(n_pt_out,rho,d1) + implicit none + BEGIN_DOC + ! Needed for the calculation of two-electron integrals. + END_DOC + include 'include/constants.F' + integer, intent(in) :: n_pt_out + double precision, intent(in) :: rho,d1(0:n_pt_out) + double precision :: u,rint1,v,val0,rint_large_n,u_inv + integer :: n,k,i + double precision :: two_rho_inv, rint_tmp, di + + + if(rho < 1.d0)then + + if(rho == 0.d0)then + rint_sum=d1(0) + else + u_inv=1.d0/dsqrt(rho) + u=rho*u_inv + rint_sum=0.5d0*u_inv*sqpi*erf(u) *d1(0) + endif + + do i=2,n_pt_out,2 + n = ishft(i,-1) + rint_sum = rint_sum + d1(i)*rint1(n,rho) + enddo + + else + + v=dexp(-rho) + u_inv=1.d0/dsqrt(rho) + u=rho*u_inv + two_rho_inv = 0.5d0*u_inv*u_inv + val0=0.5d0*u_inv*sqpi*erf(u) + rint_sum=val0*d1(0) + rint_tmp=(val0-v)*two_rho_inv + di = 3.d0 + do i=2,min(n_pt_out,40),2 + rint_sum = rint_sum + d1(i)*rint_tmp + rint_tmp = (rint_tmp*di-v)*two_rho_inv + di = di+2.d0 + enddo + do i=42,n_pt_out,2 + n = ishft(i,-1) + rint_sum = rint_sum + d1(i)*rint_large_n(n,rho) + enddo + + endif +end + +double precision function hermite(n,x) + implicit none + BEGIN_DOC +! Hermite polynomial + END_DOC + integer :: n,k + double precision :: h0,x,h1,h2 + h0=1.d0 + if(n.eq.0)then + hermite=h0 + return + endif + h1=x+x + if(n.eq.1)then + hermite=h1 + return + endif + do k=1,n-1 + h2=(x+x)*h1-dfloat(k+k)*h0 + h0=h1 + h1=h2 + enddo + hermite=h2 +end + +double precision function rint_large_n(n,rho) + implicit none + BEGIN_DOC +! Version of rint for large values of n + END_DOC + integer :: n,k,l + double precision :: rho,u,accu,eps,t1,t2,fact,alpha_k,rajout,hermite + u=dsqrt(rho) + accu=0.d0 + k=0 + eps=1.d0 + do while (eps.gt.1.d-15) + t1=1.d0 + do l=0,k + t1=t1*(n+n+l+1.d0) + enddo + t2=0.d0 + do l=0,k + t2=t2+(-1.d0)**l/fact(l+1)/fact(k-l) + enddo + alpha_k=t2*fact(k+1)*fact(k)*(-1.d0)**k + alpha_k= alpha_k/t1 + rajout=(-1.d0)**k*u**k*hermite(k,u)*alpha_k/fact(k) + accu=accu+rajout + eps=dabs(rajout)/accu + k=k+1 + enddo + rint_large_n=dexp(-rho)*accu +end + + +double precision function rint1(n,rho) + implicit none + BEGIN_DOC +! Standard version of rint + END_DOC + integer, intent(in) :: n + double precision, intent(in) :: rho + double precision, parameter :: eps=1.d-15 + double precision :: rho_tmp, diff + integer :: k + rint1=inv_int(n+n+1) + rho_tmp = 1.d0 + do k=1,20 + rho_tmp = -rho_tmp*rho + diff=rho_tmp*fact_inv(k)*inv_int(ishft(k+n,1)+1) + rint1=rint1+diff + if (dabs(diff) > eps) then + cycle + endif + return + enddo + write(*,*)'pb in rint1 k too large!' + stop 1 +end diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index ab246ff5..f3a7e0a6 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -1,135 +1,145 @@ - double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_A,power_B,dim) - implicit none -! calculates the following overlap : -! sum (x) between [-infty;+infty] of (x-A_x)^ax (x-B_x)^bx exp(-alpha(x-A_x)^2) exp(-beta(x-B_X)^2) - include 'include/constants.F' - integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynoms - double precision,intent(in) :: A_center,B_center ! center of the x1 functions - integer,intent(in) :: power_A, power_B ! power of the x1 functions - double precision :: P_new(0:max_dim),P_center,fact_p,p,alpha,beta - integer :: iorder_p - call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) - - if(fact_p.lt.0.000001d0)then - overlap_gaussian_x = 0.d0 - return - endif - - overlap_gaussian_x = 0.d0 - integer :: i - double precision :: F_integral - - do i = 0,iorder_p - overlap_gaussian_x += P_new(i) * F_integral(i,p) - enddo - +double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_A,power_B,dim) + implicit none + BEGIN_DOC + !.. math:: + ! + ! \sum_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha(x-A_x)^2) exp(-beta(x-B_X)^2) dx + ! + END_DOC + include 'include/constants.F' + integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynomials + double precision,intent(in) :: A_center,B_center ! center of the x1 functions + integer,intent(in) :: power_A, power_B ! power of the x1 functions + double precision :: P_new(0:max_dim),P_center,fact_p,p,alpha,beta + integer :: iorder_p + call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,& + beta,power_A,power_B,A_center,B_center,dim) + + if(fact_p.lt.0.000001d0)then + overlap_gaussian_x = 0.d0 + return + endif + + overlap_gaussian_x = 0.d0 + integer :: i + double precision :: F_integral + + do i = 0,iorder_p + overlap_gaussian_x += P_new(i) * F_integral(i,p) + enddo + overlap_gaussian_x*= fact_p - end +end - subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim) - implicit none -! .. math:: -! -! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\ -! S = S_x S_y S_z -! -! - include 'include/constants.F' - integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynoms - double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions - double precision, intent(in) :: alpha,beta - integer,intent(in) :: power_A(3), power_B(3) ! power of the x1 functions - double precision, intent(out) :: overlap_x,overlap_y,overlap_z,overlap - double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p - double precision :: F_integral_tab(0:max_dim) - integer :: iorder_p(3) - - call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) - if(fact_p.lt.0.000001d0)then - overlap_x = 0.d0 - overlap_y = 0.d0 - overlap_z = 0.d0 - overlap = 0.d0 - return - endif - integer :: nmax - double precision :: F_integral - nmax = maxval(iorder_p) - do i = 0,nmax - F_integral_tab(i) = F_integral(i,p) - enddo - overlap_x = P_new(0,1) * F_integral_tab(0) - overlap_y = P_new(0,2) * F_integral_tab(0) - overlap_z = P_new(0,3) * F_integral_tab(0) - - integer :: i - do i = 1,iorder_p(1) - overlap_x += P_new(i,1) * F_integral_tab(i) - enddo - call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1)) - overlap_x *= fact_p - - do i = 1,iorder_p(2) - overlap_y += P_new(i,2) * F_integral_tab(i) - enddo - call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2)) - overlap_y *= fact_p - - do i = 1,iorder_p(3) - overlap_z += P_new(i,3) * F_integral_tab(i) - enddo - call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3)) - overlap_z *= fact_p - - overlap = overlap_x * overlap_y * overlap_z - +subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& + power_B,overlap_x,overlap_y,overlap_z,overlap,dim) + implicit none + BEGIN_DOC + !.. math:: + ! + ! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\ + ! S = S_x S_y S_z + ! + END_DOC + include 'include/constants.F' + integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynomials + double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions + double precision, intent(in) :: alpha,beta + integer,intent(in) :: power_A(3), power_B(3) ! power of the x1 functions + double precision, intent(out) :: overlap_x,overlap_y,overlap_z,overlap + double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p + double precision :: F_integral_tab(0:max_dim) + integer :: iorder_p(3) + + call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) + if(fact_p.lt.0.000001d0)then + overlap_x = 0.d0 + overlap_y = 0.d0 + overlap_z = 0.d0 + overlap = 0.d0 + return + endif + integer :: nmax + double precision :: F_integral + nmax = maxval(iorder_p) + do i = 0,nmax + F_integral_tab(i) = F_integral(i,p) + enddo + overlap_x = P_new(0,1) * F_integral_tab(0) + overlap_y = P_new(0,2) * F_integral_tab(0) + overlap_z = P_new(0,3) * F_integral_tab(0) + + integer :: i + do i = 1,iorder_p(1) + overlap_x += P_new(i,1) * F_integral_tab(i) + enddo + call gaussian_product_x(alpha,A_center(1),beta,B_center(1),fact_p,p,P_center(1)) + overlap_x *= fact_p + + do i = 1,iorder_p(2) + overlap_y += P_new(i,2) * F_integral_tab(i) + enddo + call gaussian_product_x(alpha,A_center(2),beta,B_center(2),fact_p,p,P_center(2)) + overlap_y *= fact_p + + do i = 1,iorder_p(3) + overlap_z += P_new(i,3) * F_integral_tab(i) + enddo + call gaussian_product_x(alpha,A_center(3),beta,B_center(3),fact_p,p,P_center(3)) + overlap_z *= fact_p + + overlap = overlap_x * overlap_y * overlap_z + end - subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) - implicit none -! compute the following integral : -! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] - integer :: i,j,k,l - integer,intent(in) :: power_A,power_B - double precision, intent(in) :: lower_exp_val - double precision,intent(in) :: A_center, B_center,alpha,beta - double precision, intent(out) :: overlap_x,dx - integer, intent(in) :: nx - double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho - double precision :: P_center - if(power_A.lt.0.or.power_B.lt.0)then - overlap_x = 0.d0 - dx = 0.d0 - return - endif - p = alpha + beta - p_inv= 1.d0/p - rho = alpha * beta * p_inv - dist = (A_center - B_center)*(A_center - B_center) - P_center = (alpha * A_center + beta * B_center) * p_inv - factor = dexp(-rho * dist) - - double precision :: tmp - - tmp = dsqrt(lower_exp_val/p) +subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + implicit none + BEGIN_DOC + ! .. math :: + ! + ! \int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx + ! + END_DOC + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho + double precision :: P_center + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + + double precision :: tmp + + tmp = dsqrt(lower_exp_val/p) x_min = P_center - tmp x_max = P_center + tmp - domain = x_max-x_min - dx = domain/dble(nx) - overlap_x = 0.d0 - x = x_min - do i = 1, nx - x += dx - overlap_x += abs((x-A_center)**power_A * (x-B_center)**power_B) * dexp(-p * (x-P_center)*(x-P_center)) - enddo - - overlap_x *= factor * dx + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += abs((x-A_center)**power_A * (x-B_center)**power_B) * dexp(-p * (x-P_center)*(x-P_center)) + enddo + + overlap_x = factor * dx * overlap_x end - diff --git a/src/Utils/sort.irp.f b/src/Utils/sort.irp.f index e821895c..171ad2a6 100644 --- a/src/Utils/sort.irp.f +++ b/src/Utils/sort.irp.f @@ -1,42 +1,51 @@ BEGIN_TEMPLATE subroutine insertion_$Xsort (x,iorder,isize) implicit none - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer,intent(in) :: isize - $type :: xtmp - integer :: i, i0, j, jmax - + BEGIN_DOC + ! Sort array x(isize) using the insertion sort algorithm. + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer,intent(in) :: isize + $type :: xtmp + integer :: i, i0, j, jmax + do i=1,isize - xtmp = x(i) - i0 = iorder(i) - j = i-1 - do j=i-1,1,-1 - if ( x(j) > xtmp ) then - x(j+1) = x(j) - iorder(j+1) = iorder(j) - else - exit - endif - enddo - x(j+1) = xtmp - iorder(j+1) = i0 + xtmp = x(i) + i0 = iorder(i) + j = i-1 + do j=i-1,1,-1 + if ( x(j) > xtmp ) then + x(j+1) = x(j) + iorder(j+1) = iorder(j) + else + exit + endif + enddo + x(j+1) = xtmp + iorder(j+1) = i0 enddo - end subroutine insertion_$Xsort subroutine heap_$Xsort(x,iorder,isize) implicit none - $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer,intent(in) :: isize + BEGIN_DOC + ! Sort array x(isize) using the heap sort algorithm. + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer,intent(in) :: isize - integer :: i, k, j, l, i0 - $type :: xtemp - - l = isize/2+1 - k = isize - do while (.True.) + integer :: i, k, j, l, i0 + $type :: xtemp + + l = isize/2+1 + k = isize + do while (.True.) if (l>1) then l=l-1 xtemp = x(l) @@ -48,7 +57,7 @@ BEGIN_TEMPLATE iorder(k) = iorder(1) k = k-1 if (k == 1) then - x(1) = xtemp + x(1) = xtemp iorder(1) = i0 exit endif @@ -56,45 +65,52 @@ BEGIN_TEMPLATE i=l j = ishft(l,1) do while (j1) then l=l-1 xtemp = x(l) @@ -106,7 +122,7 @@ BEGIN_TEMPLATE iorder(k) = iorder(1) k = k-1 if (k == 1) then - x(1) = xtemp + x(1) = xtemp iorder(1) = i0 exit endif @@ -114,39 +130,44 @@ BEGIN_TEMPLATE i=l j = ishft(l,1) do while (j xtmp ) then - x(j+1_8) = x(j) - iorder(j+1_8) = iorder(j) - else - exit - endif - enddo - x(j+1_8) = xtmp - iorder(j+1_8) = i0 + xtmp = x(i) + i0 = iorder(i) + j = i-1_8 + do j=i-1_8,1_8,-1_8 + if ( x(j) > xtmp ) then + x(j+1_8) = x(j) + iorder(j+1_8) = iorder(j) + else + exit + endif + enddo + x(j+1_8) = xtmp + iorder(j+1_8) = i0 enddo - + end subroutine insertion_$Xsort subroutine $Xset_order_big(x,iorder,isize) implicit none - integer*8 :: isize - $type :: x(*) - $type, allocatable :: xtmp(:) - integer*8 :: iorder(*) - integer*8 :: i - allocate(xtmp(isize)) + BEGIN_DOC + ! array A has already been sorted, and iorder has contains the new order of + ! elements of A. This subroutine changes the order of x to match the new order of A. + ! This is a version for very large arrays where the indices need + ! to be in integer*8 format + END_DOC + integer*8 :: isize + $type :: x(*) + $type, allocatable :: xtmp(:) + integer*8 :: iorder(*) + integer*8 :: i + allocate(xtmp(isize)) do i=1_8,isize - xtmp(i) = x(iorder(i)) + xtmp(i) = x(iorder(i)) enddo - + do i=1_8,isize - x(i) = xtmp(i) + x(i) = xtmp(i) enddo deallocate(xtmp) end @@ -246,48 +284,54 @@ END_TEMPLATE BEGIN_TEMPLATE - recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) - implicit none - $int_type, intent(in) :: isize - $int_type, intent(inout) :: iorder(isize) - $type, intent(inout) :: x(isize) - integer, intent(in) :: iradix - integer :: iradix_new - $type, allocatable :: x2(:), x1(:) - $int_type, allocatable :: iorder1(:),iorder2(:) - $int_type :: i0, i1, i2, i3, i - integer, parameter :: integer_size=$octets - $type, parameter :: zero=$zero - $type :: mask - integer :: nthreads, omp_get_num_threads - !DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1 - - if (iradix == -1) then + recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) + implicit none + BEGIN_DOC + ! Sort integer array x(isize) using the radix sort algorithm. + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + ! iradix should be -1 in input. + END_DOC + $int_type, intent(in) :: isize + $int_type, intent(inout) :: iorder(isize) + $type, intent(inout) :: x(isize) + integer, intent(in) :: iradix + integer :: iradix_new + $type, allocatable :: x2(:), x1(:) + $int_type, allocatable :: iorder1(:),iorder2(:) + $int_type :: i0, i1, i2, i3, i + integer, parameter :: integer_size=$octets + $type, parameter :: zero=$zero + $type :: mask + integer :: nthreads, omp_get_num_threads + !DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1 + if (iradix == -1) then + ! Find most significant bit - + i0 = 0_8 i3 = -1_8 - + do i=1,isize i3 = max(i3,x(i)) enddo - + iradix_new = integer_size-1-leadz(i3) mask = ibset(zero,iradix_new) nthreads = 1 -! nthreads = 1+ishft(omp_get_num_threads(),-1) - - integer :: err + ! nthreads = 1+ishft(omp_get_num_threads(),-1) + + integer :: err allocate(x1(isize/nthreads+1),iorder1(isize/nthreads+1),x2(isize/nthreads+1),iorder2(isize/nthreads+1),stat=err) if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays' - stop + print *, irp_here, ': Unable to allocate arrays' + stop endif - + i1=1_8 i2=1_8 - + do i=1,isize if (iand(mask,x(i)) == zero) then iorder1(i1) = iorder(i) @@ -301,7 +345,7 @@ BEGIN_TEMPLATE enddo i1=i1-1_8 i2=i2-1_8 - + do i=1,i1 iorder(i0+i) = iorder1(i) x(i0+i) = x1(i) @@ -310,11 +354,11 @@ BEGIN_TEMPLATE i3 = i0 deallocate(x1,iorder1,stat=err) if (err /= 0) then - print *, irp_here, ': Unable to deallocate arrays x1, iorder1' - stop + print *, irp_here, ': Unable to deallocate arrays x1, iorder1' + stop endif - - + + do i=1,i2 iorder(i0+i) = iorder2(i) x(i0+i) = x2(i) @@ -322,80 +366,79 @@ BEGIN_TEMPLATE i0 = i0+i2 deallocate(x2,iorder2,stat=err) if (err /= 0) then - print *, irp_here, ': Unable to deallocate arrays x2, iorder2' - stop + print *, irp_here, ': Unable to deallocate arrays x2, iorder2' + stop endif - - + + if (i3>1) then call $Xradix_sort$big(x,iorder,i3,iradix_new-1) endif - + if (isize-i3>1) then call $Xradix_sort$big(x(i3+1),iorder(i3+1),isize-i3,iradix_new-1) endif - + return - endif + endif - ASSERT (iradix > 0) - if (isize < 48) then - call insertion_$Xsort$big(x,iorder,isize) - return - endif - - - allocate(x2(isize),iorder2(isize),stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays x1, iorder1' - stop - endif - - - mask = ibset(zero,iradix) - i0=1 - i1=1 - - do i=1,isize - if (iand(mask,x(i)) == zero) then - iorder(i0) = iorder(i) - x(i0) = x(i) - i0 = i0+1 - else - iorder2(i1) = iorder(i) - x2(i1) = x(i) - i1 = i1+1 - endif - enddo - i0=i0-1 - i1=i1-1 - - do i=1,i1 - iorder(i0+i) = iorder2(i) - x(i0+i) = x2(i) - enddo - - deallocate(x2,iorder2,stat=err) - if (err /= 0) then - print *, irp_here, ': Unable to allocate arrays x2, iorder2' - stop - endif - - - if (iradix == 0) then - return - endif - - - if (i1>1) then - call $Xradix_sort$big(x(i0+1),iorder(i0+1),i1,iradix-1) - endif - if (i0>1) then - call $Xradix_sort$big(x,iorder,i0,iradix-1) - endif + ASSERT (iradix > 0) + if (isize < 48) then + call insertion_$Xsort$big(x,iorder,isize) + return + endif - end + allocate(x2(isize),iorder2(isize),stat=err) + if (err /= 0) then + print *, irp_here, ': Unable to allocate arrays x1, iorder1' + stop + endif + + + mask = ibset(zero,iradix) + i0=1 + i1=1 + + do i=1,isize + if (iand(mask,x(i)) == zero) then + iorder(i0) = iorder(i) + x(i0) = x(i) + i0 = i0+1 + else + iorder2(i1) = iorder(i) + x2(i1) = x(i) + i1 = i1+1 + endif + enddo + i0=i0-1 + i1=i1-1 + + do i=1,i1 + iorder(i0+i) = iorder2(i) + x(i0+i) = x2(i) + enddo + + deallocate(x2,iorder2,stat=err) + if (err /= 0) then + print *, irp_here, ': Unable to allocate arrays x2, iorder2' + stop + endif + + + if (iradix == 0) then + return + endif + + + if (i1>1) then + call $Xradix_sort$big(x(i0+1),iorder(i0+1),i1,iradix-1) + endif + if (i0>1) then + call $Xradix_sort$big(x,iorder,i0,iradix-1) + endif + + end SUBST [ X, type, octets, is_big, big, int_type, zero ] i ; integer ; 32 ; .False. ; ; integer ; 0;; diff --git a/src/Utils/tests/Makefile b/src/Utils/tests/Makefile new file mode 100644 index 00000000..77bd84ba --- /dev/null +++ b/src/Utils/tests/Makefile @@ -0,0 +1,33 @@ +OPENMP =1 +PROFILE =0 +DEBUG = 0 + +IRPF90+= -I tests + +REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f)) + +.PHONY: clean executables serial_tests parallel_tests + +all: clean executables serial_tests parallel_tests + +parallel_tests: $(REF_FILES) + @echo ; echo " ---- Running parallel tests ----" ; echo + @OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py + +serial_tests: $(REF_FILES) + @echo ; echo " ---- Running serial tests ----" ; echo + @OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py + +executables: $(wildcard *.irp.f) veryclean + $(MAKE) -C .. + +%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables + $(QPACKAGE_ROOT)/scripts/create_test_ref.sh $* + +clean: + $(MAKE) -C .. clean + +veryclean: + $(MAKE) -C .. veryclean + + diff --git a/src/Utils/tests/test_integration.irp.f b/src/Utils/tests/test_integration.irp.f new file mode 100644 index 00000000..1b2c879e --- /dev/null +++ b/src/Utils/tests/test_integration.irp.f @@ -0,0 +1,130 @@ +program test_integration + implicit none + character*(128) :: arg + integer :: iargc + integer :: i + if (iargc() < 1) then + print *, iargc() + print *, 'usage : test_integration ' + stop 1 + endif + call getarg(1,arg) + i = len(arg) + do while (arg(i:i) == ' ') + i -= 1 + if (i == 0) then + stop 1 + endif + enddo + i -= 1 + do while (arg(i:i) /= '/') + i -= 1 + if (i == 0) then + stop 1 + endif + enddo + i += 1 + arg = arg(i:) + BEGIN_SHELL [ /bin/bash ] + for i in $(grep subroutine tests/test_integration.irp.f | cut -d ' ' -f 2 | sed 's/test_//' ) + do + echo "if (trim(arg) == '"$i"') then" + echo ' call test_'$i + echo 'endif' + done + END_SHELL +end + +subroutine test_rint1 + implicit none + integer :: i,j + double precision :: rho(10) + double precision :: rint1 + do i=1,size(rho) + rho(i) = 2.d0**(1-i) + enddo + do j=1,size(rho) + do i=0,8 + print '(I2,A,F12.8,A3,E15.8)', i, ',', rho(j), ' : ', rint1(i,rho(j)) + enddo + enddo +end + +subroutine test_rint_large_n + implicit none + integer :: i,j + double precision :: rho(10) + double precision :: rint_large_n + do i=1,size(rho) + rho(i) = 2.d0**(2-i) + enddo + do j=1,size(rho) + do i=4,20 + print '(I2,A,F12.8,A3,E15.8)', i, ',', rho(j), ' : ', rint_large_n(i,rho(j)) + enddo + enddo +end + +subroutine test_hermite + implicit none + integer :: i,j + double precision :: x(10) + double precision :: hermite + do i=1,size(x) + x(i) = (-1.d0)**i * 2.d0**(5-i) + enddo + do j=1,size(x) + do i=0,10 + print '(I2,A,F12.8,A3,E15.8)', i, ',', x(j), ' : ', hermite(i,x(j)) + enddo + enddo +end + +subroutine test_rint_sum + implicit none + integer :: i,j + double precision :: d1(0:50), rho(10) + double precision :: rint_sum + do i=0,size(d1)-1 + d1(i) = (-1.d0)**i * 2.d0**(5-i) + enddo + do i=1,size(rho) + rho(i) = abs(1.d0/d1(i)) + enddo + do j=1,size(rho) + do i=0,5 + print '(I2,A,F12.8,A3,E15.8)', 4*i+1, ',', rho(j), ' : ', rint_sum(4*i+1,rho(j),d1) + enddo + enddo +end + +subroutine test_rint + implicit none + integer :: i,j + double precision :: rho(10) + double precision :: rint + do i=1,size(rho) + rho(i) = 2.d0**(2-i) + enddo + do j=1,size(rho) + do i=0,20,3 + print '(I2,A,F12.8,A3,E15.8)', i, ',', rho(j), ' : ', rint(i,rho(j)) + enddo + enddo +end + +subroutine test_F_integral + implicit none + integer :: i,j + double precision :: rho(10) + double precision :: F_integral + do i=1,size(rho) + rho(i) = 2.d0**(2-i) + enddo + do j=1,size(rho) + do i=0,20,3 + print '(I2,A,F12.8,A3,E15.8)', i, ',', rho(j), ' : ', F_integral(i,rho(j)) + enddo + enddo + +end diff --git a/src/Utils/tests/test_integration.ref b/src/Utils/tests/test_integration.ref new file mode 100644 index 00000000..9c8643c0 --- /dev/null +++ b/src/Utils/tests/test_integration.ref @@ -0,0 +1,514 @@ +data = { + 'rint1' : { +'0, 1.00000000' : 0.74682413E+00, +'1, 1.00000000' : 0.18947235E+00, +'2, 1.00000000' : 0.10026880E+00, +'3, 1.00000000' : 0.66732275E-01, +'4, 1.00000000' : 0.49623241E-01, +'5, 1.00000000' : 0.39364865E-01, +'6, 1.00000000' : 0.32567034E-01, +'7, 1.00000000' : 0.27746002E-01, +'8, 1.00000000' : 0.24155294E-01, +'0, 0.50000000' : 0.85562439E+00, +'1, 0.50000000' : 0.24909373E+00, +'2, 0.50000000' : 0.14075054E+00, +'3, 0.50000000' : 0.97222024E-01, +'4, 0.50000000' : 0.74023511E-01, +'5, 0.50000000' : 0.59680941E-01, +'6, 0.50000000' : 0.49959693E-01, +'7, 0.50000000' : 0.42945347E-01, +'8, 0.50000000' : 0.37649547E-01, +'0, 0.25000000' : 0.92256201E+00, +'1, 0.25000000' : 0.28752246E+00, +'2, 0.25000000' : 0.16753319E+00, +'3, 0.25000000' : 0.11773034E+00, +'4, 0.25000000' : 0.90623235E-01, +'5, 0.25000000' : 0.73616662E-01, +'6, 0.25000000' : 0.61964994E-01, +'7, 0.25000000' : 0.53488276E-01, +'8, 0.25000000' : 0.47046726E-01, +'0, 0.12500000' : 0.95985044E+00, +'1, 0.12500000' : 0.30941414E+00, +'2, 0.12500000' : 0.18298209E+00, +'3, 0.12500000' : 0.12965410E+00, +'4, 0.12500000' : 0.10032732E+00, +'5, 0.12500000' : 0.81795915E-01, +'6, 0.12500000' : 0.69032643E-01, +'7, 0.12500000' : 0.59709841E-01, +'8, 0.12500000' : 0.52602850E-01, +'0, 0.06250000' : 0.97955155E+00, +'1, 0.06250000' : 0.32110789E+00, +'2, 0.06250000' : 0.19128479E+00, +'3, 0.06250000' : 0.13608717E+00, +'4, 0.06250000' : 0.10557686E+00, +'5, 0.06250000' : 0.86229246E-01, +'6, 0.06250000' : 0.72869188E-01, +'7, 0.06250000' : 0.63091082E-01, +'8, 0.06250000' : 0.55625318E-01, +'0, 0.03125000' : 0.98968027E+00, +'1, 0.03125000' : 0.32715253E+00, +'2, 0.03125000' : 0.19558951E+00, +'3, 0.03125000' : 0.13942892E+00, +'4, 0.03125000' : 0.10830743E+00, +'5, 0.03125000' : 0.88537500E-01, +'6, 0.03125000' : 0.74868200E-01, +'7, 0.03125000' : 0.64853890E-01, +'8, 0.03125000' : 0.57201824E-01, +'0, 0.01562500' : 0.99481599E+00, +'1, 0.01562500' : 0.33022570E+00, +'2, 0.01562500' : 0.19778136E+00, +'3, 0.01562500' : 0.14113208E+00, +'4, 0.01562500' : 0.10970000E+00, +'5, 0.01562500' : 0.89715269E-01, +'6, 0.01562500' : 0.75888558E-01, +'7, 0.01562500' : 0.65753944E-01, +'8, 0.01562500' : 0.58006946E-01, +'0, 0.00781250' : 0.99740193E+00, +'1, 0.00781250' : 0.33177518E+00, +'2, 0.00781250' : 0.19888731E+00, +'3, 0.00781250' : 0.14199186E+00, +'4, 0.00781250' : 0.11040323E+00, +'5, 0.00781250' : 0.90310159E-01, +'6, 0.00781250' : 0.76404035E-01, +'7, 0.00781250' : 0.66208710E-01, +'8, 0.00781250' : 0.58413795E-01, +'0, 0.00390625' : 0.99869944E+00, +'1, 0.00390625' : 0.33255317E+00, +'2, 0.00390625' : 0.19944281E+00, +'3, 0.00390625' : 0.14242381E+00, +'4, 0.00390625' : 0.11075658E+00, +'5, 0.00390625' : 0.90609118E-01, +'6, 0.00390625' : 0.76663109E-01, +'7, 0.00390625' : 0.66437288E-01, +'8, 0.00390625' : 0.58618300E-01, +'0, 0.00195312' : 0.99934934E+00, +'1, 0.00195312' : 0.33294298E+00, +'2, 0.00195312' : 0.19972119E+00, +'3, 0.00195312' : 0.14264030E+00, +'4, 0.00195312' : 0.11093370E+00, +'5, 0.00195312' : 0.90758978E-01, +'6, 0.00195312' : 0.76792981E-01, +'7, 0.00195312' : 0.66551877E-01, +'8, 0.00195312' : 0.58720824E-01, + }, + 'rint_large_n' : { + '4, 2.00000000' : 0.22769400E-01, + '5, 2.00000000' : 0.17397330E-01, + '6, 2.00000000' : 0.14008836E-01, + '7, 2.00000000' : 0.11694896E-01, + '8, 2.00000000' : 0.10022038E-01, + '9, 2.00000000' : 0.87598405E-02, +'10, 2.00000000' : 0.77754214E-02, +'11, 2.00000000' : 0.69871414E-02, +'12, 2.00000000' : 0.63422421E-02, +'13, 2.00000000' : 0.58051924E-02, +'14, 2.00000000' : 0.53512279E-02, +'15, 2.00000000' : 0.49625818E-02, +'16, 2.00000000' : 0.46261880E-02, +'17, 2.00000000' : 0.43322304E-02, +'18, 2.00000000' : 0.40731954E-02, +'19, 2.00000000' : 0.38432370E-02, +'20, 2.00000000' : 0.36377399E-02, + '4, 1.00000000' : 0.49623241E-01, + '5, 1.00000000' : 0.39364865E-01, + '6, 1.00000000' : 0.32567034E-01, + '7, 1.00000000' : 0.27746002E-01, + '8, 1.00000000' : 0.24155294E-01, + '9, 1.00000000' : 0.21380280E-01, +'10, 1.00000000' : 0.19172936E-01, +'11, 1.00000000' : 0.17376108E-01, +'12, 1.00000000' : 0.15885526E-01, +'13, 1.00000000' : 0.14629351E-01, +'14, 1.00000000' : 0.13556514E-01, +'15, 1.00000000' : 0.12629735E-01, +'16, 1.00000000' : 0.11821172E-01, +'17, 1.00000000' : 0.11109621E-01, +'18, 1.00000000' : 0.10478652E-01, +'19, 1.00000000' : 0.99153382E-02, +'20, 1.00000000' : 0.94093737E-02, + '4, 0.50000000' : 0.74131525E-01, + '5, 0.50000000' : 0.59734080E-01, + '6, 0.50000000' : 0.49988791E-01, + '7, 0.50000000' : 0.42962588E-01, + '8, 0.50000000' : 0.37660400E-01, + '9, 0.50000000' : 0.33518800E-01, +'10, 0.50000000' : 0.30195249E-01, +'11, 0.50000000' : 0.27469686E-01, +'12, 0.50000000' : 0.25194350E-01, +'13, 0.50000000' : 0.23266388E-01, +'14, 0.50000000' : 0.21612012E-01, +'15, 0.50000000' : 0.20176927E-01, +'16, 0.50000000' : 0.18920297E-01, +'17, 0.50000000' : 0.17810821E-01, +'18, 0.50000000' : 0.16824108E-01, +'19, 0.50000000' : 0.15940870E-01, +'20, 0.50000000' : 0.15145655E-01, + '4, 0.25000000' : 0.90623235E-01, + '5, 0.25000000' : 0.73616662E-01, + '6, 0.25000000' : 0.61964994E-01, + '7, 0.25000000' : 0.53488276E-01, + '8, 0.25000000' : 0.47046726E-01, + '9, 0.25000000' : 0.41987104E-01, +'10, 0.25000000' : 0.37908392E-01, +'11, 0.25000000' : 0.34550883E-01, +'12, 0.25000000' : 0.31739030E-01, +'13, 0.25000000' : 0.29349937E-01, +'14, 0.25000000' : 0.27295006E-01, +'15, 0.25000000' : 0.25508764E-01, +'16, 0.25000000' : 0.23941782E-01, +'17, 0.25000000' : 0.22556049E-01, +'18, 0.25000000' : 0.21321854E-01, +'19, 0.25000000' : 0.20215642E-01, +'20, 0.25000000' : 0.19218495E-01, + '4, 0.12500000' : 0.10032732E+00, + '5, 0.12500000' : 0.81795915E-01, + '6, 0.12500000' : 0.69032643E-01, + '7, 0.12500000' : 0.59709841E-01, + '8, 0.12500000' : 0.52602850E-01, + '9, 0.12500000' : 0.47006219E-01, +'10, 0.12500000' : 0.42485051E-01, +'11, 0.12500000' : 0.38756708E-01, +'12, 0.12500000' : 0.35629567E-01, +'13, 0.12500000' : 0.32969128E-01, +'14, 0.12500000' : 0.30678211E-01, +'15, 0.12500000' : 0.28684857E-01, +'16, 0.12500000' : 0.26934646E-01, +'17, 0.12500000' : 0.25385662E-01, +'18, 0.12500000' : 0.24005098E-01, +'19, 0.12500000' : 0.22766909E-01, +'20, 0.12500000' : 0.21650155E-01, + '4, 0.06250000' : 0.10557686E+00, + '5, 0.06250000' : 0.86229246E-01, + '6, 0.06250000' : 0.72869188E-01, + '7, 0.06250000' : 0.63091082E-01, + '8, 0.06250000' : 0.55625318E-01, + '9, 0.06250000' : 0.49738703E-01, +'10, 0.06250000' : 0.44978296E-01, +'11, 0.06250000' : 0.41049216E-01, +'12, 0.06250000' : 0.37751241E-01, +'13, 0.06250000' : 0.34943654E-01, +'14, 0.06250000' : 0.32524670E-01, +'15, 0.06250000' : 0.30418845E-01, +'16, 0.06250000' : 0.28569075E-01, +'17, 0.06250000' : 0.26931342E-01, +'18, 0.06250000' : 0.25471168E-01, +'19, 0.06250000' : 0.24161166E-01, +'20, 0.06250000' : 0.22979305E-01, + '4, 0.03125000' : 0.10830743E+00, + '5, 0.03125000' : 0.88537500E-01, + '6, 0.03125000' : 0.74868200E-01, + '7, 0.03125000' : 0.64853890E-01, + '8, 0.03125000' : 0.57201824E-01, + '9, 0.03125000' : 0.51164511E-01, +'10, 0.03125000' : 0.46279696E-01, +'11, 0.03125000' : 0.42246171E-01, +'12, 0.03125000' : 0.38859267E-01, +'13, 0.03125000' : 0.35975049E-01, +'14, 0.03125000' : 0.33489346E-01, +'15, 0.03125000' : 0.31324909E-01, +'16, 0.03125000' : 0.29423240E-01, +'17, 0.03125000' : 0.27739231E-01, +'18, 0.03125000' : 0.26237537E-01, +'19, 0.03125000' : 0.24890074E-01, +'20, 0.03125000' : 0.23674243E-01, + '4, 0.01562500' : 0.10970000E+00, + '5, 0.01562500' : 0.89715269E-01, + '6, 0.01562500' : 0.75888558E-01, + '7, 0.01562500' : 0.65753944E-01, + '8, 0.01562500' : 0.58006946E-01, + '9, 0.01562500' : 0.51892813E-01, +'10, 0.01562500' : 0.46944559E-01, +'11, 0.01562500' : 0.42857760E-01, +'12, 0.01562500' : 0.39425485E-01, +'13, 0.01562500' : 0.36502162E-01, +'14, 0.01562500' : 0.33982407E-01, +'15, 0.01562500' : 0.31788050E-01, +'16, 0.01562500' : 0.29859885E-01, +'17, 0.01562500' : 0.28152246E-01, +'18, 0.01562500' : 0.26629349E-01, +'19, 0.01562500' : 0.25262753E-01, +'20, 0.01562500' : 0.24029571E-01, + '4, 0.00781250' : 0.11040323E+00, + '5, 0.00781250' : 0.90310159E-01, + '6, 0.00781250' : 0.76404035E-01, + '7, 0.00781250' : 0.66208710E-01, + '8, 0.00781250' : 0.58413795E-01, + '9, 0.00781250' : 0.52260879E-01, +'10, 0.00781250' : 0.47280591E-01, +'11, 0.00781250' : 0.43166888E-01, +'12, 0.00781250' : 0.39711698E-01, +'13, 0.00781250' : 0.36768623E-01, +'14, 0.00781250' : 0.34231665E-01, +'15, 0.00781250' : 0.32022192E-01, +'16, 0.00781250' : 0.30080639E-01, +'17, 0.00781250' : 0.28361060E-01, +'18, 0.00781250' : 0.26827449E-01, +'19, 0.00781250' : 0.25451185E-01, +'20, 0.00781250' : 0.24209234E-01, + '4, 0.00390625' : 0.11075658E+00, + '5, 0.00390625' : 0.90609118E-01, + '6, 0.00390625' : 0.76663109E-01, + '7, 0.00390625' : 0.66437288E-01, + '8, 0.00390625' : 0.58618300E-01, + '9, 0.00390625' : 0.52445898E-01, +'10, 0.00390625' : 0.47449515E-01, +'11, 0.00390625' : 0.43322293E-01, +'12, 0.00390625' : 0.39855587E-01, +'13, 0.00390625' : 0.36902585E-01, +'14, 0.00390625' : 0.34356981E-01, +'15, 0.00390625' : 0.32139911E-01, +'16, 0.00390625' : 0.30191629E-01, +'17, 0.00390625' : 0.28466050E-01, +'18, 0.00390625' : 0.26927053E-01, +'19, 0.00390625' : 0.25545928E-01, +'20, 0.00390625' : 0.24299570E-01, +}, + 'hermite': { + '0,-16.00000000' : 0.10000000E+01, + '1,-16.00000000' : -0.32000000E+02, + '2,-16.00000000' : 0.10220000E+04, + '3,-16.00000000' : -0.32576000E+05, + '4,-16.00000000' : 0.10363000E+07, + '5,-16.00000000' : -0.32900992E+08, + '6,-16.00000000' : 0.10424687E+10, + '7,-16.00000000' : -0.32964188E+11, + '8,-16.00000000' : 0.10402595E+13, + '9,-16.00000000' : -0.32760875E+14, +'10,-16.00000000' : 0.10296233E+16, + '0, 8.00000000' : 0.10000000E+01, + '1, 8.00000000' : 0.16000000E+02, + '2, 8.00000000' : 0.25400000E+03, + '3, 8.00000000' : 0.40000000E+04, + '4, 8.00000000' : 0.62476000E+05, + '5, 8.00000000' : 0.96761600E+06, + '6, 8.00000000' : 0.14857096E+08, + '7, 8.00000000' : 0.22610214E+09, + '8, 8.00000000' : 0.34096350E+10, + '9, 8.00000000' : 0.50936525E+11, +'10, 8.00000000' : 0.75361097E+12, + '0, -4.00000000' : 0.10000000E+01, + '1, -4.00000000' : -0.80000000E+01, + '2, -4.00000000' : 0.62000000E+02, + '3, -4.00000000' : -0.46400000E+03, + '4, -4.00000000' : 0.33400000E+04, + '5, -4.00000000' : -0.23008000E+05, + '6, -4.00000000' : 0.15066400E+06, + '7, -4.00000000' : -0.92921600E+06, + '8, -4.00000000' : 0.53244320E+07, + '9, -4.00000000' : -0.27728000E+08, +'10, -4.00000000' : 0.12598422E+09, + '0, 2.00000000' : 0.10000000E+01, + '1, 2.00000000' : 0.40000000E+01, + '2, 2.00000000' : 0.14000000E+02, + '3, 2.00000000' : 0.40000000E+02, + '4, 2.00000000' : 0.76000000E+02, + '5, 2.00000000' : -0.16000000E+02, + '6, 2.00000000' : -0.82400000E+03, + '7, 2.00000000' : -0.31040000E+04, + '8, 2.00000000' : -0.88000000E+03, + '9, 2.00000000' : 0.46144000E+05, +'10, 2.00000000' : 0.20041600E+06, + '0, -1.00000000' : 0.10000000E+01, + '1, -1.00000000' : -0.20000000E+01, + '2, -1.00000000' : 0.20000000E+01, + '3, -1.00000000' : 0.40000000E+01, + '4, -1.00000000' : -0.20000000E+02, + '5, -1.00000000' : 0.80000000E+01, + '6, -1.00000000' : 0.18400000E+03, + '7, -1.00000000' : -0.46400000E+03, + '8, -1.00000000' : -0.16480000E+04, + '9, -1.00000000' : 0.10720000E+05, +'10, -1.00000000' : 0.82240000E+04, + '0, 0.50000000' : 0.10000000E+01, + '1, 0.50000000' : 0.10000000E+01, + '2, 0.50000000' : -0.10000000E+01, + '3, 0.50000000' : -0.50000000E+01, + '4, 0.50000000' : 0.10000000E+01, + '5, 0.50000000' : 0.41000000E+02, + '6, 0.50000000' : 0.31000000E+02, + '7, 0.50000000' : -0.46100000E+03, + '8, 0.50000000' : -0.89500000E+03, + '9, 0.50000000' : 0.64810000E+04, +'10, 0.50000000' : 0.22591000E+05, + '0, -0.25000000' : 0.10000000E+01, + '1, -0.25000000' : -0.50000000E+00, + '2, -0.25000000' : -0.17500000E+01, + '3, -0.25000000' : 0.28750000E+01, + '4, -0.25000000' : 0.90625000E+01, + '5, -0.25000000' : -0.27531250E+02, + '6, -0.25000000' : -0.76859375E+02, + '7, -0.25000000' : 0.36880469E+03, + '8, -0.25000000' : 0.89162891E+03, + '9, -0.25000000' : -0.63466895E+04, +'10, -0.25000000' : -0.12875976E+05, + '0, 0.12500000' : 0.10000000E+01, + '1, 0.12500000' : 0.25000000E+00, + '2, 0.12500000' : -0.19375000E+01, + '3, 0.12500000' : -0.14843750E+01, + '4, 0.12500000' : 0.11253906E+02, + '5, 0.12500000' : 0.14688477E+02, + '6, 0.12500000' : -0.10886694E+03, + '7, 0.12500000' : -0.20347845E+03, + '8, 0.12500000' : 0.14732676E+04, + '9, 0.12500000' : 0.36239722E+04, +'10, 0.12500000' : -0.25612824E+05, + '0, -0.06250000' : 0.10000000E+01, + '1, -0.06250000' : -0.12500000E+00, + '2, -0.06250000' : -0.19843750E+01, + '3, -0.06250000' : 0.74804688E+00, + '4, -0.06250000' : 0.11812744E+02, + '5, -0.06250000' : -0.74609680E+01, + '6, -0.06250000' : -0.11719482E+03, + '7, -0.06250000' : 0.10418097E+03, + '8, -0.06250000' : 0.16277049E+04, + '9, -0.06250000' : -0.18703586E+04, +'10, -0.06250000' : -0.29064893E+05, + '0, 0.03125000' : 0.10000000E+01, + '1, 0.03125000' : 0.62500000E-01, + '2, 0.03125000' : -0.19960938E+01, + '3, 0.03125000' : -0.37475586E+00, + '4, 0.03125000' : 0.11953140E+02, + '5, 0.03125000' : 0.37451181E+01, + '6, 0.03125000' : -0.11929733E+03, + '7, 0.03125000' : -0.52397501E+02, + '8, 0.03125000' : 0.16668878E+04, + '9, 0.03125000' : 0.94254050E+03, +'10, 0.03125000' : -0.29945072E+05, +}, + 'rint_sum' : { + '1, 0.06250000' : 0.31345650E+02, + '5, 0.06250000' : 0.34297082E+02, + '9, 0.06250000' : 0.34378323E+02, +'13, 0.06250000' : 0.34381587E+02, +'17, 0.06250000' : 0.34381737E+02, +'21, 0.06250000' : 0.34381745E+02, + '1, 0.12500000' : 0.30715214E+02, + '5, 0.12500000' : 0.33556491E+02, + '9, 0.12500000' : 0.33633859E+02, +'13, 0.12500000' : 0.33636955E+02, +'17, 0.12500000' : 0.33637097E+02, +'21, 0.12500000' : 0.33637104E+02, + '1, 0.25000000' : 0.29521984E+02, + '5, 0.25000000' : 0.32157230E+02, + '9, 0.25000000' : 0.32227424E+02, +'13, 0.25000000' : 0.32230208E+02, +'17, 0.25000000' : 0.32230336E+02, +'21, 0.25000000' : 0.32230342E+02, + '1, 0.50000000' : 0.27379981E+02, + '5, 0.50000000' : 0.29654231E+02, + '9, 0.50000000' : 0.29712095E+02, +'13, 0.50000000' : 0.29714351E+02, +'17, 0.50000000' : 0.29714453E+02, +'21, 0.50000000' : 0.29714458E+02, + '1, 1.00000000' : 0.23898372E+02, + '5, 1.00000000' : 0.25614689E+02, + '9, 1.00000000' : 0.25654258E+02, +'13, 1.00000000' : 0.25655742E+02, +'17, 1.00000000' : 0.25655808E+02, +'21, 1.00000000' : 0.25655811E+02, + '1, 2.00000000' : 0.19140608E+02, + '5, 2.00000000' : 0.20172111E+02, + '9, 2.00000000' : 0.20191130E+02, +'13, 2.00000000' : 0.20191783E+02, +'17, 2.00000000' : 0.20191811E+02, +'21, 2.00000000' : 0.20191812E+02, + '1, 4.00000000' : 0.14113302E+02, + '5, 4.00000000' : 0.14571079E+02, + '9, 4.00000000' : 0.14576072E+02, +'13, 4.00000000' : 0.14576208E+02, +'17, 4.00000000' : 0.14576213E+02, +'21, 4.00000000' : 0.14576214E+02, + '1, 8.00000000' : 0.10025878E+02, + '5, 8.00000000' : 0.10189658E+02, + '9, 8.00000000' : 0.10190276E+02, +'13, 8.00000000' : 0.10190285E+02, +'17, 8.00000000' : 0.10190285E+02, +'21, 8.00000000' : 0.10190285E+02, + '1, 16.00000000' : 0.70898153E+01, + '5, 16.00000000' : 0.71465026E+01, + '9, 16.00000000' : 0.71465561E+01, +'13, 16.00000000' : 0.71465563E+01, +'17, 16.00000000' : 0.71465563E+01, +'21, 16.00000000' : 0.71465563E+01, + '1, 32.00000000' : 0.50132565E+01, + '5, 32.00000000' : 0.50330691E+01, + '9, 32.00000000' : 0.50330737E+01, +'13, 32.00000000' : 0.50330737E+01, +'17, 32.00000000' : 0.50330737E+01, +'21, 32.00000000' : 0.50330737E+01, +}, + 'rint' : { + '0, 2.00000000' : 0.59814401E+00, + '3, 2.00000000' : 0.32344698E-01, + '6, 2.00000000' : 0.14008836E-01, + '9, 2.00000000' : 0.87598405E-02, +'12, 2.00000000' : 0.63422421E-02, +'15, 2.00000000' : 0.49625795E-02, +'18, 2.00000000' : 0.40719014E-02, + '0, 1.00000000' : 0.74682413E+00, + '3, 1.00000000' : 0.66732275E-01, + '6, 1.00000000' : 0.32567034E-01, + '9, 1.00000000' : 0.21380280E-01, +'12, 1.00000000' : 0.15885521E-01, +'15, 1.00000000' : 0.12617897E-01, +'18, 1.00000000' : -0.42503468E-01, + '0, 0.50000000' : 0.85562439E+00, + '3, 0.50000000' : 0.97222024E-01, + '6, 0.50000000' : 0.49959693E-01, + '9, 0.50000000' : 0.33511631E-01, +'12, 0.50000000' : 0.25191806E-01, +'15, 0.50000000' : 0.20175809E-01, +'18, 0.50000000' : 0.16823542E-01, + '0, 0.25000000' : 0.92256201E+00, + '3, 0.25000000' : 0.11773034E+00, + '6, 0.25000000' : 0.61964994E-01, + '9, 0.25000000' : 0.41987104E-01, +'12, 0.25000000' : 0.31739030E-01, +'15, 0.25000000' : 0.25508764E-01, +'18, 0.25000000' : 0.21321854E-01, + '0, 0.12500000' : 0.95985044E+00, + '3, 0.12500000' : 0.12965410E+00, + '6, 0.12500000' : 0.69032643E-01, + '9, 0.12500000' : 0.47006219E-01, +'12, 0.12500000' : 0.35629567E-01, +'15, 0.12500000' : 0.28684857E-01, +'18, 0.12500000' : 0.24005098E-01, + '0, 0.06250000' : 0.97955155E+00, + '3, 0.06250000' : 0.13608717E+00, + '6, 0.06250000' : 0.72869188E-01, + '9, 0.06250000' : 0.49738703E-01, +'12, 0.06250000' : 0.37751241E-01, +'15, 0.06250000' : 0.30418845E-01, +'18, 0.06250000' : 0.25471168E-01, + '0, 0.03125000' : 0.98968027E+00, + '3, 0.03125000' : 0.13942892E+00, + '6, 0.03125000' : 0.74868200E-01, + '9, 0.03125000' : 0.51164511E-01, +'12, 0.03125000' : 0.38859267E-01, +'15, 0.03125000' : 0.31324909E-01, +'18, 0.03125000' : 0.26237537E-01, + '0, 0.01562500' : 0.99481599E+00, + '3, 0.01562500' : 0.14113208E+00, + '6, 0.01562500' : 0.75888558E-01, + '9, 0.01562500' : 0.51892813E-01, +'12, 0.01562500' : 0.39425485E-01, +'15, 0.01562500' : 0.31788050E-01, +'18, 0.01562500' : 0.26629349E-01, + '0, 0.00781250' : 0.99740193E+00, + '3, 0.00781250' : 0.14199186E+00, + '6, 0.00781250' : 0.76404035E-01, + '9, 0.00781250' : 0.52260879E-01, +'12, 0.00781250' : 0.39711698E-01, +'15, 0.00781250' : 0.32022192E-01, +'18, 0.00781250' : 0.26827449E-01, + '0, 0.00390625' : 0.99869944E+00, + '3, 0.00390625' : 0.14242381E+00, + '6, 0.00390625' : 0.76663109E-01, + '9, 0.00390625' : 0.52445898E-01, +'12, 0.00390625' : 0.39855587E-01, +'15, 0.00390625' : 0.32139911E-01, +'18, 0.00390625' : 0.26927053E-01, +}, +'F_integral' : { +} +} diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 90d4603f..00ea2cd4 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -1,66 +1,76 @@ double precision function binom_func(i,j) - implicit none - integer,intent(in) :: i,j - double precision :: fact, f - integer, save :: ifirst - double precision, save :: memo(0:15,0:15) - !DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: memo - integer :: k,l - if (ifirst == 0) then - ifirst = 1 - do k=0,15 - f = fact(k) - do l=0,15 - memo(k,l) = f/(fact(l)*fact(k-l)) + implicit none + BEGIN_DOC + !.. math :: + ! + ! \frac{i!}{j!(i-j)!} + ! + END_DOC + integer,intent(in) :: i,j + double precision :: fact, f + integer, save :: ifirst + double precision, save :: memo(0:15,0:15) + !DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: memo + integer :: k,l + if (ifirst == 0) then + ifirst = 1 + do k=0,15 + f = fact(k) + do l=0,15 + memo(k,l) = f/(fact(l)*fact(k-l)) + enddo enddo - enddo - endif - if ( (i<=15).and.(j<=15) ) then - binom_func = memo(i,j) - else - binom_func = fact(i)/(fact(j)*fact(i-j)) - endif + endif + if ( (i<=15).and.(j<=15) ) then + binom_func = memo(i,j) + else + binom_func = fact(i)/(fact(j)*fact(i-j)) + endif end - - BEGIN_PROVIDER [ double precision, binom, (0:20,0:20) ] &BEGIN_PROVIDER [ double precision, binom_transp, (0:20,0:20) ] - implicit none - BEGIN_DOC -! Binomial coefficients - END_DOC - integer :: k,l - double precision :: fact, f - do k=0,20 - f = fact(k) - do l=0,20 - binom(k,l) = f/(fact(l)*fact(k-l)) - binom_transp(l,k) = binom(k,l) + implicit none + BEGIN_DOC + ! Binomial coefficients + END_DOC + integer :: k,l + double precision :: fact, f + do k=0,20 + f = fact(k) + do l=0,20 + binom(k,l) = f/(fact(l)*fact(k-l)) + binom_transp(l,k) = binom(k,l) + enddo enddo - enddo END_PROVIDER integer function align_double(n) - implicit none - integer :: n - include 'include/constants.F' - if (mod(n,SIMD_vector/4) /= 0) then - align_double= n + SIMD_vector/4 - mod(n,SIMD_vector/4) - else - align_double= n - endif + implicit none + BEGIN_DOC + ! Compute 1st dimension such that it is aligned for vectorization. + END_DOC + integer :: n + include 'include/constants.F' + if (mod(n,SIMD_vector/4) /= 0) then + align_double= n + SIMD_vector/4 - mod(n,SIMD_vector/4) + else + align_double= n + endif end double precision function fact(n) implicit none - integer :: n - double precision, save :: memo(1:100) - integer, save :: memomax = 1 - + BEGIN_DOC + ! n! + END_DOC + integer :: n + double precision, save :: memo(1:100) + integer, save :: memomax = 1 + if (n<=memomax) then if (n<2) then fact = 1.d0 @@ -69,8 +79,8 @@ double precision function fact(n) endif return endif - - integer :: i + + integer :: i memo(1) = 1.d0 do i=memomax+1,min(n,100) memo(i) = memo(i-1)*float(i) @@ -85,12 +95,12 @@ end function BEGIN_PROVIDER [ double precision, fact_inv, (128) ] - implicit none - BEGIN_DOC -! 1.d0/fact(k) - END_DOC - integer :: i - double precision :: fact + implicit none + BEGIN_DOC + ! 1/n! + END_DOC + integer :: i + double precision :: fact do i=1,size(fact_inv) fact_inv(i) = 1.d0/fact(i) enddo @@ -98,10 +108,13 @@ END_PROVIDER double precision function dble_fact(n) result(fact2) implicit none - integer :: n - double precision, save :: memo(1:100) - integer, save :: memomax = 1 - + BEGIN_DOC + ! n!! + END_DOC + integer :: n + double precision, save :: memo(1:100) + integer, save :: memomax = 1 + ASSERT (iand(n,1) /= 0) if (n<=memomax) then if (n<3) then @@ -111,49 +124,55 @@ double precision function dble_fact(n) result(fact2) endif return endif - - integer :: i + + integer :: i memo(1) = 1.d0 do i=memomax+2,min(n,99),2 memo(i) = memo(i-2)* float(i) enddo memomax = min(n,99) fact2 = memo(memomax) - + do i=101,n,2 fact2 = fact2*float(i) enddo - + end function subroutine write_git_log(iunit) - implicit none - integer, intent(in) :: iunit - write(iunit,*) '----------------' - write(iunit,*) 'Last git commit:' - BEGIN_SHELL [ /bin/bash ] - git log -1 | sed "s/'//g"| sed "s/^/ write(iunit,*) '/g" | sed "s/$/'/g" - END_SHELL - write(iunit,*) '----------------' + implicit none + BEGIN_DOC + ! Write the last git commit in file iunit. + END_DOC + integer, intent(in) :: iunit + write(iunit,*) '----------------' + write(iunit,*) 'Last git commit:' + BEGIN_SHELL [ /bin/bash ] + git log -1 2>/dev/null | sed "s/'//g"| sed "s/^/ write(iunit,*) '/g" | sed "s/$/'/g" || echo "Unknown" + END_SHELL + write(iunit,*) '----------------' end BEGIN_PROVIDER [ double precision, inv_int, (128) ] - implicit none - BEGIN_DOC -! 1/i - END_DOC - integer :: i - do i=1,size(inv_int) - inv_int(i) = 1.d0/dble(i) - enddo - + implicit none + BEGIN_DOC + ! 1/i + END_DOC + integer :: i + do i=1,size(inv_int) + inv_int(i) = 1.d0/dble(i) + enddo + END_PROVIDER subroutine wall_time(t) implicit none - double precision, intent(out) :: t - integer :: c - integer, save :: rate = 0 + BEGIN_DOC + ! The equivalent of cpu_time, but for the wall time. + END_DOC + double precision, intent(out) :: t + integer :: c + integer, save :: rate = 0 if (rate == 0) then CALL SYSTEM_CLOCK(count_rate=rate) endif @@ -162,17 +181,17 @@ subroutine wall_time(t) end BEGIN_PROVIDER [ integer, nproc ] - implicit none - BEGIN_DOC -! Number of current openmp threads - END_DOC - - integer :: omp_get_num_threads - nproc = 1 - !$OMP PARALLEL - !$OMP MASTER - !$ nproc = omp_get_num_threads() - !$OMP END MASTER - !$OMP END PARALLEL + implicit none + BEGIN_DOC + ! Number of current OpenMP threads + END_DOC + + integer :: omp_get_num_threads + nproc = 1 + !$OMP PARALLEL + !$OMP MASTER + !$ nproc = omp_get_num_threads() + !$OMP END MASTER + !$OMP END PARALLEL END_PROVIDER