From f73beda4e5184d7b1ca905219646210ac3863c2a Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Mon, 10 Feb 2020 07:32:53 -0700 Subject: [PATCH 01/15] attempt at quad API --- src/stdlib_experimental_quadrature.fypp | 65 ++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp index ecb15840f..7afdd9262 100644 --- a/src/stdlib_experimental_quadrature.fypp +++ b/src/stdlib_experimental_quadrature.fypp @@ -1,5 +1,6 @@ #:include "common.fypp" - +#:set WEIGHT_FUNS = ["sin", "cos", "pole"] +#:set QUAD_OK = False module stdlib_experimental_quadrature use stdlib_experimental_kinds, only: sp, dp, qp @@ -13,6 +14,15 @@ module stdlib_experimental_quadrature public :: simps public :: simps_weights + ! automatic integration of (weighted) functions + #:if QUAD_OK + public :: quad + public :: weight_t + #:for WFUN in WEIGHT_FUNS + public :: ${WFUN}$_weight_t + #:endfor + #:endif + interface trapz #:for KIND in REAL_KINDS @@ -72,4 +82,57 @@ module stdlib_experimental_quadrature #:endfor end interface simps_weights + + ! Interface for a simple f(x)-style integrand function. + ! Could become fancier as we learn about the performance + ! ramifications of different ways to do callbacks. + abstract interface + #:for KIND in REAL_KINDS + pure function integrand_${KIND}$(x) result(f) + import :: ${KIND}$ + real(${KIND}$), intent(in) :: x + real(${KIND}$) :: f + end function integrand_${KIND}$ + #:endfor + end interface + +#:if QUAD_OK + ! Base class to avoid repeating kind parameter declaration. + type, abstract :: weight_t(kind) + integer, kind :: kind + end type weight_t + + type, extends(weight_t) :: sin_weight_t + real(kind) :: omega + end type sin_weight_t + + type, extends(weight_t) :: cos_weight_t + real(kind) :: omega + end type cos_weight_t + + type, extends(weight_t) :: pole_weight_t + real(kind) :: c + end type pole_weight_t + +! gfortran 9.2.0 chokes on ICE if I include this ("buffer overflow detected") +! Interestingly, though, the ICE happens while trying to build the trapz submodule +! PDT bug? + interface quad + #:for WFUN in WEIGHT_FUNS + #:for KIND in REAL_KINDS + module function quad_${WFUN}$_${KIND}$(f, a, b, weight, points, abstol, reltol, delta) result(integral) + procedure(integrand_${KIND}$) :: f + real(${KIND}$), intent(in) :: a + real(${KIND}$), intent(in) :: b + type(${WFUN}$_weight_t(${KIND}$)), intent(in) :: weight + real(${KIND}$), intent(in), dimension(:) :: points + real(${KIND}$), intent(in) :: abstol + real(${KIND}$), intent(in) :: reltol + real(${KIND}$), intent(out), optional :: delta + real(${KIND}$) :: integral + end function quad_${WFUN}$_${KIND}$ + #:endfor + #:endfor + end interface quad +#:endif end module stdlib_experimental_quadrature From eb3b22f57cdaeee0b923c1638991e26f9f6d84c6 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sun, 31 May 2020 17:04:32 -0600 Subject: [PATCH 02/15] improve fypp usage --- src/stdlib_experimental_quadrature.fypp | 105 ++++++++++++------------ 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp index 7afdd9262..1db6bcb52 100644 --- a/src/stdlib_experimental_quadrature.fypp +++ b/src/stdlib_experimental_quadrature.fypp @@ -25,60 +25,61 @@ module stdlib_experimental_quadrature interface trapz - #:for KIND in REAL_KINDS - pure module function trapz_dx_${KIND}$(y, dx) result(integral) - real(${KIND}$), dimension(:), intent(in) :: y - real(${KIND}$), intent(in) :: dx - real(${KIND}$) :: integral - end function trapz_dx_${KIND}$ + #:for k1, t1 in REAL_KINDS_TYPES + pure module function trapz_dx_${k1}$(y, dx) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, intent(in) :: dx + ${t1}$ :: integral + end function trapz_dx_${k1}$ #:endfor - #:for KIND in REAL_KINDS - pure module function trapz_x_${KIND}$(y, x) result(integral) - real(${KIND}$), dimension(:), intent(in) :: y - real(${KIND}$), dimension(size(y)), intent(in) :: x - real(${KIND}$) :: integral - end function trapz_x_${KIND}$ + #:for k1, t1 in REAL_KINDS_TYPES + pure module function trapz_x_${k1}$(y, x) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, dimension(size(y)), intent(in) :: x + ${t1}$ :: integral + end function trapz_x_${k1}$ #:endfor end interface trapz interface trapz_weights - #:for KIND in REAL_KINDS - pure module function trapz_weights_${KIND}$(x) result(w) - real(${KIND}$), dimension(:), intent(in) :: x - real(${KIND}$), dimension(size(x)) :: w - end function trapz_weights_${KIND}$ + #:for k1, t1 in REAL_KINDS_TYPES + pure module function trapz_weights_${k1}$(x) result(w) + ${t1}$, dimension(:), intent(in) :: x + ${t1}$, dimension(size(x)) :: w + end function trapz_weights_${k1}$ #:endfor end interface trapz_weights interface simps - #:for KIND in REAL_KINDS - pure module function simps_dx_${KIND}$(y, dx, even) result(integral) - real(${KIND}$), dimension(:), intent(in) :: y - real(${KIND}$), intent(in) :: dx + ! "recursive" is an implementation detail + #:for k1, t1 in REAL_KINDS_TYPES + pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, intent(in) :: dx integer, intent(in), optional :: even - real(${KIND}$) :: integral - end function simps_dx_${KIND}$ + ${t1}$ :: integral + end function simps_dx_${k1}$ #:endfor - #:for KIND in REAL_KINDS - pure module function simps_x_${KIND}$(y, x, even) result(integral) - real(${KIND}$), dimension(:), intent(in) :: y - real(${KIND}$), dimension(size(y)), intent(in) :: x + #:for k1, t1 in REAL_KINDS_TYPES + pure recursive module function simps_x_${k1}$(y, x, even) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, dimension(size(y)), intent(in) :: x integer, intent(in), optional :: even - real(${KIND}$) :: integral - end function simps_x_${KIND}$ + ${t1}$ :: integral + end function simps_x_${k1}$ #:endfor end interface simps interface simps_weights - #:for KIND in REAL_KINDS - pure module function simps_weights_${KIND}$(x, even) result(w) - real(${KIND}$), dimension(:), intent(in) :: x - real(${KIND}$), dimension(size(x)) :: w + #:for k1, t1 in REAL_KINDS_TYPES + pure recursive module function simps_weights_${k1}$(x, even) result(w) + ${t1}$, dimension(:), intent(in) :: x integer, intent(in), optional :: even - end function simps_weights_${KIND}$ + ${t1}$, dimension(size(x)) :: w + end function simps_weights_${k1}$ #:endfor end interface simps_weights @@ -87,12 +88,12 @@ module stdlib_experimental_quadrature ! Could become fancier as we learn about the performance ! ramifications of different ways to do callbacks. abstract interface - #:for KIND in REAL_KINDS - pure function integrand_${KIND}$(x) result(f) - import :: ${KIND}$ - real(${KIND}$), intent(in) :: x - real(${KIND}$) :: f - end function integrand_${KIND}$ + #:for k1, t1 in REAL_KINDS_TYPES + pure function integrand_${k1}$(x) result(f) + import :: ${k1}$ + ${t1}$, intent(in) :: x + ${t1}$ :: f + end function integrand_${k1}$ #:endfor end interface @@ -119,18 +120,18 @@ module stdlib_experimental_quadrature ! PDT bug? interface quad #:for WFUN in WEIGHT_FUNS - #:for KIND in REAL_KINDS - module function quad_${WFUN}$_${KIND}$(f, a, b, weight, points, abstol, reltol, delta) result(integral) - procedure(integrand_${KIND}$) :: f - real(${KIND}$), intent(in) :: a - real(${KIND}$), intent(in) :: b - type(${WFUN}$_weight_t(${KIND}$)), intent(in) :: weight - real(${KIND}$), intent(in), dimension(:) :: points - real(${KIND}$), intent(in) :: abstol - real(${KIND}$), intent(in) :: reltol - real(${KIND}$), intent(out), optional :: delta - real(${KIND}$) :: integral - end function quad_${WFUN}$_${KIND}$ + #:for k1, t1 in REAL_KINDS_TYPES + module function quad_${WFUN}$_${k1}$(f, a, b, weight, points, abstol, reltol, delta) result(integral) + procedure(integrand_${k1}$) :: f + ${t1}$, intent(in) :: a + ${t1}$, intent(in) :: b + type(${WFUN}$_weight_t(${k1}$)), intent(in) :: weight + ${t1}$, intent(in), dimension(:) :: points + ${t1}$, intent(in) :: abstol + ${t1}$, intent(in) :: reltol + ${t1}$, intent(out), optional :: delta + ${t1}$ :: integral + end function quad_${WFUN}$_${k1}$ #:endfor #:endfor end interface quad From b7b273b7e7dd70ea2f3dd2a04cab906b759d4a29 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sun, 31 May 2020 17:51:54 -0600 Subject: [PATCH 03/15] dot_product->sum in trapz_weights example, looking toward complex --- src/stdlib_experimental_quadrature.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md index c62638f64..c46dde0dd 100644 --- a/src/stdlib_experimental_quadrature.md +++ b/src/stdlib_experimental_quadrature.md @@ -71,8 +71,8 @@ program demo_trapz_weights real :: x(5) = [0., 1., 2., 3., 4.] real :: y(5) = x**2 real :: w(5) - w = trapz_weight(x) - print *, dot_product(w, y) + w = trapz_weights(x) + print *, sum(w*y) ! 22.0 end program From d1db8a7a6dcbaaa70b863837604f78c7b124808d Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sun, 31 May 2020 17:54:11 -0600 Subject: [PATCH 04/15] change semantics of `even` in simps* spec --- src/stdlib_experimental_quadrature.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/stdlib_experimental_quadrature.md b/src/stdlib_experimental_quadrature.md index c46dde0dd..ed2689392 100644 --- a/src/stdlib_experimental_quadrature.md +++ b/src/stdlib_experimental_quadrature.md @@ -82,7 +82,7 @@ end program Returns the Simpson's rule integral of an array `y` representing discrete samples of a function. The integral is computed assuming either equidistant abscissas with spacing `dx` or arbitary abscissas `x`. -Simpson's rule is defined for odd-length arrays only. For even-length arrays, an optional argument `even` may be used to specify at which index to replace Simpson's rule with Simpson's 3/8 rule. The 3/8 rule will be used for the array section `y(even:even+4)` and the ordinary Simpon's rule will be used elsewhere. +Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of `even`. If `even` is negative (positive), the 3/8 rule is used at the beginning (end) of the array. If `even` is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and these two results were averaged. ### Syntax @@ -98,7 +98,7 @@ Simpson's rule is defined for odd-length arrays only. For even-length arrays, an `dx`: Shall be a scalar of type `real` having the same kind as `y`. -`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`. +`even`: (Optional) Shall be a default-kind `integer`. ### Return value @@ -106,7 +106,7 @@ The result is a scalar of type `real` having the same kind as `y`. If the size of `y` is zero or one, the result is zero. -If the size of `y` is two, the result is the same as if `trapz` had been called instead, regardless of the value of `even`. +If the size of `y` is two, the result is the same as if `trapz` had been called instead. ### Example @@ -116,7 +116,7 @@ TBD Given an array of abscissas `x`, computes the array of weights `w` such that if `y` represented function values tabulated at `x`, then `sum(w*y)` produces a Simpson's rule approximation to the integral. -Simpson's rule is defined for odd-length arrays only. For even-length arrays, an optional argument `even` may be used to specify at which index to replace Simpson's rule with Simpson's 3/8 rule. The 3/8 rule will be used for the array section `x(even:even+4)` and the ordinary Simpon's rule will be used elsewhere. +Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length arrays, Simpson's 3/8 rule is also utilized in a way that depends on the value of `even`. If `even` is negative (positive), the 3/8 rule is used at the beginning (end) of the array and the 1/3 rule used elsewhere. If `even` is zero or not present, the result is as if the 3/8 rule were first used at the beginning of the array, then at the end of the array, and then these two results were averaged. ### Syntax @@ -126,7 +126,7 @@ Simpson's rule is defined for odd-length arrays only. For even-length arrays, an `x`: Shall be a rank-one array of type `real`. -`even`: (Optional) Shall be a scalar integer of default kind. Its default value is `1`. +`even`: (Optional) Shall be a default-kind `integer`. ### Return value From cd4ef201bab672f96f8f8969d6eac7a1b97da36a Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sun, 31 May 2020 17:54:53 -0600 Subject: [PATCH 05/15] implement Simpson's rule --- src/CMakeLists.txt | 1 + src/stdlib_experimental_quadrature_simps.fypp | 277 ++++++++++ src/tests/quadrature/CMakeLists.txt | 2 +- src/tests/quadrature/test_simps.f90 | 505 ++++++++++++++++++ 4 files changed, 784 insertions(+), 1 deletion(-) create mode 100644 src/stdlib_experimental_quadrature_simps.fypp create mode 100644 src/tests/quadrature/test_simps.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 480e2737a..75362ccaf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,7 @@ set(fppFiles stdlib_experimental_stats_mean.fypp stdlib_experimental_quadrature.fypp stdlib_experimental_quadrature_trapz.fypp + stdlib_experimental_quadrature_simps.fypp ) diff --git a/src/stdlib_experimental_quadrature_simps.fypp b/src/stdlib_experimental_quadrature_simps.fypp new file mode 100644 index 000000000..858d79e1b --- /dev/null +++ b/src/stdlib_experimental_quadrature_simps.fypp @@ -0,0 +1,277 @@ +#:include "common.fypp" + +submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_simps + + implicit none + + ! internal use only + interface simps38 + #:for k1, t1 in REAL_KINDS_TYPES + module procedure simps38_dx_${k1}$ + module procedure simps38_x_${k1}$ + #:endfor + end interface simps38 + + ! internal use only + interface simps38_weights + #:for k1, t1 in REAL_KINDS_TYPES + module procedure simps38_weights_${k1}$ + #:endfor + end interface simps38_weights + +contains + +#:for k1, t1 in REAL_KINDS_TYPES + + pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, intent(in) :: dx + integer, intent(in), optional :: even + ${t1}$ :: integral + + integer :: n + + n = size(y) + + select case (n) + case (0:1) + integral = 0.0_${k1}$ + case (2) + integral = 0.5_${k1}$*dx*(y(1) + y(2)) + case (3) + integral = dx/3.0_${k1}$*(y(1) + 4*y(2) + y(3)) + case (4) + integral = simps38(y, dx) + ! case (5) not needed; handled by default + case (6) ! needs special handling because of averaged 3/8's rule case + if (present(even)) then + if (even < 0) then + ! 3/8 rule on left + integral = simps38(y(1:4), dx) + simps(y(4:6), dx) + return + else if (even > 0) then + ! 3/8 rule on right + integral = simps(y(1:3), dx) + simps38(y(3:6), dx) + return + else + ! fall through + end if + end if + ! either `even` not present or is zero + ! equivalent to averaging left and right + integral = dx/48.0_${k1}$ * (17*(y(1) + y(6)) + 59*(y(2) + y(5)) + 44*(y(3) + y(4))) + case default + if (mod(n, 2) == 1) then + integral = dx/3.0_${k1}$*(y(1) + 4*sum(y(2:n-1:2)) + 2*sum(y(3:n-2:2)) + y(n)) + else + if (present(even)) then + if (even < 0) then + ! 3/8th rule on left + integral = simps38(y(1:4), dx) + simps(y(4:n), dx) + return + else if (even > 0) then + ! 3/8 rule on right + integral = simps(y(1:n-3), dx) + simps38(y(n-3:n), dx) + return + else + ! fall through + end if + end if + ! either `even` not present or is zero + ! equivalent to averaging left and right + integral = dx/48.0_${k1}$ * (17*(y(1) + y(n)) + 59*(y(2) + y(n-1)) & + + 43*(y(3) + y(n-2)) + 49*(y(4) + y(n-3)) + 48*sum(y(5:n-4))) + end if + end select + end function simps_dx_${k1}$ + +#:endfor +#:for k1, t1 in REAL_KINDS_TYPES + + pure recursive module function simps_x_${k1}$(y, x, even) result(integral) + ${t1}$, dimension(:), intent(in) :: y + ${t1}$, dimension(size(y)), intent(in) :: x + integer, intent(in), optional :: even + ${t1}$ :: integral + + integer :: i + integer :: n + + ${t1}$ :: h1, h2, h3 + ${t1}$ :: a, b, c, d + + n = size(y) + + select case (n) + case (0:1) + integral = 0.0_${k1}$ + case (2) + integral = 0.5_${k1}$*(x(2) - x(1))*(y(1) + y(2)) + case (3) + h1 = x(2) - x(1) + h2 = x(3) - x(2) + a = (2*h1**2 + h1*h2 - h2**2)/(6*h1) + b = (h1+h2)**3/(6*h1*h2) + c = (2*h2**2 + h1*h2 - h1**2)/(6*h2) + integral = a*y(1) + b*y(2) + c*y(3) + case (4) + integral = simps38(y, x) + ! case (6) unneeded? + case default + if (mod(n, 2) == 1) then + integral = 0.0_${k1}$ + do i = 1, n-2, 2 + h1 = x(i+1) - x(i) + h2 = x(i+2) - x(i+1) + a = (2*h1**2 + h1*h2 - h2**2)/(6*h1) + b = (h1+h2)**3/(6*h1*h2) + c = (2*h2**2 + h1*h2 - h1**2)/(6*h2) + integral = integral + a*y(i) + b*y(i+1) + c*y(i+2) + end do + else + if (present(even)) then + if (even < 0) then + ! 3/8 rule on left + integral = simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n)) + return + else if (even > 0) then + ! 3/8 rule on right + integral = simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n)) + return + else + ! fall through + end if + end if + ! either `even` not present or is zero + integral = 0.5_${k1}$ * ( simps38(y(1:4), x(1:4)) + simps(y(4:n), x(4:n)) & + + simps(y(1:n-3), x(1:n-3)) + simps38(y(n-3:n), x(n-3:n)) ) + end if + end select + end function simps_x_${k1}$ + +#:endfor +#:for k1, t1 in REAL_KINDS_TYPES + + pure recursive module function simps_weights_${k1}$(x, even) result(w) + ${t1}$, dimension(:), intent(in) :: x + integer, intent(in), optional :: even + ${t1}$, dimension(size(x)) :: w + + integer :: i, n + ${t1}$ :: h1, h2, h3 + + n = size(x) + + select case (n) + case (0) + ! no action needed + case (1) + w(1) = 0.0_${k1}$ + case (2) + w = 0.5_${k1}$*(x(2) - x(1)) + case (3) + h1 = x(2) - x(1) + h2 = x(3) - x(2) + w(1) = (2*h1**2 + h1*h2 - h2**2)/(6*h1) + w(2) = (h1+h2)**3/(6*h1*h2) + w(3) = (2*h2**2 + h1*h2 - h1**2)/(6*h2) + case (4) + w = simps38_weights(x) + case default + if (mod(n, 2) == 1) then + w = 0.0_${k1}$ + do i = 1, n-2, 2 + h1 = x(i+1) - x(i) + h2 = x(i+2) - x(i+1) + w(i) = w(i) + (2*h1**2 + h1*h2 - h2**2)/(6*h1) + w(i+1) = w(i+1) + (h1+h2)**3/(6*h1*h2) + w(i+2) = w(i+2) + (2*h2**2 + h1*h2 - h1**2)/(6*h2) + end do + else + if (present(even)) then + if (even < 0) then + ! 3/8 rule on left + w = 0.0_${k1}$ + w(1:4) = simps38_weights(x(1:4)) + w(4:n) = w(4:n) + simps_weights(x(4:n)) ! position 4 needs both rules + return + else if (even > 0) then + ! 3/8 rule on right + w = 0.0_${k1}$ + w(1:n-3) = simps_weights(x(1:n-3)) + w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n)) ! position n-3 needs both rules + return + else + ! fall through + end if + end if + ! either `even` not present or is zero + w = 0.0_${k1}$ + ! 3/8 rule on left + w(1:4) = simps38_weights(x(1:4)) + w(4:n) = w(4:n) + simps_weights(x(4:n)) + ! 3/8 rule on right + w(1:n-3) = w(1:n-3) + simps_weights(x(1:n-3)) + w(n-3:n) = w(n-3:n) + simps38_weights(x(n-3:n)) + ! average + w = 0.5_${k1}$ * w + end if + end select + end function simps_weights_${k1}$ + +#:endfor +#:for k1, t1 in REAL_KINDS_TYPES + + pure function simps38_dx_${k1}$(y, dx) result (integral) + ${t1}$, dimension(4), intent(in) :: y + ${t1}$, intent(in) :: dx + ${t1}$ :: integral + + integral = 3.0_${k1}$*dx/8.0_${k1}$ * (y(1) + y(4) + 3*(y(2) + y(3))) + end function simps38_dx_${k1}$ + +#:endfor +#: for k1, t1 in REAL_KINDS_TYPES + + pure function simps38_x_${k1}$(y, x) result(integral) + ${t1}$, dimension(4), intent(in) :: y + ${t1}$, dimension(4), intent(in) :: x + ${t1}$ :: integral + + ${t1}$ :: h1, h2, h3 + ${t1}$ :: a, b, c, d + + h1 = x(2) - x(1) + h2 = x(3) - x(2) + h3 = x(4) - x(3) + + a = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2)) + b = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3)) + c = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2)) + d = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3)) + + integral = a*y(1) + b*y(2) + c*y(3) + d*y(4) + end function simps38_x_${k1}$ + +#:endfor +#:for k1, t1 in REAL_KINDS_TYPES + + pure function simps38_weights_${k1}$(x) result(w) + ${t1}$, intent(in) :: x(4) + ${t1}$ :: w(size(x)) + + ${t1}$ :: h1, h2, h3 + + h1 = x(2) - x(1) + h2 = x(3) - x(2) + h3 = x(4) - x(3) + + w(1) = (h1+h2+h3)*(3*h1**2 + 2*h1*h2 - 2*h1*h3 - h2**2 + h3**2)/(12*h1*(h1+h2)) + w(2) = (h1+h2-h3)*(h1+h2+h3)**3/(12*h1*h2*(h2+h3)) + w(3) = (h2+h3-h1)*(h1+h2+h3)**3/(12*h2*h3*(h1+h2)) + w(4) = (h1+h2+h3)*(3*h3**2 + 2*h2*h3 - 2*h1*h3 - h2**2 + h1**2)/(12*h3*(h2+h3)) + end function simps38_weights_${k1}$ + +#:endfor + +end submodule stdlib_experimental_quadrature_simps diff --git a/src/tests/quadrature/CMakeLists.txt b/src/tests/quadrature/CMakeLists.txt index 08f1ce4c6..d8069fb86 100644 --- a/src/tests/quadrature/CMakeLists.txt +++ b/src/tests/quadrature/CMakeLists.txt @@ -1,2 +1,2 @@ ADDTEST(trapz) - +ADDTEST(simps) diff --git a/src/tests/quadrature/test_simps.f90 b/src/tests/quadrature/test_simps.f90 new file mode 100644 index 000000000..800f676bb --- /dev/null +++ b/src/tests/quadrature/test_simps.f90 @@ -0,0 +1,505 @@ +program test_simps + use stdlib_experimental_kinds, only: sp, dp, qp + use stdlib_experimental_error, only: assert + use stdlib_experimental_quadrature, only: simps, simps_weights + + implicit none + + real(sp) :: tol_sp = epsilon(1.0_sp) + real(dp) :: tol_dp = epsilon(1.0_dp) + real(qp) :: tol_qp = epsilon(1.0_sp) + + call test_simps_sp + call test_simps_dp + call test_simps_qp + + call test_simps_weights_sp + call test_simps_weights_dp + call test_simps_weights_qp + + call test_simps_zero_sp + call test_simps_zero_dp + call test_simps_zero_qp + + call test_simps_even_sp + call test_simps_even_dp + call test_simps_even_qp + + call test_simps_weights_even_sp + call test_simps_weights_even_dp + call test_simps_weights_even_qp + + call test_simps_six_sp + call test_simps_six_dp + call test_simps_six_qp + +contains + + subroutine test_simps_sp + integer, parameter :: n = 13 + real(sp), dimension(n) :: y + real(sp), dimension(n) :: x + real(sp) :: val + real(sp) :: ans + integer :: i + + print *, "test_simps_sp" + + y = [(real(i-1, sp)**2, i = 1, n)] + + val = simps(y, 1.0_sp) + ans = 576.0_sp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_sp) + + val = simps(y, 0.5_sp) + ans = 288.0_sp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_sp) + + x = [(0.25_sp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 144.0_sp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_sp) + end subroutine test_simps_sp + + + subroutine test_simps_dp + integer, parameter :: n = 13 + real(dp), dimension(n) :: y + real(dp), dimension(n) :: x + real(dp) :: val + real(dp) :: ans + integer :: i + + print *, "test_simps_dp" + + y = [(real(i-1, dp)**2, i = 1, n)] + + val = simps(y, 1.0_dp) + ans = 576.0_dp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_dp) + + val = simps(y, 0.5_dp) + ans = 288.0_dp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_dp) + + x = [(0.25_dp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 144.0_dp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_dp) + end subroutine test_simps_dp + + + subroutine test_simps_qp + integer, parameter :: n = 13 + real(qp), dimension(n) :: y + real(qp), dimension(n) :: x + real(qp) :: val + real(qp) :: ans + integer :: i + + print *, "test_simps_qp" + + y = [(real(i-1, qp)**2, i = 1, n)] + + val = simps(y, 1.0_qp) + ans = 576.0_qp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_qp) + + val = simps(y, 0.5_qp) + ans = 288.0_qp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_qp) + + x = [(0.25_qp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 144.0_qp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_qp) + end subroutine test_simps_qp + + + subroutine test_simps_weights_sp + integer, parameter :: n = 17 + real(sp), dimension(n) :: y + real(sp), dimension(n) :: x + real(sp), dimension(n) :: w + integer :: i + real(sp) :: val + real(sp) :: ans + + print *, "test_simps_weights_sp" + + y = [(real(i-1, sp), i = 1, n)] + + x = y + w = simps_weights(x) + val = sum(w*y) + ans = simps(y, x) + print *, " ", val, ans + call assert(abs(val - ans) < tol_sp) + end subroutine test_simps_weights_sp + + + subroutine test_simps_weights_dp + integer, parameter :: n = 17 + real(dp), dimension(n) :: y + real(dp), dimension(n) :: x + real(dp), dimension(n) :: w + integer :: i + real(dp) :: val + real(dp) :: ans + + print *, "test_simps_weights_dp" + + y = [(real(i-1, dp), i = 1, n)] + + x = y + w = simps_weights(x) + val = sum(w*y) + ans = simps(y, x) + print *, " ", val, ans + call assert(abs(val - ans) < tol_dp) + end subroutine test_simps_weights_dp + + + subroutine test_simps_weights_qp + integer, parameter :: n = 17 + real(qp), dimension(n) :: y + real(qp), dimension(n) :: x + real(qp), dimension(n) :: w + integer :: i + real(qp) :: val + real(qp) :: ans + + print *, "test_simps_weights_qp" + + y = [(real(i-1, qp), i = 1, n)] + + x = y + w = simps_weights(x) + val = sum(w*y) + ans = simps(y, x) + print *, " ", val, ans + call assert(abs(val - ans) < tol_qp) + end subroutine test_simps_weights_qp + + + subroutine test_simps_zero_sp + real(sp), dimension(0) :: a + + print *, "test_simps_zero_sp" + + call assert(abs(simps(a, 1.0_sp)) < epsilon(0.0_sp)) + call assert(abs(simps([1.0_sp], 1.0_sp)) < epsilon(0.0_sp)) + call assert(abs(simps(a, a)) < epsilon(0.0_sp)) + call assert(abs(simps([1.0_sp], [1.0_sp])) < epsilon(0.0_sp)) + end subroutine test_simps_zero_sp + + + subroutine test_simps_zero_dp + real(dp), dimension(0) :: a + + print *, "test_simps_zero_dp" + + call assert(abs(simps(a, 1.0_dp)) < epsilon(0.0_dp)) + call assert(abs(simps([1.0_dp], 1.0_dp)) < epsilon(0.0_dp)) + call assert(abs(simps(a, a)) < epsilon(0.0_dp)) + call assert(abs(simps([1.0_dp], [1.0_dp])) < epsilon(0.0_dp)) + end subroutine test_simps_zero_dp + + + subroutine test_simps_zero_qp + real(qp), dimension(0) :: a + + print *, "test_simps_zero_qp" + + call assert(abs(simps(a, 1.0_qp)) < epsilon(0.0_qp)) + call assert(abs(simps([1.0_qp], 1.0_qp)) < epsilon(0.0_qp)) + call assert(abs(simps(a, a)) < epsilon(0.0_qp)) + call assert(abs(simps([1.0_qp], [1.0_qp])) < epsilon(0.0_qp)) + end subroutine test_simps_zero_qp + + + subroutine test_simps_even_sp + integer, parameter :: n = 11 + real(sp), dimension(n) :: y + real(sp), dimension(n) :: x + real(sp) :: val + real(sp) :: ans + integer :: i + integer :: even + + print *, "test_simps_even_sp" + + y = [(3.0_sp*real(i-1, sp)**2, i = 1, n)] + + do even = -1, 1 + print *, "even=", even + + val = simps(y, 1.0_sp) + ans = 1000.0_sp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_sp) + + val = simps(y, 0.5_sp) + ans = 500.0_sp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_sp) + + x = [(0.25_sp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 250.0_sp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_sp) + end do + end subroutine test_simps_even_sp + + + subroutine test_simps_even_dp + integer, parameter :: n = 11 + real(dp), dimension(n) :: y + real(dp), dimension(n) :: x + real(dp) :: val + real(dp) :: ans + integer :: i + + print *, "test_simps_even_dp" + + y = [(3.0_dp*real(i-1, dp)**2, i = 1, n)] + + val = simps(y, 1.0_dp) + ans = 1000.0_dp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_dp) + + val = simps(y, 0.5_dp) + ans = 500.0_dp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_dp) + + x = [(0.25_dp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 250.0_dp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_dp) + end subroutine test_simps_even_dp + + + subroutine test_simps_even_qp + integer, parameter :: n = 11 + real(qp), dimension(n) :: y + real(qp), dimension(n) :: x + real(qp) :: val + real(qp) :: ans + integer :: i + integer :: even + + print *, "test_simps_even_qp" + + y = [(3.0_qp*real(i-1, qp)**2, i = 1, n)] + + do even = -1, 1 + print *, " even=", even + + val = simps(y, 1.0_qp) + ans = 1000.0_qp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_qp) + + val = simps(y, 0.5_qp) + ans = 500.0_qp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_qp) + + x = [(0.25_qp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 250.0_qp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_qp) + end do + end subroutine test_simps_even_qp + + + subroutine test_simps_weights_even_sp + integer, parameter :: n = 16 + real(sp), dimension(n) :: y + real(sp), dimension(n) :: x + real(sp), dimension(n) :: w + integer :: i + real(sp) :: val + real(sp) :: ans + integer :: even + + print *, "test_simps_weights_even_sp" + + y = [(real(i-1, sp), i = 1, n)] + x = y + + do even = -1, 1 + w = simps_weights(x) + val = sum(w*y) + ans = simps(y, x) + print *, " even=", even, val, ans + call assert(abs(val - ans) < tol_sp) + end do + end subroutine test_simps_weights_even_sp + + + subroutine test_simps_weights_even_dp + integer, parameter :: n = 16 + real(dp), dimension(n) :: y + real(dp), dimension(n) :: x + real(dp), dimension(n) :: w + integer :: i + real(dp) :: val + real(dp) :: ans + integer :: even + + print *, "test_simps_weights_even_dp" + + y = [(real(i-1, dp), i = 1, n)] + x = y + + do even = -1, 1 + w = simps_weights(x) + val = sum(w*y) + ans = simps(y, x) + print *, " even=", even, val, ans + call assert(abs(val - ans) < tol_dp) + end do + end subroutine test_simps_weights_even_dp + + + subroutine test_simps_weights_even_qp + integer, parameter :: n = 16 + real(qp), dimension(n) :: y + real(qp), dimension(n) :: x + real(qp), dimension(n) :: w + integer :: i + real(qp) :: val + real(qp) :: ans + integer :: even + + print *, "test_simps_weights_even_qp" + + y = [(real(i-1, qp), i = 1, n)] + + x = y + + do even = -1, 1 + w = simps_weights(x) + val = sum(w*y) + ans = simps(y, x) + print *, " even=", even, val, ans + call assert(abs(val - ans) < tol_qp) + end do + end subroutine test_simps_weights_even_qp + + + subroutine test_simps_six_sp + integer, parameter :: n = 6 + real(sp), dimension(n) :: y + real(sp), dimension(n) :: x + real(sp) :: val + real(sp) :: ans + integer :: i + integer :: even + + print *, "test_simps_six_sp" + + y = [(3.0_sp*real(i-1, sp)**2, i = 1, n)] + + do even = -1, 1 + print *, "even=", even + + val = simps(y, 1.0_sp) + ans = 125.0_sp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_sp) + + val = simps(y, 0.5_sp) + ans = 62.5_sp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_sp) + + x = [(0.25_sp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 31.25_sp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_sp) + end do + end subroutine test_simps_six_sp + + + subroutine test_simps_six_dp + integer, parameter :: n = 6 + real(dp), dimension(n) :: y + real(dp), dimension(n) :: x + real(dp) :: val + real(dp) :: ans + integer :: i + + print *, "test_simps_six_dp" + + y = [(3.0_dp*real(i-1, dp)**2, i = 1, n)] + + val = simps(y, 1.0_dp) + ans = 125.0_dp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_dp) + + val = simps(y, 0.5_dp) + ans = 62.5_dp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_dp) + + x = [(0.25_dp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 31.25_dp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_dp) + end subroutine test_simps_six_dp + + + subroutine test_simps_six_qp + integer, parameter :: n = 6 + real(qp), dimension(n) :: y + real(qp), dimension(n) :: x + real(qp) :: val + real(qp) :: ans + integer :: i + integer :: even + + print *, "test_simps_six_qp" + + y = [(3.0_qp*real(i-1, qp)**2, i = 1, n)] + + do even = -1, 1 + print *, " even=", even + + val = simps(y, 1.0_qp) + ans = 125.0_qp + print *, " dx=1", val, ans + call assert(abs(val - ans) < tol_qp) + + val = simps(y, 0.5_qp) + ans = 62.5_qp + print *, " dx=0.5", val, ans + call assert(abs(val - ans) < tol_qp) + + x = [(0.25_qp*(i-1), i = 1, n)] + val = simps(y, x) + ans = 31.25_qp + print *, " x=0,0.25,0.5,...", val, ans + call assert(abs(val - ans) < tol_qp) + end do + end subroutine test_simps_six_qp + +end program From 178a8972cabecdfd81a57e09892842384eee4e24 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Tue, 2 Jun 2020 21:08:17 -0600 Subject: [PATCH 06/15] update assert -> check --- src/tests/quadrature/test_simps.f90 | 92 ++++++++++++++--------------- 1 file changed, 46 insertions(+), 46 deletions(-) diff --git a/src/tests/quadrature/test_simps.f90 b/src/tests/quadrature/test_simps.f90 index 800f676bb..8434e5419 100644 --- a/src/tests/quadrature/test_simps.f90 +++ b/src/tests/quadrature/test_simps.f90 @@ -1,6 +1,6 @@ program test_simps use stdlib_experimental_kinds, only: sp, dp, qp - use stdlib_experimental_error, only: assert + use stdlib_experimental_error, only: check use stdlib_experimental_quadrature, only: simps, simps_weights implicit none @@ -50,18 +50,18 @@ subroutine test_simps_sp val = simps(y, 1.0_sp) ans = 576.0_sp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) val = simps(y, 0.5_sp) ans = 288.0_sp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) x = [(0.25_sp*(i-1), i = 1, n)] val = simps(y, x) ans = 144.0_sp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) end subroutine test_simps_sp @@ -80,18 +80,18 @@ subroutine test_simps_dp val = simps(y, 1.0_dp) ans = 576.0_dp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) val = simps(y, 0.5_dp) ans = 288.0_dp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) x = [(0.25_dp*(i-1), i = 1, n)] val = simps(y, x) ans = 144.0_dp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) end subroutine test_simps_dp @@ -110,18 +110,18 @@ subroutine test_simps_qp val = simps(y, 1.0_qp) ans = 576.0_qp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) val = simps(y, 0.5_qp) ans = 288.0_qp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) x = [(0.25_qp*(i-1), i = 1, n)] val = simps(y, x) ans = 144.0_qp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) end subroutine test_simps_qp @@ -143,7 +143,7 @@ subroutine test_simps_weights_sp val = sum(w*y) ans = simps(y, x) print *, " ", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) end subroutine test_simps_weights_sp @@ -165,7 +165,7 @@ subroutine test_simps_weights_dp val = sum(w*y) ans = simps(y, x) print *, " ", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) end subroutine test_simps_weights_dp @@ -187,7 +187,7 @@ subroutine test_simps_weights_qp val = sum(w*y) ans = simps(y, x) print *, " ", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) end subroutine test_simps_weights_qp @@ -196,10 +196,10 @@ subroutine test_simps_zero_sp print *, "test_simps_zero_sp" - call assert(abs(simps(a, 1.0_sp)) < epsilon(0.0_sp)) - call assert(abs(simps([1.0_sp], 1.0_sp)) < epsilon(0.0_sp)) - call assert(abs(simps(a, a)) < epsilon(0.0_sp)) - call assert(abs(simps([1.0_sp], [1.0_sp])) < epsilon(0.0_sp)) + call check(abs(simps(a, 1.0_sp)) < epsilon(0.0_sp)) + call check(abs(simps([1.0_sp], 1.0_sp)) < epsilon(0.0_sp)) + call check(abs(simps(a, a)) < epsilon(0.0_sp)) + call check(abs(simps([1.0_sp], [1.0_sp])) < epsilon(0.0_sp)) end subroutine test_simps_zero_sp @@ -208,10 +208,10 @@ subroutine test_simps_zero_dp print *, "test_simps_zero_dp" - call assert(abs(simps(a, 1.0_dp)) < epsilon(0.0_dp)) - call assert(abs(simps([1.0_dp], 1.0_dp)) < epsilon(0.0_dp)) - call assert(abs(simps(a, a)) < epsilon(0.0_dp)) - call assert(abs(simps([1.0_dp], [1.0_dp])) < epsilon(0.0_dp)) + call check(abs(simps(a, 1.0_dp)) < epsilon(0.0_dp)) + call check(abs(simps([1.0_dp], 1.0_dp)) < epsilon(0.0_dp)) + call check(abs(simps(a, a)) < epsilon(0.0_dp)) + call check(abs(simps([1.0_dp], [1.0_dp])) < epsilon(0.0_dp)) end subroutine test_simps_zero_dp @@ -220,10 +220,10 @@ subroutine test_simps_zero_qp print *, "test_simps_zero_qp" - call assert(abs(simps(a, 1.0_qp)) < epsilon(0.0_qp)) - call assert(abs(simps([1.0_qp], 1.0_qp)) < epsilon(0.0_qp)) - call assert(abs(simps(a, a)) < epsilon(0.0_qp)) - call assert(abs(simps([1.0_qp], [1.0_qp])) < epsilon(0.0_qp)) + call check(abs(simps(a, 1.0_qp)) < epsilon(0.0_qp)) + call check(abs(simps([1.0_qp], 1.0_qp)) < epsilon(0.0_qp)) + call check(abs(simps(a, a)) < epsilon(0.0_qp)) + call check(abs(simps([1.0_qp], [1.0_qp])) < epsilon(0.0_qp)) end subroutine test_simps_zero_qp @@ -246,18 +246,18 @@ subroutine test_simps_even_sp val = simps(y, 1.0_sp) ans = 1000.0_sp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) val = simps(y, 0.5_sp) ans = 500.0_sp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) x = [(0.25_sp*(i-1), i = 1, n)] val = simps(y, x) ans = 250.0_sp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) end do end subroutine test_simps_even_sp @@ -277,18 +277,18 @@ subroutine test_simps_even_dp val = simps(y, 1.0_dp) ans = 1000.0_dp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) val = simps(y, 0.5_dp) ans = 500.0_dp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) x = [(0.25_dp*(i-1), i = 1, n)] val = simps(y, x) ans = 250.0_dp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) end subroutine test_simps_even_dp @@ -311,18 +311,18 @@ subroutine test_simps_even_qp val = simps(y, 1.0_qp) ans = 1000.0_qp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) val = simps(y, 0.5_qp) ans = 500.0_qp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) x = [(0.25_qp*(i-1), i = 1, n)] val = simps(y, x) ans = 250.0_qp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) end do end subroutine test_simps_even_qp @@ -347,7 +347,7 @@ subroutine test_simps_weights_even_sp val = sum(w*y) ans = simps(y, x) print *, " even=", even, val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) end do end subroutine test_simps_weights_even_sp @@ -372,7 +372,7 @@ subroutine test_simps_weights_even_dp val = sum(w*y) ans = simps(y, x) print *, " even=", even, val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) end do end subroutine test_simps_weights_even_dp @@ -398,7 +398,7 @@ subroutine test_simps_weights_even_qp val = sum(w*y) ans = simps(y, x) print *, " even=", even, val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) end do end subroutine test_simps_weights_even_qp @@ -422,18 +422,18 @@ subroutine test_simps_six_sp val = simps(y, 1.0_sp) ans = 125.0_sp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) val = simps(y, 0.5_sp) ans = 62.5_sp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) x = [(0.25_sp*(i-1), i = 1, n)] val = simps(y, x) ans = 31.25_sp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_sp) + call check(abs(val - ans) < tol_sp) end do end subroutine test_simps_six_sp @@ -453,18 +453,18 @@ subroutine test_simps_six_dp val = simps(y, 1.0_dp) ans = 125.0_dp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) val = simps(y, 0.5_dp) ans = 62.5_dp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) x = [(0.25_dp*(i-1), i = 1, n)] val = simps(y, x) ans = 31.25_dp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_dp) + call check(abs(val - ans) < tol_dp) end subroutine test_simps_six_dp @@ -487,18 +487,18 @@ subroutine test_simps_six_qp val = simps(y, 1.0_qp) ans = 125.0_qp print *, " dx=1", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) val = simps(y, 0.5_qp) ans = 62.5_qp print *, " dx=0.5", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) x = [(0.25_qp*(i-1), i = 1, n)] val = simps(y, x) ans = 31.25_qp print *, " x=0,0.25,0.5,...", val, ans - call assert(abs(val - ans) < tol_qp) + call check(abs(val - ans) < tol_qp) end do end subroutine test_simps_six_qp From 688c81a82e9cf2911e92b24a38714056796444f0 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Wed, 3 Jun 2020 08:20:24 -0600 Subject: [PATCH 07/15] relax tolerance in simps tests Co-authored-by: Jeremie Vandenplas --- src/tests/quadrature/test_simps.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tests/quadrature/test_simps.f90 b/src/tests/quadrature/test_simps.f90 index 8434e5419..b4cdc4287 100644 --- a/src/tests/quadrature/test_simps.f90 +++ b/src/tests/quadrature/test_simps.f90 @@ -5,9 +5,9 @@ program test_simps implicit none - real(sp) :: tol_sp = epsilon(1.0_sp) - real(dp) :: tol_dp = epsilon(1.0_dp) - real(qp) :: tol_qp = epsilon(1.0_sp) + real(sp), parameter :: tol_sp = 1000 * epsilon(1.0_sp) + real(dp), parameter :: tol_dp = 1000 * epsilon(1.0_dp) + real(qp), parameter :: tol_qp = 1000 *epsilon(1.0_sp) call test_simps_sp call test_simps_dp From 655693eec43c6b585fbfa5c672315663a988761f Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Wed, 3 Jun 2020 08:34:49 -0600 Subject: [PATCH 08/15] remove experimental weighted quad stuff --- src/stdlib_experimental_quadrature.fypp | 54 +------------------------ 1 file changed, 2 insertions(+), 52 deletions(-) diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp index 4dfc490a6..d36e19d8f 100644 --- a/src/stdlib_experimental_quadrature.fypp +++ b/src/stdlib_experimental_quadrature.fypp @@ -1,6 +1,4 @@ #:include "common.fypp" -#:set WEIGHT_FUNS = ["sin", "cos", "pole"] -#:set QUAD_OK = False module stdlib_experimental_quadrature !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description)) use stdlib_experimental_kinds, only: sp, dp, qp @@ -15,19 +13,10 @@ module stdlib_experimental_quadrature public :: simps public :: simps_weights - ! automatic integration of (weighted) functions - #:if QUAD_OK - public :: quad - public :: weight_t - #:for WFUN in WEIGHT_FUNS - public :: ${WFUN}$_weight_t - #:endfor - #:endif - interface trapz - !! Integrates sampled values using trapezoidal rule - !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description)) + !! Integrates sampled values using trapezoidal rule + !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description)) #:for k1, t1 in REAL_KINDS_TYPES pure module function trapz_dx_${k1}$(y, dx) result(integral) ${t1}$, dimension(:), intent(in) :: y @@ -102,43 +91,4 @@ module stdlib_experimental_quadrature #:endfor end interface -#:if QUAD_OK - ! Base class to avoid repeating kind parameter declaration. - type, abstract :: weight_t(kind) - integer, kind :: kind - end type weight_t - - type, extends(weight_t) :: sin_weight_t - real(kind) :: omega - end type sin_weight_t - - type, extends(weight_t) :: cos_weight_t - real(kind) :: omega - end type cos_weight_t - - type, extends(weight_t) :: pole_weight_t - real(kind) :: c - end type pole_weight_t - -! gfortran 9.2.0 chokes on ICE if I include this ("buffer overflow detected") -! Interestingly, though, the ICE happens while trying to build the trapz submodule -! PDT bug? - interface quad - #:for WFUN in WEIGHT_FUNS - #:for k1, t1 in REAL_KINDS_TYPES - module function quad_${WFUN}$_${k1}$(f, a, b, weight, points, abstol, reltol, delta) result(integral) - procedure(integrand_${k1}$) :: f - ${t1}$, intent(in) :: a - ${t1}$, intent(in) :: b - type(${WFUN}$_weight_t(${k1}$)), intent(in) :: weight - ${t1}$, intent(in), dimension(:) :: points - ${t1}$, intent(in) :: abstol - ${t1}$, intent(in) :: reltol - ${t1}$, intent(out), optional :: delta - ${t1}$ :: integral - end function quad_${WFUN}$_${k1}$ - #:endfor - #:endfor - end interface quad -#:endif end module stdlib_experimental_quadrature From 031ecf75407b910cd151eb624f91c0f54440c65e Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Wed, 3 Jun 2020 08:35:38 -0600 Subject: [PATCH 09/15] revise case (6) comment in simps_x_* --- src/stdlib_experimental_quadrature_simps.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_experimental_quadrature_simps.fypp b/src/stdlib_experimental_quadrature_simps.fypp index 858d79e1b..0bd8095f8 100644 --- a/src/stdlib_experimental_quadrature_simps.fypp +++ b/src/stdlib_experimental_quadrature_simps.fypp @@ -116,7 +116,7 @@ contains integral = a*y(1) + b*y(2) + c*y(3) case (4) integral = simps38(y, x) - ! case (6) unneeded? + ! case (6) unneeded; handled by default case default if (mod(n, 2) == 1) then integral = 0.0_${k1}$ From 321c9d0cbdd54c8ef0f2d26ec22c3aa8a9e862b8 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Thu, 4 Jun 2020 07:43:57 -0600 Subject: [PATCH 10/15] remove "to be implemented" from simps specs --- doc/specs/stdlib_experimental_quadrature.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_experimental_quadrature.md b/doc/specs/stdlib_experimental_quadrature.md index e01abb332..4b58d9c78 100644 --- a/doc/specs/stdlib_experimental_quadrature.md +++ b/doc/specs/stdlib_experimental_quadrature.md @@ -83,7 +83,7 @@ end program demo_trapz_weights ``` -## `simps` - integrate sampled values using Simpson's rule (to be implemented) +## `simps` - integrate sampled values using Simpson's rule ### Description @@ -119,7 +119,7 @@ If the size of `y` is two, the result is the same as if `trapz` had been called TBD -## `simps_weights` - Simpson's rule weights for given abscissas (to be implemented) +## `simps_weights` - Simpson's rule weights for given abscissas ### Description From a57b0cb2ed811c7d2a1ca1f2a7b4e8aabc88b356 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Thu, 4 Jun 2020 08:11:31 -0600 Subject: [PATCH 11/15] add examples for simps and simps_weights --- doc/specs/stdlib_experimental_quadrature.md | 26 +++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_experimental_quadrature.md b/doc/specs/stdlib_experimental_quadrature.md index 4b58d9c78..dd67195af 100644 --- a/doc/specs/stdlib_experimental_quadrature.md +++ b/doc/specs/stdlib_experimental_quadrature.md @@ -117,7 +117,18 @@ If the size of `y` is two, the result is the same as if `trapz` had been called ### Example -TBD +```fortran +program demo_simps + use stdlib_experimental_quadrature, only: simps + implicit none + real :: x(5) = [0., 1., 2., 3., 4.] + real :: y(5) = 3.*x**2 + print *, simps(y, x) +! 64.0 + print *, simps(y, 0.5) +! 32.0 +end program demo_simps +``` ## `simps_weights` - Simpson's rule weights for given abscissas @@ -147,4 +158,15 @@ If the size of `x` is two, then the result is the same as if `trapz_weights` had ### Example -TBD +```fortran +program demo_simps_weights + use stdlib_experimental_quadrature, only: simps_weights + implicit none + real :: x(5) = [0., 1., 2., 3., 4.] + real :: y(5) = 3.*x**2 + real :: w(5) + w = simps_weights(x) + print *, sum(w*y) +! 64.0 +end program demo_simps_weights +``` From 59abe88131482d6947064524910780995add84b5 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Thu, 4 Jun 2020 08:46:05 -0600 Subject: [PATCH 12/15] add doc links --- doc/specs/stdlib_experimental_quadrature.md | 6 +++--- src/stdlib_experimental_quadrature.fypp | 4 ++++ 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/doc/specs/stdlib_experimental_quadrature.md b/doc/specs/stdlib_experimental_quadrature.md index dd67195af..9bde409ce 100644 --- a/doc/specs/stdlib_experimental_quadrature.md +++ b/doc/specs/stdlib_experimental_quadrature.md @@ -93,9 +93,9 @@ Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length a ### Syntax -`result = simps(y, x [, even])` +`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, x [, even])` -`result = simps(y, dx [, even])` +`result = [[stdlib_experimental_quadrature(module):simps(interface)]](y, dx [, even])` ### Arguments @@ -140,7 +140,7 @@ Simpson's ordinary ("1/3") rule is used for odd-length arrays. For even-length a ### Syntax -`result = simps_weights(x [, even])` +`result = [[stdlib_experimental_quadrature(module):simps_weights(interface)]](x [, even])` ### Arguments diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp index d36e19d8f..a5151603e 100644 --- a/src/stdlib_experimental_quadrature.fypp +++ b/src/stdlib_experimental_quadrature.fypp @@ -47,6 +47,8 @@ module stdlib_experimental_quadrature interface simps + !! Integrates sampled values using Simpson's rule + !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3)) ! "recursive" is an implementation detail #:for k1, t1 in REAL_KINDS_TYPES pure recursive module function simps_dx_${k1}$(y, dx, even) result(integral) @@ -68,6 +70,8 @@ module stdlib_experimental_quadrature interface simps_weights + !! Integrates sampled values using trapezoidal rule weights for given abscissas + !! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_3)) #:for k1, t1 in REAL_KINDS_TYPES pure recursive module function simps_weights_${k1}$(x, even) result(w) ${t1}$, dimension(:), intent(in) :: x From f83a935e1939211e0cd2a68d121a6590f309d445 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sat, 6 Jun 2020 07:39:47 -0600 Subject: [PATCH 13/15] make x arrays assumed-shape w/ runtime size check --- src/stdlib_experimental_quadrature.fypp | 4 ++-- src/stdlib_experimental_quadrature_simps.fypp | 3 ++- src/stdlib_experimental_quadrature_trapz.fypp | 3 ++- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp index a5151603e..fdef16e64 100644 --- a/src/stdlib_experimental_quadrature.fypp +++ b/src/stdlib_experimental_quadrature.fypp @@ -27,7 +27,7 @@ module stdlib_experimental_quadrature #:for k1, t1 in REAL_KINDS_TYPES pure module function trapz_x_${k1}$(y, x) result(integral) ${t1}$, dimension(:), intent(in) :: y - ${t1}$, dimension(size(y)), intent(in) :: x + ${t1}$, dimension(:), intent(in) :: x ${t1}$ :: integral end function trapz_x_${k1}$ #:endfor @@ -61,7 +61,7 @@ module stdlib_experimental_quadrature #:for k1, t1 in REAL_KINDS_TYPES pure recursive module function simps_x_${k1}$(y, x, even) result(integral) ${t1}$, dimension(:), intent(in) :: y - ${t1}$, dimension(size(y)), intent(in) :: x + ${t1}$, dimension(:), intent(in) :: x integer, intent(in), optional :: even ${t1}$ :: integral end function simps_x_${k1}$ diff --git a/src/stdlib_experimental_quadrature_simps.fypp b/src/stdlib_experimental_quadrature_simps.fypp index 0bd8095f8..d8c9faeaf 100644 --- a/src/stdlib_experimental_quadrature_simps.fypp +++ b/src/stdlib_experimental_quadrature_simps.fypp @@ -90,7 +90,7 @@ contains pure recursive module function simps_x_${k1}$(y, x, even) result(integral) ${t1}$, dimension(:), intent(in) :: y - ${t1}$, dimension(size(y)), intent(in) :: x + ${t1}$, dimension(:), intent(in) :: x integer, intent(in), optional :: even ${t1}$ :: integral @@ -101,6 +101,7 @@ contains ${t1}$ :: a, b, c, d n = size(y) + if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size." select case (n) case (0:1) diff --git a/src/stdlib_experimental_quadrature_trapz.fypp b/src/stdlib_experimental_quadrature_trapz.fypp index 0ea66fe5e..31ff666ff 100644 --- a/src/stdlib_experimental_quadrature_trapz.fypp +++ b/src/stdlib_experimental_quadrature_trapz.fypp @@ -32,13 +32,14 @@ contains pure module function trapz_x_${KIND}$(y, x) result(integral) real(${KIND}$), dimension(:), intent(in) :: y - real(${KIND}$), dimension(size(y)), intent(in) :: x + real(${KIND}$), dimension(:), intent(in) :: x real(${KIND}$) :: integral integer :: i integer :: n n = size(y) + if (size(x) /= n) error stop "trapz: Arguments `x` and `y` must be the same size." select case (n) case (0:1) From 248d92ae3f72c48ff9e3e60e5dcf6fce8d56c96f Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sun, 7 Jun 2020 06:45:28 -0600 Subject: [PATCH 14/15] remove unused variables --- src/stdlib_experimental_quadrature_simps.fypp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stdlib_experimental_quadrature_simps.fypp b/src/stdlib_experimental_quadrature_simps.fypp index d8c9faeaf..748e5d81b 100644 --- a/src/stdlib_experimental_quadrature_simps.fypp +++ b/src/stdlib_experimental_quadrature_simps.fypp @@ -97,8 +97,8 @@ contains integer :: i integer :: n - ${t1}$ :: h1, h2, h3 - ${t1}$ :: a, b, c, d + ${t1}$ :: h1, h2 + ${t1}$ :: a, b, c n = size(y) if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size." @@ -159,7 +159,7 @@ contains ${t1}$, dimension(size(x)) :: w integer :: i, n - ${t1}$ :: h1, h2, h3 + ${t1}$ :: h1, h2 n = size(x) From 9b6c3d073eec9c667ab4d3049354f7def4a07371 Mon Sep 17 00:00:00 2001 From: Nathaniel Shaffer Date: Sun, 7 Jun 2020 17:17:54 -0600 Subject: [PATCH 15/15] {simps,trapz}_x_*: `error stop` -> `check`; drop `pure` --- src/stdlib_experimental_quadrature.fypp | 4 ++-- src/stdlib_experimental_quadrature_simps.fypp | 6 +++--- src/stdlib_experimental_quadrature_trapz.fypp | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/stdlib_experimental_quadrature.fypp b/src/stdlib_experimental_quadrature.fypp index fdef16e64..42270dfe8 100644 --- a/src/stdlib_experimental_quadrature.fypp +++ b/src/stdlib_experimental_quadrature.fypp @@ -25,7 +25,7 @@ module stdlib_experimental_quadrature end function trapz_dx_${k1}$ #:endfor #:for k1, t1 in REAL_KINDS_TYPES - pure module function trapz_x_${k1}$(y, x) result(integral) + module function trapz_x_${k1}$(y, x) result(integral) ${t1}$, dimension(:), intent(in) :: y ${t1}$, dimension(:), intent(in) :: x ${t1}$ :: integral @@ -59,7 +59,7 @@ module stdlib_experimental_quadrature end function simps_dx_${k1}$ #:endfor #:for k1, t1 in REAL_KINDS_TYPES - pure recursive module function simps_x_${k1}$(y, x, even) result(integral) + recursive module function simps_x_${k1}$(y, x, even) result(integral) ${t1}$, dimension(:), intent(in) :: y ${t1}$, dimension(:), intent(in) :: x integer, intent(in), optional :: even diff --git a/src/stdlib_experimental_quadrature_simps.fypp b/src/stdlib_experimental_quadrature_simps.fypp index 748e5d81b..8ff4c5437 100644 --- a/src/stdlib_experimental_quadrature_simps.fypp +++ b/src/stdlib_experimental_quadrature_simps.fypp @@ -1,7 +1,7 @@ #:include "common.fypp" submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_simps - + use stdlib_experimental_error, only: check implicit none ! internal use only @@ -88,7 +88,7 @@ contains #:endfor #:for k1, t1 in REAL_KINDS_TYPES - pure recursive module function simps_x_${k1}$(y, x, even) result(integral) + recursive module function simps_x_${k1}$(y, x, even) result(integral) ${t1}$, dimension(:), intent(in) :: y ${t1}$, dimension(:), intent(in) :: x integer, intent(in), optional :: even @@ -101,7 +101,7 @@ contains ${t1}$ :: a, b, c n = size(y) - if (size(x) /= n) error stop "simps: Arguments `x` and `y` must be the same size." + call check(size(x) == n, "simps: Arguments `x` and `y` must be the same size.") select case (n) case (0:1) diff --git a/src/stdlib_experimental_quadrature_trapz.fypp b/src/stdlib_experimental_quadrature_trapz.fypp index 31ff666ff..af445b1cb 100644 --- a/src/stdlib_experimental_quadrature_trapz.fypp +++ b/src/stdlib_experimental_quadrature_trapz.fypp @@ -1,7 +1,7 @@ #:include "common.fypp" submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz - + use stdlib_experimental_error, only: check implicit none contains @@ -29,8 +29,8 @@ contains #:endfor #:for KIND in REAL_KINDS - - pure module function trapz_x_${KIND}$(y, x) result(integral) + + module function trapz_x_${KIND}$(y, x) result(integral) real(${KIND}$), dimension(:), intent(in) :: y real(${KIND}$), dimension(:), intent(in) :: x real(${KIND}$) :: integral @@ -39,7 +39,7 @@ contains integer :: n n = size(y) - if (size(x) /= n) error stop "trapz: Arguments `x` and `y` must be the same size." + call check(size(x) == n, "trapz: Arguments `x` and `y` must be the same size.") select case (n) case (0:1)