Program e02jdfe
! E02JDF Example Program Text
! Mark 30.1 Release. NAG Copyright 2024.
! .. Use Statements ..
Use nag_library, Only: e02jdf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: liopts = 100, lopts = 100, nin = 5, &
nout = 6
! .. Local Scalars ..
Integer :: gsmoothness, ifail, lcoefs, lsmaxp, &
lsminp, n, nxcels, nycels
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: coefs(:), f(:), x(:), y(:)
Real (Kind=nag_wp) :: opts(lopts), pmax(2), pmin(2)
Integer :: iopts(liopts)
! .. Intrinsic Procedures ..
Intrinsic :: maxval, minval
! .. Executable Statements ..
Write (nout,*) 'E02JDF Example Program Results'
! Generate the data to fit and set the compulsory algorithmic control
! parameters.
Call generate_data(n,x,y,f,lsminp,lsmaxp,nxcels,nycels,lcoefs,coefs, &
gsmoothness)
! Initialize the options arrays and set/get some options.
Call handle_options(iopts,liopts,opts,lopts)
! Compute the spline coefficients.
ifail = 0
Call e02jdf(n,x,y,f,lsminp,lsmaxp,nxcels,nycels,lcoefs,coefs,iopts,opts, &
ifail)
! pmin and pmax form the bounding box of the spline. We must not attempt
! to evaluate the spline outside this box.
pmin(:) = (/minval(x),minval(y)/)
pmax(:) = (/maxval(x),maxval(y)/)
Deallocate (x,y,f)
! Evaluate the approximation at a vector of values.
Call evaluate_at_vector(coefs,iopts,opts,pmin,pmax)
! Evaluate the approximation on a mesh.
Call evaluate_on_mesh(coefs,iopts,opts,pmin,pmax)
Contains
Subroutine generate_data(n,x,y,f,lsminp,lsmaxp,nxcels,nycels,lcoefs, &
coefs,gsmoothness)
! Reads n from a data file and then generates an x and a y vector of n
! pseudorandom uniformly-distributed values on (0,1]. These are passed
! to the bivariate function of R. Franke to create the data set to fit.
! The remaining input data for E02JDF are set to suitable values for
! this problem, as discussed by Davydov and Zeilfelder.
! Reads the global smoothing level from a data file. This value
! determines the minimum required length of the array of spline
! coefficients, coefs.
! .. Use Statements ..
Use nag_library, Only: g05kff, g05saf
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: lseed = 4, mstate = 21
! .. Scalar Arguments ..
Integer, Intent (Out) :: gsmoothness, lcoefs, lsmaxp, lsminp, &
n, nxcels, nycels
! .. Array Arguments ..
Real (Kind=nag_wp), Allocatable, Intent (Out) :: coefs(:), f(:), x(:), &
y(:)
! .. Local Scalars ..
Integer :: genid, ifail, lstate, subid
! .. Local Arrays ..
Integer :: seed(lseed), state(mstate)
! .. Intrinsic Procedures ..
Intrinsic :: exp
! .. Executable Statements ..
! Read the size of the data set to be generated and fitted.
! (Skip the heading in the data file.)
Read (nin,*)
Read (nin,*) n
Allocate (x(n),y(n),f(n))
! Initialize the random number generator and then generate the data.
genid = 2
subid = 53
seed(:) = (/32958,39838,881818,45812/)
lstate = mstate
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
ifail = 0
Call g05saf(n,state,x,ifail)
ifail = 0
Call g05saf(n,state,y,ifail)
! Ensure that the bounding box stretches all the way to (0,0) and (1,1).
x(1) = 0.0_nag_wp
y(1) = 0.0_nag_wp
x(n) = 1.0_nag_wp
y(n) = 1.0_nag_wp
f(:) = 0.75_nag_wp*exp(-((9._nag_wp*x(:)-2._nag_wp)**2+(9._nag_wp*y(:) &
-2._nag_wp)**2)/4._nag_wp) + 0.75_nag_wp*exp(-(9._nag_wp*x(:)+ &
1._nag_wp)**2/49._nag_wp-(9._nag_wp*y(:)+1._nag_wp)/10._nag_wp) + &
0.5_nag_wp*exp(-((9._nag_wp*x(:)-7._nag_wp)**2+(9._nag_wp*y(:)- &
3._nag_wp)**2)/4._nag_wp) - 0.2_nag_wp*exp(-(9._nag_wp*x(:)- &
4._nag_wp)**2-(9._nag_wp*y(:)-7._nag_wp)**2)
! Set the grid size for the approximation.
nxcels = 6
nycels = nxcels
! Read the required level of global smoothing.
Read (nin,*) gsmoothness
! Identify the computation.
Write (nout,*)
Write (nout,99998) 'Computing the coefficients of a C^', gsmoothness, &
' spline approximation to Franke''s function'
Write (nout,99999) 'Using a ', nxcels, ' by ', nycels, ' grid'
! Set the local-approximation control parameters.
lsminp = 3
lsmaxp = 100
! Set up the array to hold the computed spline coefficients.
Select Case (gsmoothness)
Case (1)
lcoefs = (((nxcels+2)*(nycels+2)+1)/2)*10 + 1
Case (2)
lcoefs = 28*(nxcels+2)*(nycels+2)*4 + 1
Case Default
lcoefs = 0
End Select
Allocate (coefs(lcoefs))
Return
99999 Format (1X,A,I2,A,I2,A)
99998 Format (1X,A,I1,A)
End Subroutine generate_data
Subroutine handle_options(iopts,liopts,opts,lopts)
! Auxiliary routine for initializing the options arrays and
! for demonstrating how to set and get optional parameters.
! .. Use Statements ..
Use nag_library, Only: e02zkf, e02zlf
! .. Scalar Arguments ..
Integer, Intent (In) :: liopts, lopts
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: opts(lopts)
Integer, Intent (Out) :: iopts(liopts)
! .. Local Scalars ..
Real (Kind=nag_wp) :: rvalue
Integer :: ifail, ivalue, optype
Logical :: supersmooth
Character (16) :: cvalue
Character (80) :: optstr
! .. Executable Statements ..
ifail = 0
Call e02zkf('Initialize = E02JDF',iopts,liopts,opts,lopts,ifail)
! Configure the global approximation method.
Write (optstr,99998) 'Global Smoothing Level = ', gsmoothness
ifail = 0
Call e02zkf(optstr,iopts,liopts,opts,lopts,ifail)
! If C^2 smoothing is requested, compute the spline using additional
! super-smoothness constraints?
! (The default is 'No'.)
Read (nin,*) supersmooth
If (gsmoothness==2 .And. supersmooth) Then
ifail = 0
Call e02zkf('Supersmooth C2 = Yes',iopts,liopts,opts,lopts,ifail)
End If
ifail = 0
Call e02zkf('Averaged Spline = Yes',iopts,liopts,opts,lopts,ifail)
! Configure the local approximation method.
! (The default is 'Polynomial'.)
ifail = 0
Call e02zkf('Local Method = Polynomial',iopts,liopts,opts,lopts,ifail)
Write (optstr,99999) 'Minimum Singular Value LPA = ', &
1._nag_wp/32._nag_wp
ifail = 0
Call e02zkf(optstr,iopts,liopts,opts,lopts,ifail)
Select Case (gsmoothness)
Case (1)
optstr = 'Polynomial Starting Degree = 3'
Case (2)
If (supersmooth) Then
! We can benefit from starting with local polynomials of greater
! degree than with regular C^2 smoothing.
Write (nout,*) 'Using super-smoothing'
optstr = 'Polynomial Starting Degree = 6'
Else
optstr = 'Polynomial Starting Degree = 5'
End If
End Select
ifail = 0
Call e02zkf(optstr,iopts,liopts,opts,lopts,ifail)
! As an example of how to get the value of an optional parameter,
! display whether averaging of local approximations is in operation.
ifail = 0
Call e02zlf('Averaged Spline',ivalue,rvalue,cvalue,optype,iopts,opts, &
ifail)
If (cvalue=='YES') Then
Write (nout,*) 'Using an averaged local approximation'
End If
Return
99999 Format (A,E16.9)
99998 Format (A,I1)
End Subroutine handle_options
Subroutine evaluate_at_vector(coefs,iopts,opts,pmin,pmax)
! Evaluates the approximation at a vector of values.
! .. Use Statements ..
Use nag_library, Only: e02jef
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: coefs(*), opts(*), pmax(2), pmin(2)
Integer, Intent (In) :: iopts(*)
! .. Local Scalars ..
Integer :: i, ifail, nevalv
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: fevalv(:), xevalv(:), yevalv(:)
! .. Intrinsic Procedures ..
Intrinsic :: max, min
! .. Executable Statements ..
Read (nin,*) nevalv
Allocate (xevalv(nevalv),yevalv(nevalv),fevalv(nevalv))
Read (nin,*)(xevalv(i),yevalv(i),i=1,nevalv)
! Force the points to be within the bounding box of the spline.
Do i = 1, nevalv
xevalv(i) = max(xevalv(i),pmin(1))
xevalv(i) = min(xevalv(i),pmax(1))
yevalv(i) = max(yevalv(i),pmin(2))
yevalv(i) = min(yevalv(i),pmax(2))
End Do
ifail = 0
Call e02jef(nevalv,xevalv,yevalv,coefs,fevalv,iopts,opts,ifail)
Write (nout,*)
Write (nout,*) 'Values of computed spline at (x_i,y_i):'
Write (nout,*)
Write (nout,99999) 'x_i', 'y_i', 'f(x_i,y_i)'
Write (nout,99998)(xevalv(i),yevalv(i),fevalv(i),i=1,nevalv)
Return
99999 Format (1X,3A12)
99998 Format (1X,3F12.2)
End Subroutine evaluate_at_vector
Subroutine evaluate_on_mesh(coefs,iopts,opts,pmin,pmax)
! Evaluates the approximation on a mesh of n_x * n_y values.
! .. Use Statements ..
Use nag_library, Only: e02jff
! .. Implicit None Statement ..
Implicit None
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: coefs(*), opts(*), pmax(2), pmin(2)
Integer, Intent (In) :: iopts(*)
! .. Local Scalars ..
Integer :: i, ifail, j, nxeval, nyeval
Logical :: print_mesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: fevalm(:,:), xevalm(:), yevalm(:)
Real (Kind=nag_wp) :: h(2), ll_corner(2), ur_corner(2)
! .. Intrinsic Procedures ..
Intrinsic :: max, min, real
! .. Executable Statements ..
Read (nin,*) nxeval, nyeval
Allocate (xevalm(nxeval),yevalm(nyeval),fevalm(nxeval,nyeval))
! Define the mesh by its lower-left and upper-right corners.
Read (nin,*) ll_corner(1:2)
Read (nin,*) ur_corner(1:2)
! Set the mesh spacing and the evaluation points.
! Force the points to be within the bounding box of the spline.
h(1) = (ur_corner(1)-ll_corner(1))/real(nxeval-1,nag_wp)
h(2) = (ur_corner(2)-ll_corner(2))/real(nyeval-1,nag_wp)
Do i = 1, nxeval
xevalm(i) = ll_corner(1) + real(i-1,nag_wp)*h(1)
xevalm(i) = max(xevalm(i),pmin(1))
xevalm(i) = min(xevalm(i),pmax(1))
End Do
Do j = 1, nyeval
yevalm(j) = ll_corner(2) + real(j-1,nag_wp)*h(2)
yevalm(j) = max(yevalm(j),pmin(2))
yevalm(j) = min(yevalm(j),pmax(2))
End Do
! Evaluate.
ifail = 0
Call e02jff(nxeval,nyeval,xevalm,yevalm,coefs,fevalm,iopts,opts,ifail)
! Output the computed function values?
Read (nin,*) print_mesh
If (.Not. print_mesh) Then
Write (nout,*)
Write (nout,*) &
'Outputting of the function values on the mesh is disabled'
Else
Write (nout,*)
Write (nout,*) 'Values of computed spline at (x_i,y_j):'
Write (nout,*)
Write (nout,99999) 'x_i', 'y_j', 'f(x_i,y_j)'
Write (nout,99998)((xevalm(i),yevalm(j),fevalm(i, &
j),i=1,nxeval),j=1,nyeval)
End If
Return
99999 Format (1X,3A12)
99998 Format (1X,3F12.2)
End Subroutine evaluate_on_mesh
End Program e02jdfe