Commit a2d94600 authored by Luca's avatar Luca

cleanup, test cases, check convergence

parent f33001b6
MODULE functions
implicit none
......@@ -26,9 +27,19 @@ CONTAINS
j = (/(i, i=1,n,1)/)
polylog = sum(x**j / j**m)
END FUNCTION polylog
FUNCTION multiple_polylog_converges(x)
! checks if the MPL converges
complex(kind=prec) :: x(:)
logical :: multiple_polylog_converges
if(abs(product(x)) < 1) then
multiple_polylog_converges = .true.
else
multiple_polylog_converges = .false.
end if
END FUNCTION multiple_polylog_converges
recursive FUNCTION multiple_polylog(m, x, n_passed) result(res)
! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
integer :: m(:)
......@@ -59,5 +70,12 @@ CONTAINS
end if
END FUNCTION multiple_polylog
END MODULE functions
! PROGRAM test
! use functions
! logical :: result
! result = multiple_polylog_converges( dcmplx((/10.1d0,.7d0,.3d0/)) )
! print*, result
! end PROGRAM test
......@@ -3,7 +3,7 @@ UNAME_S := $(shell uname -s)
SHA1=sha1sum
FC=gfortran
FFLAGS=-fdefault-real-8 -fdefault-real-8 -cpp
FFLAGS=-fdefault-real-8 -cpp
LD=gfortran
......@@ -16,7 +16,6 @@ OBJ=functions.o
@$(FC) $(FFLAGS) -c $<
# Rules for linking
test: $(OBJ) test.o
@echo "LD $@"
......
! im terminal kann man den exit code bekommen via echo $?
PROGRAM TEST
use functions
implicit none
complex(kind=prec) :: result, x(4)
! print*, dilog((0.8,0),20) ! should be 1.07479 + 0i
! print*, dilog((0.2,0.5),20) ! should be 0.133909 + 0.537628i
real, parameter :: tol = 1.0e-15
integer :: m(2)
complex(kind=prec) :: x(2)
complex(kind=prec) :: res, ref
! print*, polylog(2,(0.2,0.5),20) ! should be 0.133909 + 0.537628i
! result = polylog(5,(0.2d0,0.5d0),20) ! should be 0.192872345 + 0.505898833i
print*, 'testing multiple polylog...'
m = (/ 1,1 /)
x = cmplx((/ 0.3156498673740053, 0.3431255827785649/))
ref = cmplx(0.022696600480693277651633)
call check_multiple_polylog(m,x,ref)
m = (/ 1,1 /)
x = cmplx((/ 0.3156498673740053, 0.3431255827785649/))
ref = cmplx(0.022696600480693277651633)
call check_multiple_polylog(m,x,ref)
! result = multiple_polylog((/ 5 /),(/ (0.2,0.5) /),10)
! print*, 'result = ', result
! result = multiple_polylog((/ 5, 5 /),(/ (0.8,0),(0.3,0.5) /), 20)
CONTAINS
result = multiple_polylog((/ 0, 3, 4, 5 /),cmplx((/ 0.3,0.8,0.3,0.5 /))) ! 6.929162361968684E-7
subroutine check_multiple_polylog(m,x,ref)
! checks multiple_polylog(m,x) against reference value ref
integer :: m(2)
complex(kind=prec) :: x(2)
complex(kind=prec) :: res, ref
res = multiple_polylog(m,x)
if(abs((res-ref)/ref) < tol) then
print*, 'passed'
else
print*, 'm=',m,'x=',x,'failed'
print*, 'res=',res,'ref=',ref
stop 1
end if
print*, 'result ', result
print*, 'compare', cmplx(6.9291623623242096179E-7)
end subroutine check_multiple_polylog
END PROGRAM TEST
\ No newline at end of file
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