test updates to avoid aliased arguments to ink routines

jacobwilliams · jacobwilliams · commit 2f617bc9e0e6 · 2023-10-09T21:36:28.000-05:00
Fixes #89 also updated the linear test.
diff --git a/test/bspline_extrap_test.f90 b/test/bspline_extrap_test.f90
@@ -22,7 +22,7 @@ program bspline_extrap_test
     real(wp) :: x(nx)
     real(wp) :: xval(nxv),fval(nxv)
     real(wp) :: tx(nx+kx)
-    real(wp) :: fcn_1d(nx)
+    real(wp) :: fcn_1d(nx), bcoef(nx)
     real(wp) :: val,tru,err,errmax
     integer(ip) :: i,j,idx,iflag,inbvx,iloy
     logical :: extrap
@@ -47,7 +47,7 @@ program bspline_extrap_test
     iloy  = 1
 
     ! initialize
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag)
+    call db1ink(x,nx,fcn_1d,kx,iknot,tx,bcoef,iflag)
 
     if (iflag/=0) then
         write(*,*) 'Error initializing 1D spline: '//get_status_message(iflag)
@@ -67,7 +67,7 @@ program bspline_extrap_test
 
         errmax = 0.0_wp
         do i=1,nxv
-            call db1val(xval(i),idx,tx,nx,kx,fcn_1d,val,iflag,inbvx,w1_1d,extrap=extrap)
+            call db1val(xval(i),idx,tx,nx,kx,bcoef,val,iflag,inbvx,w1_1d,extrap=extrap)
             fval(i) = val  ! save it for plot
             tru    = f1(xval(i))
             err    = abs(tru-val)
diff --git a/test/bspline_linear_test.f90 b/test/bspline_linear_test.f90
@@ -12,29 +12,20 @@ program bspline_linear_test
     implicit none
 
     integer(ip),parameter :: nx    = 10    !! number of points in x
-    integer(ip),parameter :: nxv   = 111 !100   !! number of points to evaluate interpolant
-
-    integer(ip),parameter :: kx    = 2    !! order in x [linear]
-    !integer(ip),parameter :: kx    = bspline_order_cubic    !! order in x [cubic]
-    !integer(ip),parameter :: kx    = bspline_order_quartic
+    integer(ip),parameter :: nxv   = 111   !! number of points to evaluate interpolant
     integer(ip),parameter :: iknot = 0    !! automatically select the knots
+    integer,parameter :: idx = 0 !! compute value for the spline interpolation
 
     real(wp) :: x(nx)
-    real(wp) :: xval(nxv),fval(nxv),fval_linear(nxv)
-    real(wp) :: tx(nx+kx)
+    real(wp) :: xval(nxv),fval(nxv),fval_linear(nxv),fval3(nxv),fval4(nxv)
     real(wp) :: fcn_1d(nx)
-    real(wp) :: val,tru,err,errmax
-    integer(ip) :: i,idx,iflag,inbvx,iloy
-    type(pyplot) :: plt
+    real(wp) :: val,err,errmax
+    integer(ip) :: i,iflag,inbvx,iloy
     integer :: istat  !! pyplot-fortran status flag
-    real(wp),dimension(3*kx) :: w1_1d !! work array
+    type(pyplot) :: plt
+    type(bspline_1d) :: b2, b3, b4
     type(linear_interp_1d) :: s1
 
-    idx = 0
-
-    x = huge(1.0_wp)
-    xval = huge(1.0_wp)
-
     do i=1,nx
         x(i) = real(i-1,wp) * 10.0_wp
         fcn_1d(i) = f1(x(i))
@@ -47,59 +38,69 @@ program bspline_linear_test
     inbvx = 1
     iloy  = 1
 
+    ! the low-level routines are given the wrong result? check this...
+
     ! initialize
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag)
-    if (iflag/=0) then
-        write(*,*) 'Error initializing 1D spline: '//get_status_message(iflag)
-    end if
+    call b2%initialize(x,fcn_1d,2,iflag,extrap=.true.) ! linear
+    if (iflag/=0) error stop 'Error initializing 1D linear spline: '//get_status_message(iflag)
+    call b3%initialize(x,fcn_1d,3,iflag,extrap=.true.) ! quadratic
+    if (iflag/=0) error stop 'Error initializing 1D quadratic spline: '//get_status_message(iflag)
+    call b4%initialize(x,fcn_1d,4,iflag,extrap=.true.) ! cubic
+    if (iflag/=0) error stop 'Error initializing 1D cubic spline: '//get_status_message(iflag)
     call s1%initialize(x,fcn_1d,iflag)
-    if (iflag/=0) then
-        write(*,*) 'Error initializing 1D linear interpolator: ', iflag
-    end if
+    if (iflag/=0) error stop 'Error initializing 1D linear interpolator'
 
     !initialize the plot:
     call plt%initialize(grid=.true.,xlabel='x (deg)',ylabel='f(x)',&
-                        title='Linear Test',legend=.true.)
+                        title='Linear Test',legend=.true.,figsize=[10,5])
     call plt%add_plot(x,fcn_1d,&
                         label='Function $f(x) = \sin(x)$',&
                         linestyle='ko',markersize=5,linewidth=2,istat=istat)
 
     errmax = 0.0_wp
     do i=1,nxv
-        call db1val(xval(i),idx,tx,nx,kx,fcn_1d,val,iflag,inbvx,w1_1d,extrap=.true.)
+        call b2%evaluate(xval(i),idx,val,iflag)
         fval(i) = val  ! save it for plot
-        if (iflag/=0) error stop 'error calling db1val: '//get_status_message(iflag)
-
-        tru    = f1(xval(i))
-        err    = abs(tru-val)
-        errmax = max(err,errmax)
-        !write(*,*) xval(i), val, tru, err, iflag
+        if (iflag/=0) error stop 'error evaluating linear spline: '//get_status_message(iflag)
 
         ! also compute linear interpolation:
         call s1%evaluate(xval(i),val)
         fval_linear(i) = val ! linear vector for plot
-        write(*,*) "error : ", xval(i), fval(i) - fval_linear(i)
+
+        err    = abs(fval(i) - fval_linear(i))
+        errmax = max(err,errmax)
+
+        ! also others:
+        call b3%evaluate(xval(i),idx,fval3(i),iflag)
+        if (iflag/=0) error stop 'error evaluating quadratic spline: '//get_status_message(iflag)
+        call b4%evaluate(xval(i),idx,fval4(i),iflag)
+        if (iflag/=0) error stop 'error evaluating cubic spline: '//get_status_message(iflag)
 
     end do
 
-    ! check max error against tolerance
     write(*,*) ''
-    write(*,*) '1D: max error:', errmax
+    write(*,*) 'Max difference (spline - linear):', errmax
     write(*,*) ''
 
-    call plt%add_plot(xval,fval,&
-            label='k=2 Spline Interpolated',&
-            linestyle='g.-',linewidth=1,istat=istat)
     call plt%add_plot(xval,fval_linear,&
             label='Linear Interpolated',&
-            linestyle='r.-',linewidth=1,istat=istat)
+            linestyle='r.',linewidth=1,istat=istat)
+    call plt%add_plot(xval,fval,&
+            label='Linear ($k=2$) Spline Interpolated',&
+            linestyle='r-',linewidth=1,istat=istat)
+    call plt%add_plot(xval,fval3,&
+            label='Quadratic ($k=3$) Spline Interpolated',&
+            linestyle='k.-',linewidth=1,istat=istat)
+    call plt%add_plot(xval,fval4,&
+            label='Cubic ($k=4$) Spline Interpolated',&
+            linestyle='c.-',linewidth=1,istat=istat)
     call plt%savefig('bspline_linear_test.png',istat=istat)
 
     contains
 
         real(wp) function f1(x) !! 1d test function
         implicit none
-        real(wp) :: x
+        real(wp),intent(in) :: x
         real(wp),parameter :: a = acos(-1.0_wp)/18.0_wp
         f1 = sin(a*x)
         end function f1
diff --git a/test/bspline_test.f90 b/test/bspline_test.f90
@@ -27,12 +27,12 @@ program bspline_test
 
     real(wp) :: x(nx),y(ny),z(nz),q(nq),r(nr),s(ns)
     real(wp) :: tx(nx+kx),ty(ny+ky),tz(nz+kz),tq(nq+kq),tr(nr+kr),ts(ns+ks)
-    real(wp) :: fcn_1d(nx)
-    real(wp) :: fcn_2d(nx,ny)
-    real(wp) :: fcn_3d(nx,ny,nz)
-    real(wp) :: fcn_4d(nx,ny,nz,nq)
-    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)
-    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns)
+    real(wp) :: fcn_1d(nx)               , bcoef_1d(nx)
+    real(wp) :: fcn_2d(nx,ny)            , bcoef_2d(nx,ny)
+    real(wp) :: fcn_3d(nx,ny,nz)         , bcoef_3d(nx,ny,nz)
+    real(wp) :: fcn_4d(nx,ny,nz,nq)      , bcoef_4d(nx,ny,nz,nq)
+    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)   , bcoef_5d(nx,ny,nz,nq,nr)
+    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns), bcoef_6d(nx,ny,nz,nq,nr,ns)
 
     real(wp),dimension(3*kx)                     :: w1_1d
     real(wp),dimension(ky)                       :: w1_2d
@@ -229,12 +229,12 @@ program bspline_test
     ilos  = 1
 
     ! initialize
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag(1))
-    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,fcn_2d,iflag(2))
-    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,fcn_3d,iflag(3))
-    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,fcn_4d,iflag(4))
-    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,fcn_5d,iflag(5))
-    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,fcn_6d,iflag(6))
+    call db1ink(x,nx,fcn_1d,kx,iknot,tx,bcoef_1d,iflag(1))
+    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,bcoef_2d,iflag(2))
+    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,bcoef_3d,iflag(3))
+    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,bcoef_4d,iflag(4))
+    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,bcoef_5d,iflag(5))
+    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,bcoef_6d,iflag(6))
 
     if (any(iflag/=0)) then
         do i=1,6
@@ -249,46 +249,46 @@ program bspline_test
      errmax = 0.0_wp
      do i=1,nx
                         call db1val(x(i),idx,&
-                                            tx,nx,kx,fcn_1d,val(1),iflag(1),inbvx,&
+                                            tx,nx,kx,bcoef_1d,val(1),iflag(1),inbvx,&
                                             w1_1d)
                         tru(1)    = f1(x(i))
                         err(1)    = abs(tru(1)-val(1))
                         errmax(1) = max(err(1),errmax(1))
         do j=1,ny
                         call db2val(x(i),y(j),idx,idy,&
-                                            tx,ty,nx,ny,kx,ky,fcn_2d,val(2),iflag(2),&
+                                            tx,ty,nx,ny,kx,ky,bcoef_2d,val(2),iflag(2),&
                                             inbvx,inbvy,iloy,&
                                             w1_2d,w2_2d)
                         tru(2)    = f2(x(i),y(j))
                         err(2)    = abs(tru(2)-val(2))
                         errmax(2) = max(err(2),errmax(2))
            do k=1,nz
                         call db3val(x(i),y(j),z(k),idx,idy,idz,&
-                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,fcn_3d,val(3),iflag(3),&
+                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,bcoef_3d,val(3),iflag(3),&
                                             inbvx,inbvy,inbvz,iloy,iloz,&
                                             w1_3d,w2_3d,w3_3d)
                         tru(3)    = f3(x(i),y(j),z(k))
                         err(3)    = abs(tru(3)-val(3))
                         errmax(3) = max(err(3),errmax(3))
               do l=1,nq
                         call db4val(x(i),y(j),z(k),q(l),idx,idy,idz,idq,&
-                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,fcn_4d,val(4),iflag(4),&
+                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,bcoef_4d,val(4),iflag(4),&
                                             inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,&
                                             w1_4d,w2_4d,w3_4d,w4_4d)
                         tru(4)    = f4(x(i),y(j),z(k),q(l))
                         err(4)    = abs(tru(4)-val(4))
                         errmax(4) = max(err(4),errmax(4))
                 do m=1,nr
                         call db5val(x(i),y(j),z(k),q(l),r(m),idx,idy,idz,idq,idr,&
-                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,fcn_5d,val(5),iflag(5),&
+                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef_5d,val(5),iflag(5),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,&
                                             w1_5d,w2_5d,w3_5d,w4_5d,w5_5d)
                         tru(5)    = f5(x(i),y(j),z(k),q(l),r(m))
                         err(5)    = abs(tru(5)-val(5))
                         errmax(5) = max(err(5),errmax(5))
                     do n=1,ns
                         call db6val(x(i),y(j),z(k),q(l),r(m),s(n),idx,idy,idz,idq,idr,ids,&
-                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,fcn_6d,val(6),iflag(6),&
+                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,bcoef_6d,val(6),iflag(6),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,&
                                             w1_6d,w2_6d,w3_6d,w4_6d,w5_6d,w6_6d)
                         tru(6)    = f6(x(i),y(j),z(k),q(l),r(m),s(n))
diff --git a/test/bspline_test_2.f90 b/test/bspline_test_2.f90
@@ -30,12 +30,12 @@ program bspline_test_2
     real(wp) :: x(nx),y(ny),z(nz),q(nq),r(nr),s(ns)
     real(wp) :: xval(nx+2),yval(ny+2),zval(nz+2),qval(nq+2),rval(nr+2),sval(ns+2)
     real(wp) :: tx(nx+kx),ty(ny+ky),tz(nz+kz),tq(nq+kq),tr(nr+kr),ts(ns+ks)
-    real(wp) :: fcn_1d(nx)
-    real(wp) :: fcn_2d(nx,ny)
-    real(wp) :: fcn_3d(nx,ny,nz)
-    real(wp) :: fcn_4d(nx,ny,nz,nq)
-    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)
-    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns)
+    real(wp) :: fcn_1d(nx)               , bcoef_1d(nx)
+    real(wp) :: fcn_2d(nx,ny)            , bcoef_2d(nx,ny)
+    real(wp) :: fcn_3d(nx,ny,nz)         , bcoef_3d(nx,ny,nz)
+    real(wp) :: fcn_4d(nx,ny,nz,nq)      , bcoef_4d(nx,ny,nz,nq)
+    real(wp) :: fcn_5d(nx,ny,nz,nq,nr)   , bcoef_5d(nx,ny,nz,nq,nr)
+    real(wp) :: fcn_6d(nx,ny,nz,nq,nr,ns), bcoef_6d(nx,ny,nz,nq,nr,ns)
 
     real(wp),dimension(3*kx)                     :: w1_1d
     real(wp),dimension(ky)                       :: w1_2d
@@ -146,12 +146,12 @@ program bspline_test_2
 
     ! initialize
     iflag = 0
-    call db1ink(x,nx,fcn_1d,kx,iknot,tx,fcn_1d,iflag(1))
-    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,fcn_2d,iflag(2))
-    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,fcn_3d,iflag(3))
-    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,fcn_4d,iflag(4))
-    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,fcn_5d,iflag(5))
-    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,fcn_6d,iflag(6))
+    call db1ink(x,nx,fcn_1d,kx,iknot,tx,bcoef_1d,iflag(1))
+    call db2ink(x,nx,y,ny,fcn_2d,kx,ky,iknot,tx,ty,bcoef_2d,iflag(2))
+    call db3ink(x,nx,y,ny,z,nz,fcn_3d,kx,ky,kz,iknot,tx,ty,tz,bcoef_3d,iflag(3))
+    call db4ink(x,nx,y,ny,z,nz,q,nq,fcn_4d,kx,ky,kz,kq,iknot,tx,ty,tz,tq,bcoef_4d,iflag(4))
+    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,fcn_5d,kx,ky,kz,kq,kr,iknot,tx,ty,tz,tq,tr,bcoef_5d,iflag(5))
+    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,fcn_6d,kx,ky,kz,kq,kr,ks,iknot,tx,ty,tz,tq,tr,ts,bcoef_6d,iflag(6))
 
     if (any(iflag/=0)) then
         do i=1,6
@@ -168,30 +168,30 @@ program bspline_test_2
      err = -99999.9_wp
      do i=1,nx+2
                         call db1val(xval(i),idx,&
-                                            tx,nx,kx,fcn_1d,val(1),iflag(1),inbvx,w1_1d,extrap=.true.)
+                                            tx,nx,kx,bcoef_1d,val(1),iflag(1),inbvx,w1_1d,extrap=.true.)
                         tru(1)    = f1(xval(i))
                         err(1)    = abs(tru(1)-val(1))
                         errmax(1) = max(err(1),errmax(1))
                         if (iflag(1)/=0) write(*,*) 'Error evaluating 1D spline: '//get_status_message(iflag(1))
         do j=1,ny+2
                         call db2val(xval(i),yval(j),idx,idy,&
-                                            tx,ty,nx,ny,kx,ky,fcn_2d,val(2),iflag(2),&
+                                            tx,ty,nx,ny,kx,ky,bcoef_2d,val(2),iflag(2),&
                                             inbvx,inbvy,iloy,w1_2d,w2_2d,extrap=.true.)
                         tru(2)    = f2(xval(i),yval(j))
                         err(2)    = abs(tru(2)-val(2))
                         errmax(2) = max(err(2),errmax(2))
                         if (iflag(2)/=0) write(*,*) 'Error evaluating 2D spline: '//get_status_message(iflag(2))
            do k=1,nz+2
                         call db3val(xval(i),yval(j),zval(k),idx,idy,idz,&
-                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,fcn_3d,val(3),iflag(3),&
+                                            tx,ty,tz,nx,ny,nz,kx,ky,kz,bcoef_3d,val(3),iflag(3),&
                                             inbvx,inbvy,inbvz,iloy,iloz,w1_3d,w2_3d,w3_3d,extrap=.true.)
                         tru(3)    = f3(xval(i),yval(j),zval(k))
                         err(3)    = abs(tru(3)-val(3))
                         errmax(3) = max(err(3),errmax(3))
                         if (iflag(3)/=0) write(*,*) 'Error evaluating 3D spline: '//get_status_message(iflag(3))
               do l=1,nq+2
                         call db4val(xval(i),yval(j),zval(k),qval(l),idx,idy,idz,idq,&
-                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,fcn_4d,val(4),iflag(4),&
+                                            tx,ty,tz,tq,nx,ny,nz,nq,kx,ky,kz,kq,bcoef_4d,val(4),iflag(4),&
                                             inbvx,inbvy,inbvz,inbvq,iloy,iloz,iloq,&
                                             w1_4d,w2_4d,w3_4d,w4_4d,extrap=.true.)
                         tru(4)    = f4(xval(i),yval(j),zval(k),qval(l))
@@ -200,7 +200,7 @@ program bspline_test_2
                         if (iflag(4)/=0) write(*,*) 'Error evaluating 4D spline: '//get_status_message(iflag(4))
                 do m=1,nr+2
                         call db5val(xval(i),yval(j),zval(k),qval(l),rval(m),idx,idy,idz,idq,idr,&
-                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,fcn_5d,val(5),iflag(5),&
+                                            tx,ty,tz,tq,tr,nx,ny,nz,nq,nr,kx,ky,kz,kq,kr,bcoef_5d,val(5),iflag(5),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,iloy,iloz,iloq,ilor,&
                                             w1_5d,w2_5d,w3_5d,w4_5d,w5_5d,extrap=.true.)
                         tru(5)    = f5(xval(i),yval(j),zval(k),qval(l),rval(m))
@@ -209,7 +209,7 @@ program bspline_test_2
                         if (iflag(5)/=0) write(*,*) 'Error evaluating 5D spline: '//get_status_message(iflag(5))
                     do n=1,ns+2
                         call db6val(xval(i),yval(j),zval(k),qval(l),rval(m),sval(n),idx,idy,idz,idq,idr,ids,&
-                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,fcn_6d,val(6),iflag(6),&
+                                            tx,ty,tz,tq,tr,ts,nx,ny,nz,nq,nr,ns,kx,ky,kz,kq,kr,ks,bcoef_6d,val(6),iflag(6),&
                                             inbvx,inbvy,inbvz,inbvq,inbvr,inbvs,iloy,iloz,iloq,ilor,ilos,&
                                             w1_6d,w2_6d,w3_6d,w4_6d,w5_6d,w6_6d,extrap=.true.)
                         tru(6)    = f6(xval(i),yval(j),zval(k),qval(l),rval(m),sval(n))
diff --git a/test/test_integrate.f90 b/test/test_integrate.f90
@@ -41,6 +41,7 @@ program bspline_integrate_test
     real(wp)                          :: f_true !! the true integral of
                                                 !! the analytic function
     real(wp),dimension(:),allocatable :: w1_1d  !! work array
+    real(wp),dimension(nx)            :: bcoef  !! spline coefficient array
 
     do imeth = 1,2  ! the two methods
 
@@ -69,17 +70,17 @@ program bspline_integrate_test
 
             ! initialize:
             ! [note we are overwriting fcn here with the b coeffs]
-            call db1ink(x,nx,fcn,kx,iknot,tx,fcn,iflag)
+            call db1ink(x,nx,fcn,kx,iknot,tx,bcoef,iflag)
             if (iflag/=0) error stop 'error calling db1ink'
 
             ! integrate:
             if (imeth==1) then
                 if (kx>20) cycle
                 meth = 'db1sqad'
-                call db1sqad(tx,fcn,nx,kx,x1,x2,f,iflag,w1_1d)
+                call db1sqad(tx,bcoef,nx,kx,x1,x2,f,iflag,w1_1d)
             else
                 meth = 'db1fqad'
-                call db1fqad(test_function,tx,fcn,nx,kx,0_ip,x1,x2,tol,f,iflag,w1_1d)
+                call db1fqad(test_function,tx,bcoef,nx,kx,0_ip,x1,x2,tol,f,iflag,w1_1d)
             end if
 
             ! display results:
