Example description
    Program e02jdfe

!     E02JDF Example Program Text

!     Mark 27.0 Release. NAG Copyright 2019.

!     .. 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 ..
        Continue

!       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