Skip to content

Commit 97ce9d3

Browse files
nshaffernshaffer
committed
initial commit -- trapz
1 parent ecfbfb0 commit 97ce9d3

File tree

7 files changed

+368
-0
lines changed

7 files changed

+368
-0
lines changed

‎src/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,8 @@ set(fppFiles
55
stdlib_experimental_io.fypp
66
stdlib_experimental_stats.fypp
77
stdlib_experimental_stats_mean.fypp
8+
stdlib_experimental_quadrature.fypp
9+
stdlib_experimental_quadrature_trapz.fypp
810
)
911

1012

‎src/stdlib_experimental_quadrature.fypp

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
#:include "common.fypp"
2+
3+
module stdlib_experimental_quadrature
4+
use stdlib_experimental_kinds, only: sp, dp, qp
5+
6+
implicit none
7+
8+
private
9+
public :: trapz, trapz_weights
10+
11+
interface trapz
12+
#:for KIND in REAL_KINDS
13+
pure module function trapz_dx_${KIND}$(y, dx) result(integral)
14+
real(${KIND}$), dimension(:), intent(in) :: y
15+
real(${KIND}$), intent(in) :: dx
16+
real(${KIND}$) :: integral
17+
end function trapz_dx_${KIND}$
18+
#:endfor
19+
#:for KIND in REAL_KINDS
20+
pure module function trapz_x_${KIND}$(y, x) result(integral)
21+
real(${KIND}$), dimension(:), intent(in) :: y
22+
real(${KIND}$), dimension(size(y)), intent(in) :: x
23+
real(${KIND}$) :: integral
24+
end function trapz_x_${KIND}$
25+
#:endfor
26+
end interface trapz
27+
28+
interface trapz_weights
29+
#:for KIND in REAL_KINDS
30+
pure module function trapz_weights_${KIND}$(x) result(w)
31+
real(${KIND}$), dimension(:), intent(in) :: x
32+
real(${KIND}$), dimension(size(x)) :: w
33+
end function trapz_weights_${KIND}$
34+
#:endfor
35+
end interface trapz_weights
36+
37+
end module stdlib_experimental_quadrature

‎src/stdlib_experimental_quadrature_trapz.fypp

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
#:include "common.fypp"
2+
3+
submodule (stdlib_experimental_quadrature) stdlib_experimental_quadrature_trapz
4+
use stdlib_experimental_kinds, only: sp, dp, qp
5+
6+
implicit none
7+
8+
contains
9+
10+
#:for KIND in REAL_KINDS
11+
12+
pure module function trapz_dx_${KIND}$(y, dx) result(integral)
13+
real(${KIND}$), dimension(:), intent(in) :: y
14+
real(${KIND}$), intent(in) :: dx
15+
real(${KIND}$) :: integral
16+
17+
integer :: n
18+
19+
n = size(y)
20+
21+
select case (n)
22+
case (0:1)
23+
integral = 0.0_${KIND}$
24+
case (2)
25+
integral = 0.5_${KIND}$*dx*(y(1) + y(2))
26+
case default
27+
integral = dx*(sum(y(2:n-1)) + 0.5_${KIND}$*(y(1) + y(n)))
28+
end select
29+
end function trapz_dx_${KIND}$
30+
31+
#:endfor
32+
#:for KIND in REAL_KINDS
33+
34+
pure module function trapz_x_${KIND}$(y, x) result(integral)
35+
real(${KIND}$), dimension(:), intent(in) :: y
36+
real(${KIND}$), dimension(size(y)), intent(in) :: x
37+
real(${KIND}$) :: integral
38+
39+
integer :: i
40+
integer :: n
41+
42+
n = size(y)
43+
44+
select case (n)
45+
case (0:1)
46+
integral = 0.0_${KIND}$
47+
case (2)
48+
integral = 0.5_${KIND}$*(x(2) - x(1))*(y(1) + y(2))
49+
case default
50+
integral = 0.0_${KIND}$
51+
do i = 2, n
52+
integral = integral + (x(i) - x(i-1))*(y(i) + y(i-1))
53+
end do
54+
integral = 0.5_${KIND}$*integral
55+
end select
56+
end function trapz_x_${KIND}$
57+
58+
#:endfor
59+
#:for KIND in REAL_KINDS
60+
61+
pure module function trapz_weights_${KIND}$(x) result(w)
62+
real(${KIND}$), dimension(:), intent(in) :: x
63+
real(${KIND}$), dimension(size(x)) :: w
64+
65+
integer :: i
66+
integer :: n
67+
68+
n = size(x)
69+
70+
select case (n)
71+
case (0)
72+
! no action needed
73+
case (1)
74+
w(1) = 0.0_${KIND}$
75+
case (2)
76+
w = 0.5_${KIND}$*(x(2) - x(1))
77+
case default
78+
w(1) = 0.5_${KIND}$*(x(2) - x(1))
79+
w(n) = 0.5_${KIND}$*(x(n) - x(n-1))
80+
do i = 2, size(x)-1
81+
w(i) = 0.5_${KIND}$*(x(i+1) - x(i-1))
82+
end do
83+
end select
84+
end function trapz_weights_${KIND}$
85+
86+
#:endfor
87+
end submodule stdlib_experimental_quadrature_trapz

‎src/tests/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ add_subdirectory(io)
1111
add_subdirectory(optval)
1212
add_subdirectory(stats)
1313
add_subdirectory(system)
14+
add_subdirectory(quadrature)
1415

1516
ADDTEST(always_skip)
1617
set_tests_properties(always_skip PROPERTIES SKIP_RETURN_CODE 77)

‎src/tests/quadrature/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
ADDTEST(trapz)
2+

‎src/tests/quadrature/Makefile.manual

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
PROGS_SRC = test_trapz.f90
2+
3+
include ../Makefile.manual.test.mk

‎src/tests/quadrature/test_trapz.f90

Lines changed: 236 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,236 @@
1+
program test_trapz
2+
use stdlib_experimental_kinds, only: sp, dp, qp
3+
use stdlib_experimental_error, only: assert
4+
use stdlib_experimental_quadrature, only: trapz, trapz_weights
5+
6+
implicit none
7+
8+
call test_trapz_sp
9+
call test_trapz_dp
10+
call test_trapz_qp
11+
12+
call test_trapz_weights_sp
13+
call test_trapz_weights_dp
14+
call test_trapz_weights_qp
15+
16+
call test_trapz_zero_sp
17+
call test_trapz_zero_dp
18+
call test_trapz_zero_qp
19+
20+
contains
21+
22+
subroutine test_trapz_sp
23+
integer, parameter :: n = 17
24+
real(sp), dimension(n) :: y
25+
real(sp), dimension(n) :: x
26+
real(sp) :: val
27+
real(sp) :: ans
28+
integer :: i
29+
30+
print *, "test_trapz_sp"
31+
32+
y = [(real(i-1, sp), i = 1, n)]
33+
34+
val = trapz(y, 1.0_sp)
35+
ans = 128.0_sp
36+
call assert(val == ans)
37+
38+
val = trapz(y, 0.5_sp)
39+
ans = 64.0_sp
40+
call assert(val == ans)
41+
42+
x = [((i-1)*4.0_sp/real(n-1, sp), i = 1, n)]
43+
val = trapz(y, x)
44+
ans = 32.0_sp
45+
call assert(val == ans)
46+
47+
x = y**2
48+
val = trapz(y, x)
49+
ans = 2728.0_sp
50+
call assert(val == ans)
51+
end subroutine test_trapz_sp
52+
53+
subroutine test_trapz_dp
54+
integer, parameter :: n = 17
55+
real(dp), dimension(n) :: y
56+
real(dp), dimension(n) :: x
57+
real(dp) :: val
58+
real(dp) :: ans
59+
integer :: i
60+
61+
print *, "test_trapz_dp"
62+
63+
y = [(real(i-1, dp), i = 1, n)]
64+
65+
val = trapz(y, 1.0_dp)
66+
ans = 128.0_dp
67+
call assert(val == ans)
68+
69+
val = trapz(y, 0.5_dp)
70+
ans = 64.0_dp
71+
call assert(val == ans)
72+
73+
x = [((i-1)*4.0_dp/real(n-1, dp), i = 1, n)]
74+
val = trapz(y, x)
75+
ans = 32.0_dp
76+
call assert(val == ans)
77+
78+
x = y**2
79+
val = trapz(y, x)
80+
ans = 2728.0_sp
81+
call assert(val == ans)
82+
end subroutine test_trapz_dp
83+
84+
85+
subroutine test_trapz_qp
86+
integer, parameter :: n = 17
87+
real(qp), dimension(n) :: y
88+
real(qp), dimension(n) :: x
89+
real(qp) :: val
90+
real(qp) :: ans
91+
integer :: i
92+
93+
print *, "test_trapz_qp"
94+
95+
y = [(real(i-1, qp), i = 1, n)]
96+
97+
val = trapz(y, 1.0_qp)
98+
ans = 128.0_qp
99+
call assert(val == ans)
100+
101+
val = trapz(y, 0.5_qp)
102+
ans = 64.0_qp
103+
call assert(val == ans)
104+
105+
x = [((i-1)*4.0_qp/real(n-1, qp), i = 1, n)]
106+
val = trapz(y, x)
107+
ans = 32.0_qp
108+
call assert(val == ans)
109+
110+
x = y**2
111+
val = trapz(y, x)
112+
ans = 2728.0_qp
113+
call assert(val == ans)
114+
end subroutine test_trapz_qp
115+
116+
117+
subroutine test_trapz_weights_sp
118+
integer, parameter :: n = 17
119+
real(sp), dimension(n) :: y
120+
real(sp), dimension(n) :: x
121+
real(sp), dimension(n) :: w
122+
integer :: i
123+
real(sp) :: val
124+
real(sp) :: ans
125+
126+
print *, "test_trapz_weights_sp"
127+
128+
y = [(real(i-1, sp), i = 1, n)]
129+
130+
x = y
131+
w = trapz_weights(x)
132+
val = dot_product(w, y)
133+
ans = trapz(y, x)
134+
call assert(val == ans)
135+
136+
x = y**2
137+
w = trapz_weights(x)
138+
val = dot_product(w, y)
139+
ans = trapz(y, x)
140+
call assert(val == ans)
141+
142+
end subroutine test_trapz_weights_sp
143+
144+
145+
subroutine test_trapz_weights_dp
146+
integer, parameter :: n = 17
147+
real(dp), dimension(n) :: y
148+
real(dp), dimension(n) :: x
149+
real(dp), dimension(n) :: w
150+
integer :: i
151+
real(dp) :: val
152+
real(dp) :: ans
153+
154+
print *, "test_trapz_weights_dp"
155+
156+
y = [(real(i-1, dp), i = 1, n)]
157+
158+
x = y
159+
w = trapz_weights(x)
160+
val = dot_product(w, y)
161+
ans = trapz(y, x)
162+
call assert(val == ans)
163+
164+
x = y**2
165+
w = trapz_weights(x)
166+
val = dot_product(w, y)
167+
ans = trapz(y, x)
168+
call assert(val == ans)
169+
170+
end subroutine test_trapz_weights_dp
171+
172+
173+
subroutine test_trapz_weights_qp
174+
integer, parameter :: n = 17
175+
real(qp), dimension(n) :: y
176+
real(qp), dimension(n) :: x
177+
real(qp), dimension(n) :: w
178+
integer :: i
179+
real(qp) :: val
180+
real(qp) :: ans
181+
182+
print *, "test_trapz_weights_qp"
183+
184+
y = [(real(i-1, qp), i = 1, n)]
185+
186+
x = y
187+
w = trapz_weights(x)
188+
val = dot_product(w, y)
189+
ans = trapz(y, x)
190+
call assert(val == ans)
191+
192+
x = y**2
193+
w = trapz_weights(x)
194+
val = dot_product(w, y)
195+
ans = trapz(y, x)
196+
call assert(val == ans)
197+
198+
end subroutine test_trapz_weights_qp
199+
200+
201+
subroutine test_trapz_zero_sp
202+
real(sp), dimension(0) :: a
203+
204+
print *, "test_trapz_zero_sp"
205+
206+
call assert(trapz(a, 1.0_sp) == 0.0_sp)
207+
call assert(trapz([1.0_sp], 1.0_sp) == 0.0_sp)
208+
call assert(trapz(a, a) == 0.0_sp)
209+
call assert(trapz([1.0_sp], [1.0_sp]) == 0.0_sp)
210+
end subroutine test_trapz_zero_sp
211+
212+
213+
subroutine test_trapz_zero_dp
214+
real(dp), dimension(0) :: a
215+
216+
print *, "test_trapz_zero_dp"
217+
218+
call assert(trapz(a, 1.0_dp) == 0.0_dp)
219+
call assert(trapz([1.0_dp], 1.0_dp) == 0.0_dp)
220+
call assert(trapz(a, a) == 0.0_dp)
221+
call assert(trapz([1.0_dp], [1.0_dp]) == 0.0_dp)
222+
end subroutine test_trapz_zero_dp
223+
224+
225+
subroutine test_trapz_zero_qp
226+
real(qp), dimension(0) :: a
227+
228+
print *, "test_trapz_zero_qp"
229+
230+
call assert(trapz(a, 1.0_qp) == 0.0_qp)
231+
call assert(trapz([1.0_qp], 1.0_qp) == 0.0_qp)
232+
call assert(trapz(a, a) == 0.0_qp)
233+
call assert(trapz([1.0_qp], [1.0_qp]) == 0.0_qp)
234+
end subroutine test_trapz_zero_qp
235+
236+
end program test_trapz

0 commit comments

Comments
 (0)