mpl_module.f90 2.4 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 10 11 12 13 14 15 16 17 18
#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
19 20
CONTAINS 

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

ulrich_y's avatar
ulrich_y committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
#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
60
  FUNCTION MPL(m, x) result(res)
61
    integer :: m(:)
Luca's avatar
Luca committed
62 63
    complex(kind=prec) :: x(:)
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
64 65
    complex(kind=prec) :: t(size(x))
    integer :: q, j, k
ulrich_y's avatar
ulrich_y committed
66 67 68 69
#ifdef MPL_CACHE
    if (check_cache(m,x,res)) return
#endif

ulrich_y's avatar
ulrich_y committed
70
    j = size(x)
71

Luca's avatar
Luca committed
72

73 74 75
    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
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
    res=0
    q=0
    t=0
    do while(.true.)
      res = t(1)
      q = q+1
      t(j) = t(j) + x(j)**q/q**m(j)
      do k=1,j-1
        t(j-k) = t(j-k) + t(j-k+1) * x(j-k)**(k+q)/(k+q)**m(j-k)
      enddo

      if (mod(q,2) .eq. 1) then
        if (abs(t(1)-res).lt.MPLdel) exit
      endif
    enddo
    res = t(1)

Luca Naterop's avatar
Luca Naterop committed
93

ulrich_y's avatar
ulrich_y committed
94 95 96 97 98 99 100 101 102 103 104 105
#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
      cache(j,cacheind(j)) = el( reshape(m, (/cachesize(1)/), pad=[      0  ]), &
                                 reshape(x, (/cachesize(1)/), pad=[cmplx(0.)]), &
                                 res )
    endif
#endif
Luca's avatar
Luca committed
106
  END FUNCTION MPL
Luca's avatar
Luca committed
107

Luca's avatar
Luca committed
108
END MODULE mpl_module