Skip to content

Commit 8dc665d

Browse files
jacobwilliamsjacobwilliams
authored
Merge pull request #51 from jacobwilliams/50-work-arrays
Work arrays
2 parents c50fec1 + e5de5d1 commit 8dc665d

15 files changed

+2081
-1706
lines changed

‎.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,7 @@
2828
*.app
2929

3030
# Directories
31+
archive
3132
build
3233
doc
3334
lib

‎README.md

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -19,27 +19,27 @@ The core routines for the subroutine interface are:
1919
2020
!f(x)
2121
subroutine db1ink(x,nx,fcn,kx,iknot,tx,bcoef,iflag)
22-
subroutine db1val(xval,idx,tx,nx,kx,bcoef,f,iflag,inbvx,extrap)
22+
subroutine db1val(xval,idx,tx,nx,kx,bcoef,f,iflag,inbvx,w0,extrap)
2323
2424
!f(x,y)
2525
subroutine db2ink(x,nx,y,ny,fcn,kx,ky,iknot,tx,ty,bcoef,iflag)
26-
subroutine db2val(xval,yval,idx,idy,tx,ty,nx,ny,kx,ky,bcoef,f,iflag,inbvx,inbvy,iloy,extrap)
26+
subroutine db2val(xval,yval,idx,idy,tx,ty,nx,ny,kx,ky,bcoef,f,iflag,inbvx,inbvy,iloy,w1,w0,extrap)
2727
2828
!f(x,y,z)
2929
subroutine db3ink(x,nx,y,ny,z,nz,fcn,kx,ky,kz,iknot,tx,ty,tz,bcoef,iflag)
30-
subroutine db3val(xval,yval,zval,idx,idy,idz,tx,ty,tz,nx,ny,nz,kx,ky,kz,bcoef,f,iflag,inbvx,inbvy,inbvz,iloy,iloz,extrap)
30+
subroutine db3val(xval,yval,zval,idx,idy,idz,tx,ty,tz,nx,ny,nz,kx,ky,kz,bcoef,f,iflag,inbvx,inbvy,inbvz,iloy,iloz,w2,w1,w0,extrap)
3131
3232
!f(x,y,z,q)
3333
subroutine db4ink(x,nx,y,ny,z,nz,q,nq,fcn,kx,ky,kz,kq,iknot,tx,ty,tz,tq,bcoef,iflag)
34-
subroutine db4val(xval,yval,zval,qval,idx,idy,idz,idq,tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,bcoef,f,iflag,inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,extrap)
34+
subroutine db4val(xval,yval,zval,qval,idx,idy,idz,idq,tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,bcoef,f,iflag,inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,w3,w2,w1,w0,extrap)
3535
3636
!f(x,y,z,q,r)
3737
subroutine db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,bcoef,iflag)
38-
subroutine db5val(xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef,f,iflag,inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,extrap)
38+
subroutine db5val(xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef,f,iflag,inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,w4,w3,w2,w1,w0,extrap)
3939
4040
!f(x,y,z,q,r,s)
4141
subroutine db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,bcoef,iflag)
42-
subroutine db6val(xval,yval,zval,qval,rval,sval,idx,idy,idz,idq,idr,ids,tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,bcoef,f,iflag,inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,extrap)
42+
subroutine db6val(xval,yval,zval,qval,rval,sval,idx,idy,idz,idq,idr,ids,tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,bcoef,f,iflag,inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,w5,w4,w3,w2,w1,w0,extrap)
4343
```
4444

4545
The ```ink``` routines compute the interpolant coefficients, and the ```val``` routines evalute the interpolant at the specified value of each coordinate. The 2D and 3D routines are extensively refactored versions of the original routines from the [NIST Core Math Library](http://www.nist.gov/itl/math/mcsd-software.cfm). The others are new, and are simply extensions of the same algorithm into the other dimensions.
@@ -86,7 +86,7 @@ See the [examples](https://github.com/jacobwilliams/bspline-fortran/tree/master/
8686

8787
# Compiling
8888

89-
A simple bash script ```build.sh``` is provided for building bspline-fortran with gfortran using [FoBiS](https://github.com/szaghi/FoBiS). It also builds the API documentation using [FORD](https://github.com/cmacmackin/ford). The library can also be compiled with the Intel Fortran Compiler (and presumably any other Fortran compiler that supports modern standards).
89+
A simple bash script ```build.sh``` is provided for building bspline-fortran with gfortran using [FoBiS](https://github.com/szaghi/FoBiS). It also builds the API documentation using [FORD](https://github.com/Fortran-FOSS-Programmers/ford). The library can also be compiled with the Intel Fortran Compiler (and presumably any other Fortran compiler that supports modern standards).
9090

9191
A basic CMake configuration file is also included. For example, to build a static library:
9292

‎src/bspline_kinds_module.f90

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,8 @@ module bspline_kinds_module
1313

1414
private
1515

16-
integer,parameter,public :: wp = real64 !! Real precision
16+
integer,parameter,public :: wp = real64 !! Real working precision
17+
integer,parameter,public :: ip = int32 !! Integer working precision
1718

1819
!*****************************************************************************************
1920
end module bspline_kinds_module

‎src/bspline_oo_module.f90

Lines changed: 645 additions & 514 deletions
Large diffs are not rendered by default.

‎src/bspline_sub_module.f90

Lines changed: 1039 additions & 994 deletions
Large diffs are not rendered by default.

‎src/tests/bspline_extrap_test.f90

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -8,26 +8,27 @@
88
program bspline_extrap_test
99

1010
use bspline_module
11-
use bspline_kinds_module, only: wp
11+
use bspline_kinds_module, only: wp, ip
1212
use pyplot_module
1313

1414
implicit none
1515

16-
integer,parameter :: nx = 400 !! number of points in x
17-
integer,parameter :: nxv = 800 !! number of points to evaluate interpolant
16+
integer(ip),parameter :: nx = 400 !! number of points in x
17+
integer(ip),parameter :: nxv = 800 !! number of points to evaluate interpolant
1818

19-
integer,parameter :: kx = 4 !! order in x
20-
integer,parameter :: iknot = 0 !! automatically select the knots
19+
integer(ip),parameter :: kx = 4 !! order in x
20+
integer(ip),parameter :: iknot = 0 !! automatically select the knots
2121

2222
real(wp) :: x(nx)
2323
real(wp) :: xval(nxv),fval(nxv)
2424
real(wp) :: tx(nx+kx)
2525
real(wp) :: fcn_1d(nx)
2626
real(wp) :: val,tru,err,errmax
27-
integer :: i,j,idx,iflag,inbvx,iloy
27+
integer(ip) :: i,j,idx,iflag,inbvx,iloy
2828
logical :: extrap
2929
type(pyplot) :: plt
3030
integer :: istat !! pyplot-fortran status flag
31+
real(wp),dimension(3*kx) :: w1_1d !! work array
3132

3233
real(wp),parameter :: rad2deg = 180.0_wp / acos(-1.0_wp) !! deg. to radians conversion factor
3334

@@ -65,7 +66,7 @@ program bspline_extrap_test
6566

6667
errmax = 0.0_wp
6768
do i=1,nxv
68-
call db1val(xval(i),idx,tx,nx,kx,fcn_1d,val,iflag,inbvx,extrap=extrap)
69+
call db1val(xval(i),idx,tx,nx,kx,fcn_1d,val,iflag,inbvx,w1_1d,extrap=extrap)
6970
fval(i) = val ! save it for plot
7071
tru = f1(xval(i))
7172
err = abs(tru-val)

‎src/tests/bspline_test.f90

Lines changed: 58 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -5,25 +5,25 @@
55
program bspline_test
66

77
use bspline_module
8-
use bspline_kinds_module, only: wp
8+
use bspline_kinds_module, only: wp, ip
99

1010
implicit none
1111

12-
integer,parameter :: nx = 6 !! number of points in x
13-
integer,parameter :: ny = 6 !! number of points in y
14-
integer,parameter :: nz = 6 !! number of points in z
15-
integer,parameter :: nq = 6 !! number of points in q
16-
integer,parameter :: nr = 6 !! number of points in r
17-
integer,parameter :: ns = 6 !! number of points in s
12+
integer(ip),parameter :: nx = 6 !! number of points in x
13+
integer(ip),parameter :: ny = 6 !! number of points in y
14+
integer(ip),parameter :: nz = 6 !! number of points in z
15+
integer(ip),parameter :: nq = 6 !! number of points in q
16+
integer(ip),parameter :: nr = 6 !! number of points in r
17+
integer(ip),parameter :: ns = 6 !! number of points in s
1818

19-
integer,parameter :: kx = 4 !! order in x
20-
integer,parameter :: ky = 4 !! order in y
21-
integer,parameter :: kz = 4 !! order in z
22-
integer,parameter :: kq = 4 !! order in q
23-
integer,parameter :: kr = 4 !! order in r
24-
integer,parameter :: ks = 4 !! order in s
19+
integer(ip),parameter :: kx = 4 !! order in x
20+
integer(ip),parameter :: ky = 4 !! order in y
21+
integer(ip),parameter :: kz = 4 !! order in z
22+
integer(ip),parameter :: kq = 4 !! order in q
23+
integer(ip),parameter :: kr = 4 !! order in r
24+
integer(ip),parameter :: ks = 4 !! order in s
2525

26-
integer,parameter :: iknot = 0 !! automatically select the knots
26+
integer(ip),parameter :: iknot = 0 !! automatically select the knots
2727

2828
real(wp) :: x(nx),y(ny),z(nz),q(nq),r(nr),s(ns)
2929
real(wp) :: tx(nx+kx),ty(ny+ky),tz(nz+kz),tq(nq+kq),tr(nr+kr),ts(ns+ks)
@@ -34,13 +34,35 @@ program bspline_test
3434
real(wp) :: fcn_5d(nx,ny,nz,nq,nr)
3535
real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns)
3636

37+
real(wp),dimension(3*kx) :: w1_1d
38+
real(wp),dimension(ky) :: w1_2d
39+
real(wp),dimension(3*max(kx,ky)) :: w2_2d
40+
real(wp),dimension(ky,kz) :: w1_3d
41+
real(wp),dimension(kz) :: w2_3d
42+
real(wp),dimension(3*max(kx,ky,kz)) :: w3_3d
43+
real(wp),dimension(ky,kz,kq) :: w1_4d
44+
real(wp),dimension(kz,kq) :: w2_4d
45+
real(wp),dimension(kq) :: w3_4d
46+
real(wp),dimension(3*max(kx,ky,kz,kq)) :: w4_4d
47+
real(wp),dimension(ky,kz,kq,kr) :: w1_5d
48+
real(wp),dimension(kz,kq,kr) :: w2_5d
49+
real(wp),dimension(kq,kr) :: w3_5d
50+
real(wp),dimension(kr) :: w4_5d
51+
real(wp),dimension(3*max(kx,ky,kz,kq,kr)) :: w5_5d
52+
real(wp),dimension(ky,kz,kq,kr,ks) :: w1_6d
53+
real(wp),dimension(kz,kq,kr,ks) :: w2_6d
54+
real(wp),dimension(kq,kr,ks) :: w3_6d
55+
real(wp),dimension(kr,ks) :: w4_6d
56+
real(wp),dimension(ks) :: w5_6d
57+
real(wp),dimension(3*max(kx,ky,kz,kq,kr,ks)) :: w6_6d
58+
3759
real(wp) :: tol
3860
real(wp),dimension(6) :: val,tru,err,errmax
3961
logical :: fail
40-
integer :: i,j,k,l,m,n,idx,idy,idz,idq,idr,ids
41-
integer,dimension(6) :: iflag
42-
integer :: inbvx,inbvy,inbvz,inbvq,inbvr,inbvs
43-
integer :: iloy,iloz,iloq,ilor,ilos
62+
integer(ip) :: i,j,k,l,m,n,idx,idy,idz,idq,idr,ids
63+
integer(ip),dimension(6) :: iflag
64+
integer(ip) :: inbvx,inbvy,inbvz,inbvq,inbvr,inbvs
65+
integer(ip) :: iloy,iloz,iloq,ilor,ilos
4466

4567
fail = .false.
4668
tol = 1.0e-14_wp
@@ -52,22 +74,22 @@ program bspline_test
5274
ids = 0
5375

5476
do i=1,nx
55-
x(i) = dble(i-1)/dble(nx-1)
77+
x(i) = real(i-1,wp)/real(nx-1,wp)
5678
end do
5779
do j=1,ny
58-
y(j) = dble(j-1)/dble(ny-1)
80+
y(j) = real(j-1,wp)/real(ny-1,wp)
5981
end do
6082
do k=1,nz
61-
z(k) = dble(k-1)/dble(nz-1)
83+
z(k) = real(k-1,wp)/real(nz-1,wp)
6284
end do
6385
do l=1,nq
64-
q(l) = dble(l-1)/dble(nq-1)
86+
q(l) = real(l-1,wp)/real(nq-1,wp)
6587
end do
6688
do m=1,nr
67-
r(m) = dble(m-1)/dble(nr-1)
89+
r(m) = real(m-1,wp)/real(nr-1,wp)
6890
end do
6991
do n=1,ns
70-
s(n) = dble(n-1)/dble(ns-1)
92+
s(n) = real(n-1,wp)/real(ns-1,wp)
7193
end do
7294
do i=1,nx
7395
fcn_1d(i) = f1(x(i))
@@ -122,42 +144,48 @@ program bspline_test
122144
errmax = 0.0_wp
123145
do i=1,nx
124146
call db1val(x(i),idx,&
125-
tx,nx,kx,fcn_1d,val(1),iflag(1),inbvx)
147+
tx,nx,kx,fcn_1d,val(1),iflag(1),inbvx,&
148+
w1_1d)
126149
tru(1) = f1(x(i))
127150
err(1) = abs(tru(1)-val(1))
128151
errmax(1) = max(err(1),errmax(1))
129152
do j=1,ny
130153
call db2val(x(i),y(j),idx,idy,&
131154
tx,ty,nx,ny,kx,ky,fcn_2d,val(2),iflag(2),&
132-
inbvx,inbvy,iloy)
155+
inbvx,inbvy,iloy,&
156+
w1_2d,w2_2d)
133157
tru(2) = f2(x(i),y(j))
134158
err(2) = abs(tru(2)-val(2))
135159
errmax(2) = max(err(2),errmax(2))
136160
do k=1,nz
137161
call db3val(x(i),y(j),z(k),idx,idy,idz,&
138162
tx,ty,tz,nx,ny,nz,kx,ky,kz,fcn_3d,val(3),iflag(3),&
139-
inbvx,inbvy,inbvz,iloy,iloz)
163+
inbvx,inbvy,inbvz,iloy,iloz,&
164+
w1_3d,w2_3d,w3_3d)
140165
tru(3) = f3(x(i),y(j),z(k))
141166
err(3) = abs(tru(3)-val(3))
142167
errmax(3) = max(err(3),errmax(3))
143168
do l=1,nq
144169
call db4val(x(i),y(j),z(k),q(l),idx,idy,idz,idq,&
145170
tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,fcn_4d,val(4),iflag(4),&
146-
inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq)
171+
inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,&
172+
w1_4d,w2_4d,w3_4d,w4_4d)
147173
tru(4) = f4(x(i),y(j),z(k),q(l))
148174
err(4) = abs(tru(4)-val(4))
149175
errmax(4) = max(err(4),errmax(4))
150176
do m=1,nr
151177
call db5val(x(i),y(j),z(k),q(l),r(m),idx,idy,idz,idq,idr,&
152178
tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,fcn_5d,val(5),iflag(5),&
153-
inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor)
179+
inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,&
180+
w1_5d,w2_5d,w3_5d,w4_5d,w5_5d)
154181
tru(5) = f5(x(i),y(j),z(k),q(l),r(m))
155182
err(5) = abs(tru(5)-val(5))
156183
errmax(5) = max(err(5),errmax(5))
157184
do n=1,ns
158185
call db6val(x(i),y(j),z(k),q(l),r(m),s(n),idx,idy,idz,idq,idr,ids,&
159186
tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,fcn_6d,val(6),iflag(6),&
160-
inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos)
187+
inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,&
188+
w1_6d,w2_6d,w3_6d,w4_6d,w5_6d,w6_6d)
161189
tru(6) = f6(x(i),y(j),z(k),q(l),r(m),s(n))
162190
err(6) = abs(tru(6)-val(6))
163191
errmax(6) = max(err(6),errmax(6))

0 commit comments

Comments
 (0)