! Copyright 2019 Wojciech Kosior ! This is free and unencumbered software released into the public domain. ! Anyone is free to copy, modify, publish, use, compile, sell, or ! distribute this software, either in source code form or as a compiled ! binary, for any purpose, commercial or non-commercial, and by any ! means. ! In jurisdictions that recognize copyright laws, the author or authors ! of this software dedicate any and all copyright interest in the ! software to the public domain. We make this dedication for the benefit ! of the public at large and to the detriment of our heirs and ! successors. We intend this dedication to be an overt act of ! relinquishment in perpetuity of all present and future rights to this ! software under copyright law. ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. ! IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR ! OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ! ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR ! OTHER DEALINGS IN THE SOFTWARE. ! For more information, please refer to MODULE quadratures IMPLICIT none integer(kind=8) :: subintervals = 100 INTERFACE FUNCTION integrate(ibeg, iend, myfun, p) result(val) IMPLICIT none ! beginning of integration interval real(kind=8), intent(in) :: ibeg ! ending of integration interval real(kind=8), intent(in) :: iend ! function to be integrated procedure(funint) :: myfun ! polynomial order integer(kind=4), intent(in) :: p ! result of integration real(kind=8) :: val END FUNCTION integrate END INTERFACE INTERFACE FUNCTION quadrature(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val END FUNCTION quadrature END INTERFACE INTERFACE FUNCTION funint(x) result(y) IMPLICIT none real(kind=8), intent(in) :: x real(kind=8) :: y END FUNCTION funint END INTERFACE CONTAINS FUNCTION integrate_generic(ibeg, iend, fun, quad) result(val) IMPLICIT none real(kind=8), intent(in) :: ibeg real(kind=8), intent(in) :: iend procedure(funint) :: fun procedure(quadrature) :: quad real(kind=8) :: val, subinterval_width, qbeg, qend, partsums(48) logical :: partsums_mask(48) = .true. real(kind=8), allocatable :: partval[:] integer(kind=8) :: min_i, max_i, i, j, subintervals_per_thread integer(kind=4) :: im allocate(partval[*]) subintervals_per_thread = & (subintervals + num_images() - 1) / num_images() min_i = subintervals_per_thread * (this_image() - 1) + 1 max_i = min(subintervals, subintervals_per_thread * this_image()) subinterval_width = (iend - ibeg) / subintervals ! compute integral using quadrature pointed by quad partval = 0 DO i = min_i, max_i qend = ibeg + i * subinterval_width qbeg = ibeg + (i - 1) * subinterval_width val = quad(qbeg, qend, fun) DO j = 1, 48 IF (partsums_mask(j)) THEN partsums_mask(j) = .false. partsums(j) = val EXIT END IF val = val + partsums(j) partsums_mask(j) = .true. END DO END DO partval = 0 DO j = 1, 48 IF (.not. partsums_mask(j)) partval = partval + partsums(j) END DO IF (this_image() == 1 .and. num_images() > 1) THEN sync images([(im, im = 2, num_images())]) partval = partval + sum([(partval[im], im = 2, num_images())]) sync images([(im, im = 2, num_images())]) END IF IF (this_image() /= 1) THEN sync images(1) sync images(1) END IF val = partval[1] END FUNCTION integrate_generic FUNCTION newton_cotes(ibeg, iend, fun, p) result(val) USE, intrinsic :: ieee_arithmetic IMPLICIT none real(kind=8), intent(in) :: ibeg real(kind=8), intent(in) :: iend procedure(funint) :: fun integer(kind=4), intent(in) :: p real(kind=8) :: val procedure(quadrature), pointer :: quad SELECT CASE (p) CASE (0) quad => rectangle CASE (1) quad => trapeze CASE (2) quad => simpson_1_third CASE default val = ieee_value(val, ieee_quiet_nan) RETURN END SELECT val = integrate_generic(ibeg, iend, fun, quad) END FUNCTION newton_cotes FUNCTION rectangle(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val val = (qend - qbeg) * fun((qend + qbeg) * 0.5) END FUNCTION rectangle FUNCTION trapeze(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val val = (qend - qbeg) * 0.5 * (fun(qbeg) + fun(qend)) END FUNCTION trapeze FUNCTION simpson_1_third(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val val = (qend - qbeg) * (1/6.0) * & (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend)) END FUNCTION simpson_1_third FUNCTION gauss(ibeg, iend, fun, p) result(val) USE, intrinsic :: ieee_arithmetic IMPLICIT none real(kind=8), intent(in) :: ibeg real(kind=8), intent(in) :: iend procedure(funint) :: fun integer(kind=4), intent(in) :: p real(kind=8) :: val procedure(quadrature), pointer :: quad SELECT CASE (p) CASE (0) quad => gauss_n1 CASE (1) quad => gauss_n2 CASE (2) quad => gauss_n3 CASE default val = ieee_value(val, ieee_quiet_nan) RETURN END SELECT val = integrate_generic(ibeg, iend, fun, quad) END FUNCTION gauss FUNCTION gauss_generic(mid, halfwidth, fun, points, weights) & result(val) IMPLICIT none real(kind=8), intent(in) :: mid, halfwidth, points(:), weights(:) procedure(funint) :: fun real(kind=8) :: val integer(kind=4) :: i val = halfwidth * sum(weights * & [(fun(points(i) * halfwidth + mid), i = 1, size(points))]) END FUNCTION gauss_generic FUNCTION gauss_n1(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val, weights(1) = [2], points(1) = [0] val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, & fun, points, weights) END FUNCTION gauss_n1 FUNCTION gauss_n2(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val, weights(2) = [1, 1], & points(2) = [1 / sqrt(3.0), -1 / sqrt(3.0)] val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, & fun, points, weights) END FUNCTION gauss_n2 FUNCTION gauss_n3(qbeg, qend, fun) result(val) IMPLICIT none real(kind=8), intent(in) :: qbeg, qend procedure(funint) :: fun real(kind=8) :: val, weights(3) = [8 / 9.0, 5 / 9.0, 5 / 9.0],& points(3) = [0.0, sqrt(3 / 5.0), -sqrt(3 / 5.0)] val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, & fun, points, weights) END FUNCTION gauss_n3 END MODULE quadratures