AGH-fortran-course/fortran-assignment3 - This repo contains my homework for the fortran course at AGH University of Science and Technology

! 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 <http://unlicense.org/>

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