From 2f617bc9e0e671d0a0eac5fd9a268fb524d67176 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Mon, 9 Oct 2023 21:36:28 -0500 Subject: [PATCH] test updates to avoid aliased arguments to ink routines Fixes #89 also updated the linear test. --- test/bspline_extrap_test.f90 | 6 +-- test/bspline_linear_test.f90 | 79 ++++++++++++++++++------------------ test/bspline_test.f90 | 36 ++++++++-------- test/bspline_test_2.f90 | 36 ++++++++-------- test/test_integrate.f90 | 7 ++-- 5 files changed, 83 insertions(+), 81 deletions(-) diff --git a/test/bspline_extrap_test.f90 b/test/bspline_extrap_test.f90 index 8650cefe..833fe9ee 100644 --- 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 index e73f9dbc..044e4d13 100644 --- 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 index 2367f027..d03de79d 100644 --- 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,14 +249,14 @@ 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)) @@ -264,7 +264,7 @@ program bspline_test 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)) @@ -272,7 +272,7 @@ program bspline_test 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)) @@ -280,7 +280,7 @@ program bspline_test 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)) @@ -288,7 +288,7 @@ program bspline_test 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 index 96f36ec4..9100d5d8 100644 --- 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,14 +168,14 @@ 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)) @@ -183,7 +183,7 @@ program bspline_test_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)) @@ -191,7 +191,7 @@ program bspline_test_2 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 index a307cef6..4145d7c7 100644 --- 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: