mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
Remove overflow and 400% acceleration
This commit is contained in:
Thomas Applencourt
parent
9452013782
commit
2ec05fd08d
@ -430,7 +430,6 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
|
||||
do lambdap=0,lmax+ntotB
|
||||
do k=1,kmax
|
||||
do l=0,lmax
|
||||
|
||||
array_R(ktot,k,l,0,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
|
||||
enddo
|
||||
enddo
|
||||
@ -1827,58 +1826,11 @@ double precision function coef_nk(n,k)
|
||||
|
||||
double precision gam,dble_fact,fact
|
||||
|
||||
double precision, save :: store_result(0:2,0:10)
|
||||
|
||||
store_result(0, 0) = 1.00000000000000d0
|
||||
store_result(0, 1) = 0.166666666666667d0
|
||||
store_result(0, 2) = 8.333333333333333d-3
|
||||
store_result(0, 3) = 1.984126984126984d-4
|
||||
store_result(0, 4) = 2.755731922398589d-6
|
||||
store_result(0, 5) = 2.505210838544172d-8
|
||||
store_result(0, 6) = 1.605904383682161d-10
|
||||
store_result(0, 7) = 7.647163731819816d-13
|
||||
store_result(0, 8) = 2.811457254345521d-15
|
||||
store_result(0, 9) = 8.220635246624329d-018
|
||||
store_result(0,10) = 1.957294106339126d-020
|
||||
|
||||
store_result(1, 0) = 0.333333333333333d0
|
||||
store_result(1, 1) = 3.333333333333333d-002
|
||||
store_result(1, 2) = 1.190476190476191d-003
|
||||
store_result(1, 3) = 2.204585537918871d-005
|
||||
store_result(1, 4) = 2.505210838544172d-007
|
||||
store_result(1, 5) = 1.927085260418594d-009
|
||||
store_result(1, 6) = 1.070602922454774d-011
|
||||
store_result(1, 7) = 4.498331606952833d-014
|
||||
store_result(1, 8) = 1.479714344392379d-016
|
||||
store_result(1, 9) = 3.914588212678252d-019
|
||||
store_result(1,10) = 8.509974375387505d-022
|
||||
|
||||
store_result(2, 0) = 6.666666666666667d-002
|
||||
store_result(2, 1) = 4.761904761904762d-003
|
||||
store_result(2, 2) = 1.322751322751323d-004
|
||||
store_result(2, 3) = 2.004168670835337d-006
|
||||
store_result(2, 4) = 1.927085260418594d-008
|
||||
store_result(2, 5) = 1.284723506945729d-010
|
||||
store_result(2, 6) = 6.297664249733966d-013
|
||||
store_result(2, 7) = 2.367542951027807d-015
|
||||
store_result(2, 8) = 7.046258782820854d-018
|
||||
store_result(2, 9) = 1.701994875077501d-020
|
||||
store_result(2,10) = 3.403989750155002d-023
|
||||
|
||||
if (n.LE.2) then
|
||||
if (k.LE.10) then
|
||||
coef_nk = store_result(n,k)
|
||||
return
|
||||
endif
|
||||
endif
|
||||
|
||||
if ((n+k).GE.60) then
|
||||
coef_nk = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
gam=dble_fact(2*(n+k)+1)
|
||||
coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
|
||||
|
||||
! coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
|
||||
|
||||
coef_nk=1.d0/(2.d0**k*fact(k)*gam)
|
||||
|
||||
return
|
||||
|
||||
@ -1897,10 +1849,11 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
|
||||
integer n,k,m,q,l,kcp
|
||||
double precision gam,dble_fact,fact,pi,a,b
|
||||
double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
|
||||
double precision dump
|
||||
|
||||
double precision, allocatable :: temp_array_a(:)
|
||||
double precision, allocatable :: temp_array_b(:)
|
||||
integer :: n_1,m_1,nlm
|
||||
double precision :: term_A, term_B, term_rap, expo
|
||||
double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square
|
||||
double precision :: int_prod_bessel_loc
|
||||
|
||||
logical done
|
||||
|
||||
@ -1928,41 +1881,58 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
|
||||
intold=-1.d0
|
||||
int=0.d0
|
||||
done=.false.
|
||||
kcp=0
|
||||
|
||||
allocate( temp_array_a(0:100))
|
||||
allocate( temp_array_b(0+m:200+m))
|
||||
n_1 = 2*(n)+1
|
||||
m_1 = 2*m+1
|
||||
nlm = n+m+l
|
||||
pi=dacos(-1.d0)
|
||||
a_over_b_square = (a/b)**2
|
||||
|
||||
! Calcul first term of the sequence
|
||||
|
||||
term_a =dble_fact(nlm-1) / (dble_fact(n_1)*dble_fact(m_1))
|
||||
expo=0.5d0*dfloat(nlm+1)
|
||||
term_rap = term_a / (2.d0*gam)**expo
|
||||
|
||||
s_0_0=term_rap*a**(n)*b**(m)
|
||||
if(mod(nlm,2).eq.0)s_0_0=s_0_0*dsqrt(pi/2.d0)
|
||||
|
||||
! Initialise the first recurence terme for the q loop
|
||||
s_q_0 = s_0_0
|
||||
|
||||
! Loop over q for the convergence of the sequence
|
||||
do while (.not.done)
|
||||
|
||||
kcp=kcp+1
|
||||
sum=0.d0
|
||||
|
||||
temp_array_a(q) = coef_nk(n,q)*(a**(n+2*q))
|
||||
temp_array_b(q) = coef_nk(m,q)*b**(m+2*q)
|
||||
! Init
|
||||
sum=0
|
||||
s_q_k=s_q_0
|
||||
|
||||
! Iteration of k
|
||||
do k=0,q
|
||||
sum=sum+temp_array_a(k)*temp_array_b(q-2*k)
|
||||
sum=sum+s_q_k
|
||||
s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)/dble(2*(k+n)+3) ) * ( dble(q-k)/dble(k+1)) * s_q_k
|
||||
enddo
|
||||
|
||||
int=int+sum*crochet(2*q+n+m+l,gam)
|
||||
int=int+sum
|
||||
|
||||
if(dabs(int-intold).lt.1d-15)then
|
||||
done=.true.
|
||||
if(dabs(int-intold).lt.1d-15)then
|
||||
done=.true.
|
||||
else
|
||||
|
||||
!Compute the s_q+1_0
|
||||
s_q_0=s_q_0*(2.d0*q+nlm+1)*b**2/((2.d0*(m+q)+3)*4.d0*(q+1)*gam)
|
||||
|
||||
if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi/2.d0)
|
||||
! Increment q
|
||||
q=q+1
|
||||
intold=int
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
deallocate(temp_array_a)
|
||||
deallocate(temp_array_b)
|
||||
|
||||
int_prod_bessel=int*freal
|
||||
|
||||
if(kcp.gt.100)then
|
||||
if(q.gt.100)then
|
||||
print*,'**WARNING** bad convergence in int_prod_bessel'
|
||||
endif
|
||||
|
||||
@ -1971,61 +1941,15 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
|
||||
|
||||
if(a.eq.0.d0.and.b.ne.0.d0)then
|
||||
|
||||
if(n.ne.0)then
|
||||
int_prod_bessel=0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
q=0
|
||||
intold=-1.d0
|
||||
int=0.d0
|
||||
done=.false.
|
||||
kcp=0
|
||||
|
||||
do while (.not.done)
|
||||
|
||||
kcp=kcp+1
|
||||
int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
|
||||
|
||||
if(dabs(int-intold).lt.1d-15)then
|
||||
done=.true.
|
||||
else
|
||||
q=q+1
|
||||
intold=int
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
int = int_prod_bessel_loc(l,gam,m,b)
|
||||
int_prod_bessel=int*freal
|
||||
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
|
||||
return
|
||||
endif
|
||||
|
||||
if(a.ne.0.d0.and.b.eq.0.d0)then
|
||||
if(m.ne.0)then
|
||||
int_prod_bessel=0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
q=0
|
||||
intold=-1.d0
|
||||
int=0.d0
|
||||
done=.false.
|
||||
kcp=0
|
||||
|
||||
do while (.not.done)
|
||||
kcp=kcp+1
|
||||
int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
|
||||
if(dabs(int-intold).lt.1d-15)then
|
||||
done=.true.
|
||||
else
|
||||
q=q+1
|
||||
intold=int
|
||||
endif
|
||||
enddo
|
||||
|
||||
int = int_prod_bessel_loc(l,gam,n,a)
|
||||
int_prod_bessel=int*freal
|
||||
if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
|
||||
return
|
||||
|
||||
endif
|
||||
@ -2101,30 +2025,49 @@ end
|
||||
!!
|
||||
!! M_n(x) modified spherical bessel function
|
||||
!!
|
||||
double precision function int_prod_bessel_loc(l,gam,n,a)
|
||||
implicit none
|
||||
integer n,k,l,kcp
|
||||
double precision gam,a
|
||||
double precision int,intold,coef_nk,crochet
|
||||
logical done
|
||||
k=0
|
||||
intold=-1.d0
|
||||
int=0.d0
|
||||
done=.false.
|
||||
kcp=0
|
||||
do while (.not.done)
|
||||
kcp=kcp+1
|
||||
int=int+coef_nk(n,k)*a**(n+2*k)*crochet(2*k+n+l,gam)
|
||||
if(dabs(int-intold).lt.1d-15)then
|
||||
done=.true.
|
||||
else
|
||||
k=k+1
|
||||
intold=int
|
||||
endif
|
||||
enddo
|
||||
int_prod_bessel_loc=int
|
||||
if(kcp.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel'
|
||||
end
|
||||
double precision function int_prod_bessel_loc(l,gam,n,a)
|
||||
implicit none
|
||||
integer n,k,l,kcp
|
||||
double precision gam,a
|
||||
double precision int,intold,coef_nk,crochet,dble_fact, fact, pi, expo
|
||||
double precision :: f_0, f_k
|
||||
logical done
|
||||
|
||||
pi=dacos(-1.d0)
|
||||
intold=-1.d0
|
||||
done=.false.
|
||||
int=0
|
||||
|
||||
! Int f_0
|
||||
coef_nk=1.d0/dble_fact(2*n+1)
|
||||
expo=0.5d0*dfloat(n+l+1)
|
||||
crochet=dble_fact(n+l-1)/(2.d0*gam)**expo
|
||||
if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(pi/2.d0)
|
||||
|
||||
f_0 = coef_nk*a**n*crochet
|
||||
|
||||
k=0
|
||||
|
||||
f_k=f_0
|
||||
do while (.not.done)
|
||||
|
||||
int=int+f_k
|
||||
|
||||
f_k = f_k*(a**2*(2*(k+1)+n+l-1)) / (2*(k+1)*(2*(n+k+1)+1)*2*gam)
|
||||
|
||||
if(dabs(int-intold).lt.1d-15)then
|
||||
done=.true.
|
||||
else
|
||||
k=k+1
|
||||
intold=int
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
int_prod_bessel_loc=int
|
||||
if(k.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel'
|
||||
|
||||
end
|
||||
|
||||
double precision function int_prod_bessel_num(l,gam,n,m,a,b)
|
||||
implicit none
|
||||
|
||||
Loading…
Reference in New Issue
Block a user