1
+ ! *****************************************************************************************
2
+ !
3
+ ! > Speed test for 1d-6d tensor product b-spline interpolation (object-oriented version).
4
+
5
+ program bspline_speed_test
6
+
7
+ use bspline_oo_module
8
+ use ,intrinsic :: iso_fortran_env, only: wp = > real64
9
+
10
+ implicit none
11
+
12
+ integer ,parameter :: nx = 6 ! number of points
13
+ integer ,parameter :: ny = 8
14
+ integer ,parameter :: nz = 10
15
+ integer ,parameter :: nq = 9
16
+ integer ,parameter :: nr = 7
17
+ integer ,parameter :: ns = 8
18
+
19
+ integer ,parameter :: kx = 2 ! order
20
+ integer ,parameter :: ky = 3
21
+ integer ,parameter :: kz = 4
22
+ integer ,parameter :: kq = 3
23
+ integer ,parameter :: kr = 2
24
+ integer ,parameter :: ks = 3
25
+
26
+ real (wp) :: x(nx),y(ny),z(nz),q(nq),r(nr),s(ns)
27
+ real (wp) :: fcn_1d(nx)
28
+ real (wp) :: fcn_2d(nx,ny)
29
+ real (wp) :: fcn_3d(nx,ny,nz)
30
+ real (wp) :: fcn_4d(nx,ny,nz,nq)
31
+ real (wp) :: fcn_5d(nx,ny,nz,nq,nr)
32
+ real (wp) :: fcn_6d(nx,ny,nz,nq,nr,ns)
33
+
34
+ type (bspline_1d) :: s1
35
+ type (bspline_2d) :: s2
36
+ type (bspline_3d) :: s3
37
+ type (bspline_4d) :: s4
38
+ type (bspline_5d) :: s5
39
+ type (bspline_6d) :: s6
40
+
41
+ real (wp) :: val,sumval
42
+ integer :: i,j,k,l,m,n,idx,idy,idz,idq,idr,ids,iflag,n_cases
43
+ real :: tstart, tend
44
+
45
+ idx = 0
46
+ idy = 0
47
+ idz = 0
48
+ idq = 0
49
+ idr = 0
50
+ ids = 0
51
+
52
+ x = [(dble (i), i= 1 ,nx)]
53
+ y = [(dble (i), i= 1 ,ny)]
54
+ z = [(dble (i), i= 1 ,nz)]
55
+ q = [(dble (i), i= 1 ,nq)]
56
+ r = [(dble (i), i= 1 ,nr)]
57
+ s = [(dble (i), i= 1 ,ns)]
58
+
59
+ do i= 1 ,nx
60
+ fcn_1d(i) = f1(x(i))
61
+ do j= 1 ,ny
62
+ fcn_2d(i,j) = f2(x(i),y(j))
63
+ do k= 1 ,nz
64
+ fcn_3d(i,j,k) = f3(x(i),y(j),z(k))
65
+ do l= 1 ,nq
66
+ fcn_4d(i,j,k,l) = f4(x(i),y(j),z(k),q(l))
67
+ do m= 1 ,nr
68
+ fcn_5d(i,j,k,l,m) = f5(x(i),y(j),z(k),q(l),r(m))
69
+ do n= 1 ,ns
70
+ fcn_6d(i,j,k,l,m,n) = f6(x(i),y(j),z(k),q(l),r(m),s(n))
71
+ end do
72
+ end do
73
+ end do
74
+ end do
75
+ end do
76
+ end do
77
+
78
+ ! initialize:
79
+
80
+ call s1% initialize(x,fcn_1d,kx,iflag)
81
+ call s2% initialize(x,y,fcn_2d,kx,ky,iflag)
82
+ call s3% initialize(x,y,z,fcn_3d,kx,ky,kz,iflag)
83
+ call s4% initialize(x,y,z,q,fcn_4d,kx,ky,kz,kq,iflag)
84
+ call s5% initialize(x,y,z,q,r,fcn_5d,kx,ky,kz,kq,kr,iflag)
85
+ call s6% initialize(x,y,z,q,r,s,fcn_6d,kx,ky,kz,kq,kr,ks,iflag)
86
+
87
+ ! compute max error at interpolation points
88
+ sumval = 0.0_wp
89
+ n_cases = nx* ny* nz* nq* nr* ns
90
+
91
+ call cpu_time(tstart)
92
+ do i= 1 ,nx
93
+ call s1% evaluate(x(i),idx,val,iflag)
94
+ sumval = sumval + val
95
+ do j= 1 ,ny
96
+ call s2% evaluate(x(i),y(j),idx,idy,val,iflag)
97
+ sumval = sumval + val
98
+ do k= 1 ,nz
99
+ call s3% evaluate(x(i),y(j),z(k),idx,idy,idz,val,iflag)
100
+ sumval = sumval + val
101
+ do l= 1 ,nq
102
+ call s4% evaluate(x(i),y(j),z(k),q(l),idx,idy,idz,idq,val,iflag)
103
+ sumval = sumval + val
104
+ do m= 1 ,nr
105
+ call s5% evaluate(x(i),y(j),z(k),q(l),r(m),idx,idy,idz,idq,idr,val,iflag)
106
+ sumval = sumval + val
107
+ do n= 1 ,ns
108
+ call s6% evaluate(x(i),y(j),z(k),q(l),r(m),s(n),idx,idy,idz,idq,idr,ids,val,iflag)
109
+ sumval = sumval + val
110
+ end do
111
+ end do
112
+ end do
113
+ end do
114
+ end do
115
+ end do
116
+ call cpu_time(tend)
117
+
118
+ write (* ,* ) ' result :' , sumval
119
+ write (* ,* ) ' number of cases:' , n_cases
120
+ write (* ,* ) ' cases/sec :' , n_cases/ (tend- tstart)
121
+
122
+ contains
123
+
124
+ real(wp) function f1 (x ) ! ! 1d test function
125
+ implicit none
126
+ real (wp) :: x
127
+ f1 = 0.5_wp * (x* exp (- x) + sin (x) )
128
+ end function f1
129
+
130
+ real(wp) function f2 (x ,y ) ! ! 2d test function
131
+ implicit none
132
+ real (wp) x,y,piov2
133
+ piov2 = 2.0_wp * atan (1.0_wp )
134
+ f2 = 0.5_wp * (y* exp (- x) + sin (piov2* y) )
135
+ end function f2
136
+
137
+ real(wp) function f3 (x ,y ,z ) ! ! 3d test function
138
+ implicit none
139
+ real (wp) x,y,z,piov2
140
+ piov2 = 2.0_wp * atan (1.0_wp )
141
+ f3 = 0.5_wp * ( y* exp (- x) + z* sin (piov2* y) )
142
+ end function f3
143
+
144
+ real(wp) function f4 (x ,y ,z ,q ) ! ! 4d test function
145
+ implicit none
146
+ real (wp) x,y,z,q,piov2
147
+ piov2 = 2.0_wp * atan (1.0_wp )
148
+ f4 = 0.5_wp * ( y* exp (- x) + z* sin (piov2* y) + q )
149
+ end function f4
150
+
151
+ real(wp) function f5 (x ,y ,z ,q ,r ) ! ! 5d test function
152
+ implicit none
153
+ real (wp) x,y,z,q,r,piov2
154
+ piov2 = 2.0_wp * atan (1.0_wp )
155
+ f5 = 0.5_wp * ( y* exp (- x) + z* sin (piov2* y) + q* r )
156
+ end function f5
157
+
158
+ real(wp) function f6 (x ,y ,z ,q ,r ,s ) ! ! 6d test function
159
+ implicit none
160
+ real (wp) x,y,z,q,r,s,piov2
161
+ piov2 = 2.0_wp * atan (1.0_wp )
162
+ f6 = 0.5_wp * ( y* exp (- x) + z* sin (piov2* y) + q* r + 2.0_wp * s )
163
+ end function f6
164
+
165
+ end program bspline_speed_test
0 commit comments