src/quadratures.f90

   1 ! Copyright 2019 Wojciech Kosior
   2
   3 ! This is free and unencumbered software released into the public domain.
   4
   5 ! Anyone is free to copy, modify, publish, use, compile, sell, or
   6 ! distribute this software, either in source code form or as a compiled
   7 ! binary, for any purpose, commercial or non-commercial, and by any
   8 ! means.
   9
  10 ! In jurisdictions that recognize copyright laws, the author or authors
  11 ! of this software dedicate any and all copyright interest in the
  12 ! software to the public domain. We make this dedication for the benefit
  13 ! of the public at large and to the detriment of our heirs and
  14 ! successors. We intend this dedication to be an overt act of
  15 ! relinquishment in perpetuity of all present and future rights to this
  16 ! software under copyright law.
  17
  18 ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  19 ! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  20 ! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  21 ! IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  22 ! OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  23 ! ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24 ! OTHER DEALINGS IN THE SOFTWARE.
  25
  26 ! For more information, please refer to <http://unlicense.org/>
  27
  28 MODULE quadratures
  29
  30   IMPLICIT none
  31
  32   integer(kind=8) :: subintervals = 100
  33
  34   INTERFACE
  35      FUNCTION integrate(ibeg, iend, myfun, p) result(val)
  36        IMPLICIT none
  37        ! beginning of integration interval
  38        real(kind=8), intent(in) :: ibeg
  39        ! ending of integration interval
  40        real(kind=8), intent(in) :: iend
  41        ! function to be integrated
  42        procedure(funint) :: myfun
  43        ! polynomial order
  44        integer(kind=4), intent(in) :: p
  45        ! result of integration
  46        real(kind=8) :: val
  47      END FUNCTION integrate
  48   END INTERFACE
  49
  50   INTERFACE
  51      FUNCTION quadrature(qbeg, qend, fun) result(val)
  52        IMPLICIT none
  53        real(kind=8), intent(in) :: qbeg, qend
  54        procedure(funint) :: fun
  55        real(kind=8) :: val
  56      END FUNCTION quadrature
  57   END INTERFACE
  58
  59   INTERFACE
  60      FUNCTION funint(x) result(y)
  61        IMPLICIT none
  62        real(kind=8), intent(in) :: x
  63        real(kind=8) :: y
  64      END FUNCTION funint
  65   END INTERFACE
  66
  67 CONTAINS
  68
  69   FUNCTION integrate_generic(ibeg, iend, fun, quad) result(val)
  70     IMPLICIT none
  71     real(kind=8), intent(in) :: ibeg
  72     real(kind=8), intent(in) :: iend
  73     procedure(funint) :: fun
  74     procedure(quadrature) :: quad
  75
  76     real(kind=8) :: val, subinterval_width, qbeg, qend, partsums(48)
  77     logical :: partsums_mask(48) = .true.
  78     real(kind=8), allocatable :: partval[:]
  79
  80     integer(kind=8) :: min_i, max_i, i, j, subintervals_per_thread
  81     integer(kind=4) :: im
  82
  83     allocate(partval[*])
  84
  85     subintervals_per_thread =                               &
  86          (subintervals + num_images() - 1) / num_images()
  87
  88     min_i = subintervals_per_thread * (this_image() - 1) + 1
  89     max_i = min(subintervals, subintervals_per_thread * this_image())
  90
  91     subinterval_width = (iend - ibeg) / subintervals
  92
  93     ! compute integral using quadrature pointed by quad
  94     partval = 0
  95
  96     DO i = min_i, max_i
  97        qend = ibeg + i * subinterval_width
  98        qbeg = ibeg + (i - 1) * subinterval_width
  99        val = quad(qbeg, qend, fun)
 100
 101        DO j = 1, 48
 102           IF (partsums_mask(j)) THEN
 103              partsums_mask(j) = .false.
 104              partsums(j) = val
 105              EXIT
 106           END IF
 107
 108           val = val + partsums(j)
 109           partsums_mask(j) = .true.
 110        END DO
 111     END DO
 112
 113     partval = 0
 114     DO j = 1, 48
 115        IF (.not. partsums_mask(j)) partval = partval + partsums(j)
 116     END DO
 117
 118     IF (this_image() == 1 .and. num_images() > 1) THEN
 119        sync images([(im, im = 2, num_images())])
 120        partval = partval + sum([(partval[im], im = 2, num_images())])
 121        sync images([(im, im = 2, num_images())])
 122     END IF
 123
 124     IF (this_image() /= 1) THEN
 125        sync images(1)
 126        sync images(1)
 127     END IF
 128
 129     val = partval[1]
 130
 131   END FUNCTION integrate_generic
 132
 133   FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
 134     USE, intrinsic :: ieee_arithmetic
 135
 136     IMPLICIT none
 137     real(kind=8), intent(in) :: ibeg
 138     real(kind=8), intent(in) :: iend
 139     procedure(funint) :: fun
 140     integer(kind=4), intent(in) :: p
 141     real(kind=8) :: val
 142
 143     procedure(quadrature), pointer :: quad
 144
 145     SELECT CASE (p)
 146     CASE (0)
 147        quad => rectangle
 148     CASE (1)
 149        quad => trapeze
 150     CASE (2)
 151        quad => simpson_1_third
 152     CASE default
 153        val = ieee_value(val, ieee_quiet_nan)
 154        RETURN
 155     END SELECT
 156
 157     val = integrate_generic(ibeg, iend, fun, quad)
 158   END FUNCTION newton_cotes
 159
 160   FUNCTION rectangle(qbeg, qend, fun) result(val)
 161     IMPLICIT none
 162     real(kind=8), intent(in) :: qbeg, qend
 163     procedure(funint) :: fun
 164     real(kind=8) :: val
 165
 166     val = (qend - qbeg) * fun((qend + qbeg) * 0.5)
 167   END FUNCTION rectangle
 168
 169   FUNCTION trapeze(qbeg, qend, fun) result(val)
 170     IMPLICIT none
 171     real(kind=8), intent(in) :: qbeg, qend
 172     procedure(funint) :: fun
 173     real(kind=8) :: val
 174
 175     val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend))
 176   END FUNCTION trapeze
 177
 178   FUNCTION simpson_1_third(qbeg, qend, fun) result(val)
 179     IMPLICIT none
 180     real(kind=8), intent(in) :: qbeg, qend
 181     procedure(funint) :: fun
 182     real(kind=8) :: val
 183
 184     val = (qend - qbeg) * (1/6.0) *                               &
 185          (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
 186   END FUNCTION simpson_1_third
 187
 188   FUNCTION gauss(ibeg, iend, fun, p) result(val)
 189     USE, intrinsic :: ieee_arithmetic
 190
 191     IMPLICIT none
 192     real(kind=8), intent(in) :: ibeg
 193     real(kind=8), intent(in) :: iend
 194     procedure(funint) :: fun
 195     integer(kind=4), intent(in) :: p
 196     real(kind=8) :: val
 197
 198     procedure(quadrature), pointer :: quad
 199
 200     SELECT CASE (p)
 201     CASE (0)
 202        quad => gauss_n1
 203     CASE (1)
 204        quad => gauss_n2
 205     CASE (2)
 206        quad => gauss_n3
 207     CASE default
 208        val = ieee_value(val, ieee_quiet_nan)
 209        RETURN
 210     END SELECT
 211
 212     val = integrate_generic(ibeg, iend, fun, quad)
 213   END FUNCTION gauss
 214
 215   FUNCTION gauss_generic(mid, halfwidth, fun, points, weights)    &
 216        result(val)
 217     IMPLICIT none
 218     real(kind=8), intent(in) :: mid, halfwidth, points(:), weights(:)
 219     procedure(funint) :: fun
 220     real(kind=8) :: val
 221
 222     integer(kind=4) :: i
 223
 224     val = halfwidth * sum(weights * &
 225          [(fun(points(i) * halfwidth + mid), i = 1, size(points))])
 226   END FUNCTION gauss_generic
 227
 228   FUNCTION gauss_n1(qbeg, qend, fun) result(val)
 229     IMPLICIT none
 230     real(kind=8), intent(in) :: qbeg, qend
 231     procedure(funint) :: fun
 232     real(kind=8) :: val, weights(1) = [2], points(1) = [0]
 233
 234     val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
 235          fun, points, weights)
 236   END FUNCTION gauss_n1
 237
 238   FUNCTION gauss_n2(qbeg, qend, fun) result(val)
 239     IMPLICIT none
 240     real(kind=8), intent(in) :: qbeg, qend
 241     procedure(funint) :: fun
 242     real(kind=8) :: val, weights(2) = [1, 1],                     &
 243          points(2) = [1 / sqrt(3.0), -1 / sqrt(3.0)]
 244
 245     val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
 246          fun, points, weights)
 247   END FUNCTION gauss_n2
 248
 249   FUNCTION gauss_n3(qbeg, qend, fun) result(val)
 250     IMPLICIT none
 251     real(kind=8), intent(in) :: qbeg, qend
 252     procedure(funint) :: fun
 253     real(kind=8) :: val, weights(3) = [8 / 9.0, 5 / 9.0, 5 / 9.0],&
 254          points(3) = [0.0, sqrt(3 / 5.0), -sqrt(3 / 5.0)]
 255
 256     val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
 257          fun, points, weights)
 258   END FUNCTION gauss_n3
 259
 260 END MODULE quadratures