Commit db028962 authored by ulrich_y's avatar ulrich_y

Added polylog cache

parent 9ff6a342
......@@ -28,6 +28,10 @@ ifeq ($(HAVE_GINAC),1)
FFLAGS += -DHAVE_GINAC
endif
ifeq ($(DISABLE_CACHE),1)
FFLAGS += -DNOCACHE
endif
else
FFLAGS=-autodouble -module build -fpp -stand f03 -O3 -xHost -fast
endif
......
......@@ -2,6 +2,9 @@ MODULE GPL
use globals, only: prec, GPLopts=>set_options
use ieps, only: inum, toinum, di0, abs, real, aimag
use gpl_module, only: G
#ifndef NOCACHE
use maths_functions, only: clearcache
#endif
implicit none
END MODULE GPL
......@@ -11,6 +11,8 @@ MODULE globals
real(kind=prec), protected :: MPLdel = 1e-15 ! if the MPL sum changes less then del it is truncated.
integer, protected :: PolylogInfinity = 1000 ! expansion order for Polylogs
real(kind=prec), protected :: HoelderCircle = 1.1 ! when to apply Hoelder convolution?
integer, parameter :: PolyLogCacheSize(2) = (/ 5, 100 /)
! = (/ mmax, n /). At most n polylogs with weight mmax will be cached
integer :: verb = 0
......
......@@ -11,6 +11,13 @@ MODULE maths_functions
1.03692775514337, 1.0173430619844488, 1.008349277381923, &
1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
type el
type(inum) :: c
complex(kind=prec) ans
end type el
type(el) :: cache(PolyLogCacheSize(1),PolyLogCacheSize(2))
integer :: plcachesize(PolyLogCacheSize(1)) = 0
CONTAINS
FUNCTION naive_polylog(m,x) result(res)
......@@ -292,30 +299,46 @@ CONTAINS
integer :: m
type(inum) :: x, inv
complex(kind=prec) :: res
integer i
#ifdef DEBUG
if(verb >= 70) print*, 'called polylog(',m,',',x%c,x%i0,')'
#endif
#ifndef NOCACHE
if (m.le.5) then
do i=1,plcachesize(m)
if( abs(cache(m,i)%c%c-x%c).lt.zero ) then
res = cache(m,i)%ans
return
endif
enddo
endif
#endif
if ((m.le.9).and.(abs(x%c-1.).lt.zero)) then
res = zeta(m)
return
else if ((m.le.9).and.(abs(x%c+1.).lt.zero)) then
res = -(1. - 2.**(1-m))*zeta(m)
return
else if (abs(x) .gt. 1) then
inv = inum(1./x%c, x%i0)
res = (-1)**(m-1)*polylog(m,inv) &
- cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*conjg(log(-x%c))/2/pi) / factorial(m)
return
endif
if(m == 2) then
else if(m == 2) then
res = dilog(x%c)
else if(m == 3) then
res = trilog(x%c)
else
res = naive_polylog(m,x%c)
end if
#ifndef NOCACHE
if (m.le.PolyLogCacheSize(1)) then
if (plcachesize(m).lt.PolyLogCacheSize(2)) then
plcachesize(m) = plcachesize(m) + 1
cache(m,plcachesize(m)) = el(x,res)
endif
endif
#endif
END FUNCTION polylog1
......@@ -339,6 +362,12 @@ CONTAINS
plog1 = log(1.-a%c/b%c)
END FUNCTION
#ifndef NOCACHE
SUBROUTINE CLEARCACHE
plcachesize=0
END SUBROUTINE
#endif
END MODULE maths_functions
! PROGRAM test
......
......@@ -1302,9 +1302,11 @@ CONTAINS
endif
end subroutine
subroutine do_muone_tests(x,y,msg)
use maths_functions, only:clearcache
complex(kind=prec) x, y
real(kind=prec) tstart, tend
character(len=*) :: msg
call clearcache
call cpu_time(tstart)
call test_one_ginac([(-1.,0.),x],'6.1')
call test_one_ginac([(-1.,0.),(-1.,0.),x],'6.2')
......@@ -1507,7 +1509,7 @@ CONTAINS
call cpu_time(tend)
write(*,900) msg,198./(tend-tstart)
900 format("Evaluating ",A," at ",F8.2,"G/s")
900 format("Evaluating ",A," at ",F9.2,"G/s")
end subroutine
#endif
! subroutine do_shuffle_tests()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment