mpl_module.f90 2.98 KB
Newer Older
Luca's avatar
Luca committed
1 2
 
MODULE mpl_module
3
  use globals
ulrich_y's avatar
ulrich_y committed
4
  use ieps
5
  implicit none
Luca's avatar
Luca committed
6

ulrich_y's avatar
ulrich_y committed
7 8 9
  real(kind=prec), parameter :: underflowalert = 1.e-250_prec
  real(kind=prec), parameter ::  overflowalert = 1.e+250_prec

ulrich_y's avatar
ulrich_y committed
10 11 12 13 14 15 16 17 18 19 20 21
#ifdef MPL_CACHE
  ! Max weight, no. cache
  integer,parameter :: cachesize(2) = (/ 4, 2500 /)
  type el
    integer :: m(cachesize(1))
    complex(kind=prec) :: c(cachesize(1))
    complex(kind=prec) :: ans
  end type el

  type(el) cache(cachesize(1),abs(cachesize(2)))
  integer cacheind(cachesize(1))
#endif
22 23
CONTAINS 

24
  FUNCTION MPL_converges(m,x)
Luca's avatar
Luca committed
25 26
    ! checks if the MPL converges 
    complex(kind=prec) :: x(:)
27
    integer :: m(:)
Luca's avatar
Luca committed
28
    logical :: MPL_converges
29
    MPL_converges = .false.
Luca's avatar
Luca committed
30
    if(abs(product(x)) < 1) then
31
      if(m(1) /= 1 .or. abs(x(1) - 1) < zero) then
32 33
        MPL_converges = .true.
      end if
Luca's avatar
Luca committed
34
    end if
Luca's avatar
Luca committed
35
  END FUNCTION MPL_converges
Luca's avatar
Luca committed
36

ulrich_y's avatar
ulrich_y committed
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
#ifdef MPL_CACHE
  FUNCTION CHECK_CACHE(m, x, res)
    integer m(:)
    complex(kind=prec) :: x(:), res
    logical check_cache
    integer k, j, i
    j = size(x)
    if (j<=4) then
      do k=1,cacheind(j)
        if (all(cache(j,k)%m(:j) == m)) then
          do i=1,j
            if( abs(cache(j,k)%c(i)-x(i)).gt.zero ) goto 123
          enddo
          res = cache(j,k) % ans
          check_cache = .true.
          return
        endif
123     continue
      enddo
    endif
    check_cache = .false.
  END FUNCTION
#endif



ulrich_y's avatar
ulrich_y committed
63
  FUNCTION MPL(m, x) result(res)
64
    integer :: m(:)
Luca's avatar
Luca committed
65 66
    complex(kind=prec) :: x(:)
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
67 68
    complex(kind=prec) :: t(size(x)), cpow
    integer(kind=ikin) :: q, j, k, ipow
ulrich_y's avatar
ulrich_y committed
69 70 71 72
#ifdef MPL_CACHE
    if (check_cache(m,x,res)) return
#endif

ulrich_y's avatar
ulrich_y committed
73
    j = size(x)
74

Luca's avatar
Luca committed
75

76 77 78
#ifdef DEBUG
    if(verb >= 70) print*, 'called MPL(',m,',',x,')'
#endif
79 80 81
    if(size(m) /= size(x)) then
      print*, 'Error: m and x must have the same length'
    end if
ulrich_y's avatar
ulrich_y committed
82 83 84 85 86 87
    res=0
    q=0
    t=0
    do while(.true.)
      res = t(1)
      q = q+1
88
      if (q < 0) exit
ulrich_y's avatar
ulrich_y committed
89 90 91 92 93 94 95 96 97

      cpow = x(j)**q
      ipow = q**m(j)

      if(ipow .eq. 0) exit
      if (abs(cpow).lt.underflowalert) exit
      if (abs(cpow).gt. overflowalert) exit

      t(j) = t(j) + cpow/ipow
ulrich_y's avatar
ulrich_y committed
98
      do k=1,j-1
ulrich_y's avatar
ulrich_y committed
99 100 101 102 103 104
        ipow = (k+q)**m(j-k)
        cpow = x(j-k)**(k+q)
        if(ipow .eq. 0) exit
        if (abs(cpow).lt.underflowalert) exit
        if (abs(cpow).gt. overflowalert) exit
        t(j-k) = t(j-k) + t(j-k+1) * cpow/ipow
ulrich_y's avatar
ulrich_y committed
105 106
      enddo

107
      if (mod(q,2_ikin) .eq. 1) then
ulrich_y's avatar
ulrich_y committed
108
        if (abs(t(1)-res).lt.MPLdelta/100.) exit
ulrich_y's avatar
ulrich_y committed
109 110 111 112
      endif
    enddo
    res = t(1)

Luca Naterop's avatar
Luca Naterop committed
113

ulrich_y's avatar
ulrich_y committed
114 115 116 117 118 119 120
#ifdef MPL_CACHE
    if (j<=cachesize(1)) then
      if (cacheind(j)+1 > abs(cachesize(2))) then
        if (cachesize(2) < 0) print*,"MPL cache for depth=",j," is full. Try increasing cachesize(2)"
        return
      endif
      cacheind(j) = cacheind(j) + 1
ulrich_y's avatar
ulrich_y committed
121 122
      cache(j,cacheind(j)) = el( reshape(m, (/cachesize(1)/), pad=[      0            ]), &
                                 reshape(x, (/cachesize(1)/), pad=[cmplx(0.,kind=prec)]), &
ulrich_y's avatar
ulrich_y committed
123 124 125
                                 res )
    endif
#endif
Luca's avatar
Luca committed
126
  END FUNCTION MPL
Luca's avatar
Luca committed
127

Luca's avatar
Luca committed
128
END MODULE mpl_module