From ec34d430dcc100e15fb47349fc0c88239b8287ad Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@users.sourceforge.net>
Date: Tue, 22 Dec 2009 00:07:34 +0100
Subject: [PATCH] 75% acceleration

---
 src/eplf_function.irp.f | 31 +++++++++++++++++++++++--------
 1 file changed, 23 insertions(+), 8 deletions(-)

diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f
index 9d83ab3..857ebc4 100644
--- a/src/eplf_function.irp.f
+++ b/src/eplf_function.irp.f
@@ -282,8 +282,8 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  integer :: ii, jj, kk, ll
  real :: xp1,xp
  real :: p1,p
- double precision :: S00, xpa, xpb
- double precision :: inv_p,c,c1
+ double precision :: xpa, xpb
+ double precision :: inv_p(2),S00, c
  double precision :: ObaraS
 
  ASSERT (a>0.)
@@ -291,15 +291,30 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  ASSERT (i>=0)
  ASSERT (j>=0)
 
- ! Gaussian product
- call gaussian_product(a,xa,b,xb,c1,p1,xp1)
- call gaussian_product(p1,xp1,gmma,xr,c,p,xp)
+ ! Gaussian products
+!double precision :: t(2), xab(2), ab(2)
+ real :: t(2), xab(2), ab(2)
+!call gaussian_product(a,xa,b,xb,c1,p1,xp1)
+ inv_p(1) = 1.d0/(a+b)
+ p1 = a+b
+ ab(1) = a*b
+ inv_p(2) = 1.d0/(p1+gmma)
+ t(1) = (a*xa+b*xb)
+ xab(1) = xa-xb
+ xp1 = t(1)*inv_p(1)
+ p = p1+gmma
+ ab(2) = p1*gmma
+ t(2) = (p1*xp1+gmma*xr)
+ xab(2) = xp1-xr
+ xp = t(2)*inv_p(2)
+ c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
+             ab(2)*inv_p(2)*xab(2)**2) )
+
 
- inv_p = 1.d0/p
- S00 = dsqrt(pi*inv_p)*c*c1
  xpa = xp-xa
  xpb = xp-xb
- ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p,S00)
+ S00 = dsqrt(pi*inv_p(2))*c
+ ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
 
 end function