Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Quadrature: spec tweaks & Simpson's rule #209

Merged
merged 20 commits into from
Jun 15, 2020
+879 −54
Merged

Quadrature: spec tweaks & Simpson's rule #209

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
f73beda
attempt at quad API
nshaffer Feb 10, 2020
eb3b22f
improve fypp usage
nshaffer May 31, 2020
b7b273b
dot_product->sum in trapz_weights example, looking toward complex
nshaffer May 31, 2020
d1db8a7
change semantics of `even` in simps* spec
nshaffer May 31, 2020
cd4ef20
implement Simpson's rule
nshaffer May 31, 2020
b2c6264
Merge branch 'dev-quadrature' of https://github.com/nshaffer/stdlib i…
nshaffer Jun 1, 2020
85f3df4
Merge branch 'master' of https://github.com/fortran-lang/stdlib
nshaffer Jun 3, 2020
178a897
update assert -> check
nshaffer Jun 3, 2020
688c81a
relax tolerance in simps tests
nshaffer Jun 3, 2020
655693e
remove experimental weighted quad stuff
nshaffer Jun 3, 2020
031ecf7
revise case (6) comment in simps_x_*
nshaffer Jun 3, 2020
15eee77
Merge branch 'dev-quadrature' of https://github.com/nshaffer/stdlib i…
nshaffer Jun 3, 2020
321c9d0
remove "to be implemented" from simps specs
nshaffer Jun 4, 2020
a57b0cb
add examples for simps and simps_weights
nshaffer Jun 4, 2020
59abe88
add doc links
nshaffer Jun 4, 2020
f83a935
make x arrays assumed-shape w/ runtime size check
nshaffer Jun 6, 2020
7961658
Merge branch 'master' of https://github.com/fortran-lang/stdlib
nshaffer Jun 6, 2020
248d92a
remove unused variables
nshaffer Jun 7, 2020
856e9be
Merge branch 'master' into dev-quadrature
nshaffer Jun 7, 2020
9b6c3d0
{simps,trapz}_x_*: `error stop` -> `check`; drop `pure`
nshaffer Jun 7, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions doc/specs/stdlib_experimental_quadrature.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,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 demo_trapz_weights

Expand All @@ -89,7 +89,7 @@ end program demo_trapz_weights

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

Expand All @@ -105,15 +105,15 @@ 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

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

Expand All @@ -125,7 +125,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

Expand All @@ -135,7 +135,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

Expand Down
1 change: 1 addition & 0 deletions src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ set(fppFiles
stdlib_experimental_stats_var.fypp
stdlib_experimental_quadrature.fypp
stdlib_experimental_quadrature_trapz.fypp
stdlib_experimental_quadrature_simps.fypp
)


Expand Down
138 changes: 101 additions & 37 deletions src/stdlib_experimental_quadrature.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#:include "common.fypp"

#:set WEIGHT_FUNS = ["sin", "cos", "pole"]
#:set QUAD_OK = False
nshaffer marked this conversation as resolved.
Show resolved Hide resolved
module stdlib_experimental_quadrature
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description))
use stdlib_experimental_kinds, only: sp, dp, qp
Expand All @@ -14,67 +15,130 @@ 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))
#: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}$
!! 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
${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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this a private auxiliary function where you know that the actual argument for x will be contiguous in practice, or is this a public user-callable function? (I haven't looked carefully at the whole component to know.) If it is the latter, I don't think you want to be declaring x as an explicit-shape array, dimension(size(y)). This implies F77-style argument passing, and if the actual argument isn't contiguous it will do an implicit copy-in of that argument to a contiguous temporary and pass it instead. Clearly x and y need to have the same size, but that style of declaration does nothing whatsoever to ensure that that is the case (in fact it hides the true length of the actual argument). This really requires a run-time check or assertion, with x declared as an assumed-shape array.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This routine is public by way of the generic trapz interface, so your comment is germane. I have no complaints about using assumed-shape declaration here, but I need to educate myself more on array argument passing to understand why it matters. It is a little surprising to learn that the explicit-shape declaration introduces a temporary for non-contiguous actual arguments, whereas the assumed-shape declaration does not.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is described is just the implementation detail of (a) specific compiler(s). Although I didn't have the time to check, I strongly doubt that any of such considerations are present in the standard.

Copy link
Member

@nncarlson nncarlson Jun 3, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's quite true that this isn't explicitly stated in the standard, however it is implied. The array ${t1}$, dimension(size(y)), intent(in) :: x is simply contiguous (6.5.4) whereas the corresponding actual argument need not be. The only way for that to work is for the array data to be reorganized into contiguous storage.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure that the only possible way to match a non-contiguous array with a contiguous one is by implicit copy-in? I don't know the answer.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the recommendation! I like the approach of making the in-argument assumed shape, and then later benchmarking against making it contiguous as well.

A related questions is if function results should also be assumed-shape. E.g.,

function trapz_weights_dp(x) result(w)
real(dp), dimension(:), intent(in) :: x
real(dp), dimension(size(x)) :: w ! explicit-shape result
! real(dp), dimension(:), allocatable :: w ! assumed-shape result
...
end function trapz_weights_dp

The discussion in #114 suggests that allocatable function results are undesirable. Do I understand right that an explicit-shape result is still expected to be the better option?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure there was any conclusion in #114 as to which performed better. The issue probably being whether the array was allocated on the stack or the heap. I haven't looked to see how you're using this particular function, but personally I've switched to generally avoiding array-valued functions altogether, as I feel this results in a lot of unnecessary allocations/deallocations of memory that passing the result as an argument avoids. I'm not really confident that compiler optimizers are able to eliminate these.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Strictly speaking, there is no such thing as an "assumed-shape" function result. It's just an allocatable.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@klausler, thanks for that excellent advice and explanation. @milancurcic we should put this recommendation somewhere in our tutorial section or FAQ, as this question will keep coming up often.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have just pushed a commit that makes the x argument of trapz_x_* and simps_x_* assumed-shape with a runtime check on its size. Thanks, @nncarlson, for raising the point, and @epagone and @klausler for the discussion.

${t1}$ :: integral
end function trapz_x_${k1}$
#:endfor
end interface trapz


interface trapz_weights
!! Integrates sampled values using trapezoidal rule weights for given abscissas
!! ([Specification](../page/specs/stdlib_experimental_quadrature.html#description_1))
#: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


! 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 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

#: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?
nshaffer marked this conversation as resolved.
Show resolved Hide resolved
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
Loading