src/quadratures.f90

   1 MODULE quadratures
   2
   3   IMPLICIT none
   4
   5   integer(kind=8) :: subintervals = 100
   6
   7   INTERFACE
   8      FUNCTION integrate(ibeg, iend, myfun, p) result(val)
   9        IMPLICIT none
  10        ! beginning of integration interval
  11        real(kind=8), intent(in) :: ibeg
  12        ! ending of integration interval
  13        real(kind=8), intent(in) :: iend
  14        ! function to be integrated
  15        procedure(funint) :: myfun
  16        ! polynomial order
  17        integer(kind=4), intent(in) :: p
  18        ! result of integration
  19        real(kind=8) :: val
  20      END FUNCTION integrate
  21   END INTERFACE
  22
  23   INTERFACE
  24      FUNCTION quadrature(qbeg, qend, fun) result(val)
  25        IMPLICIT none
  26        real(kind=8), intent(in) :: qbeg, qend
  27        procedure(funint) :: fun
  28        real(kind=8) :: val
  29      END FUNCTION quadrature
  30   END INTERFACE
  31
  32   INTERFACE
  33      FUNCTION funint(x) result(y)
  34        IMPLICIT none
  35        real(kind=8), intent(in) :: x
  36        real(kind=8) :: y
  37      END FUNCTION funint
  38   END INTERFACE
  39
  40 CONTAINS
  41
  42   FUNCTION integrate_generic(ibeg, iend, fun, quad) result(val)
  43     IMPLICIT none
  44     real(kind=8), intent(in) :: ibeg
  45     real(kind=8), intent(in) :: iend
  46     procedure(funint) :: fun
  47     procedure(quadrature) :: quad
  48
  49     real(kind=8) :: val, subinterval_width, qbeg, qend, partsums(48)
  50     logical :: partsums_mask(48) = .true.
  51     real(kind=8), allocatable :: partval[:]
  52
  53     integer(kind=8) :: min_i, max_i, i, j, subintervals_per_thread
  54     integer(kind=4) :: im
  55
  56     allocate(partval[*])
  57
  58     subintervals_per_thread =                               &
  59          (subintervals + num_images() - 1) / num_images()
  60
  61     min_i = subintervals_per_thread * (this_image() - 1) + 1
  62     max_i = min(subintervals, subintervals_per_thread * this_image())
  63
  64     subinterval_width = (iend - ibeg) / subintervals
  65
  66     ! compute integral using quadrature pointed by quad
  67     partval = 0
  68
  69     DO i = min_i, max_i
  70        qend = ibeg + i * subinterval_width
  71        qbeg = ibeg + (i - 1) * subinterval_width
  72        val = quad(qbeg, qend, fun)
  73
  74        DO j = 1, 48
  75           IF (partsums_mask(j)) THEN
  76              partsums_mask(j) = .false.
  77              partsums(j) = val
  78              EXIT
  79           END IF
  80
  81           val = val + partsums(j)
  82           partsums_mask(j) = .true.
  83        END DO
  84     END DO
  85
  86     partval = 0
  87     DO j = 1, 48
  88        IF (.not. partsums_mask(j)) partval = partval + partsums(j)
  89     END DO
  90
  91     IF (this_image() == 1 .and. num_images() > 1) THEN
  92        sync images([(im, im = 2, num_images())])
  93        partval = partval + sum([(partval[im], im = 2, num_images())])
  94        sync images([(im, im = 2, num_images())])
  95     END IF
  96
  97     IF (this_image() /= 1) THEN
  98        sync images(1)
  99        sync images(1)
 100     END IF
 101
 102     val = partval[1]
 103
 104   END FUNCTION integrate_generic
 105
 106   FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
 107     IMPLICIT none
 108     real(kind=8), intent(in) :: ibeg
 109     real(kind=8), intent(in) :: iend
 110     procedure(funint) :: fun
 111     integer(kind=4), intent(in) :: p
 112     real(kind=8) :: val
 113
 114     procedure(quadrature), pointer :: quad
 115
 116     SELECT CASE (p)
 117     CASE (:-1)
 118        STOP "negative interpolationg polynomial order passed"
 119     CASE (0)
 120        quad => rectangle
 121     CASE (1)
 122        quad => trapeze
 123     CASE (2)
 124        quad => simpson_1_third
 125     CASE default
 126        STOP "Newton-Cotes quadratures only implemented for order < 3"
 127     END SELECT
 128
 129     val = integrate_generic(ibeg, iend, fun, quad)
 130   END FUNCTION newton_cotes
 131
 132   FUNCTION rectangle(qbeg, qend, fun) result(val)
 133     IMPLICIT none
 134     real(kind=8), intent(in) :: qbeg, qend
 135     procedure(funint) :: fun
 136     real(kind=8) :: val
 137
 138     val = (qend - qbeg) * fun((qend + qbeg) * 0.5)
 139   END FUNCTION rectangle
 140
 141   FUNCTION trapeze(qbeg, qend, fun) result(val)
 142     IMPLICIT none
 143     real(kind=8), intent(in) :: qbeg, qend
 144     procedure(funint) :: fun
 145     real(kind=8) :: val
 146
 147     val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend))
 148   END FUNCTION trapeze
 149
 150   FUNCTION simpson_1_third(qbeg, qend, fun) result(val)
 151     IMPLICIT none
 152     real(kind=8), intent(in) :: qbeg, qend
 153     procedure(funint) :: fun
 154     real(kind=8) :: val
 155
 156     val = (qend - qbeg) * (1/6.0) *                               &
 157          (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
 158   END FUNCTION simpson_1_third
 159
 160   FUNCTION gauss(ibeg, iend, fun, p) result(val)
 161     IMPLICIT none
 162     real(kind=8), intent(in) :: ibeg
 163     real(kind=8), intent(in) :: iend
 164     procedure(funint) :: fun
 165     integer(kind=4), intent(in) :: p
 166     real(kind=8) :: val
 167
 168     procedure(quadrature), pointer :: quad
 169
 170     SELECT CASE (p)
 171     CASE (:0)
 172        STOP "non-positive Legendre polynomial order passed"
 173     CASE (1)
 174        quad => gauss_n1
 175     CASE (2)
 176        quad => gauss_n2
 177     CASE (3)
 178        quad => gauss_n3
 179     CASE default
 180        STOP "Gauss quadratures only implemented with roots of" &
 181             // " Legendgre polynomial of order <= 3"
 182     END SELECT
 183
 184     val = integrate_generic(ibeg, iend, fun, quad)
 185   END FUNCTION gauss
 186
 187   FUNCTION gauss_generic(mid, halfwidth, fun, points, weights)    &
 188        result(val)
 189     IMPLICIT none
 190     real(kind=8), intent(in) :: mid, halfwidth, points(:), weights(:)
 191     procedure(funint) :: fun
 192     real(kind=8) :: val
 193
 194     integer(kind=4) :: i
 195
 196     val = halfwidth * sum(weights * &
 197          [(fun(points(i) * halfwidth + mid), i = 1, size(points))])
 198   END FUNCTION gauss_generic
 199
 200   FUNCTION gauss_n1(qbeg, qend, fun) result(val)
 201     IMPLICIT none
 202     real(kind=8), intent(in) :: qbeg, qend
 203     procedure(funint) :: fun
 204     real(kind=8) :: val, weights(1) = [2], points(1) = [0]
 205
 206     val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
 207          fun, points, weights)
 208   END FUNCTION gauss_n1
 209
 210   FUNCTION gauss_n2(qbeg, qend, fun) result(val)
 211     IMPLICIT none
 212     real(kind=8), intent(in) :: qbeg, qend
 213     procedure(funint) :: fun
 214     real(kind=8) :: val, weights(2) = [1, 1],                     &
 215          points(2) = [1 / sqrt(3.0), -1 / sqrt(3.0)]
 216
 217     val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
 218          fun, points, weights)
 219   END FUNCTION gauss_n2
 220
 221   FUNCTION gauss_n3(qbeg, qend, fun) result(val)
 222     IMPLICIT none
 223     real(kind=8), intent(in) :: qbeg, qend
 224     procedure(funint) :: fun
 225     real(kind=8) :: val, weights(3) = [8 / 9.0, 5 / 9.0, 5 / 9.0],&
 226          points(3) = [0.0, sqrt(3 / 5.0), -sqrt(3 / 5.0)]
 227
 228     val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
 229          fun, points, weights)
 230   END FUNCTION gauss_n3
 231
 232 END MODULE quadratures