!   G05ZSF Example Program Text

!   Mark 26.1 Release. NAG Copyright 2016.

    Program g05zsfe

!     G05ZSF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: g05zrf, g05zsf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: lenst = 17, nin = 5, nout = 6,       &
                                          npmax = 4
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: rho, var, xmax, xmin, ymax, ymin
      Integer                          :: approx, icorr, icount, icov2, ifail, &
                                          norm, np, pad, s
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: eig(3), params(npmax)
      Real (Kind=nag_wp), Allocatable  :: lam(:), xx(:), yy(:), z(:,:)
      Integer                          :: m(2), maxm(2), ns(2), state(lenst)
!     .. Executable Statements ..
      Write (nout,*) 'G05ZSF Example Program Results'
      Write (nout,*)
      Flush (nout)

!     Get problem specifications from data file
      Call read_input_data(icov2,np,params,norm,var,xmin,xmax,ymin,ymax,ns,    &
        maxm,icorr,pad,s)

      Allocate (lam(maxm(1)*maxm(2)),xx(ns(1)),yy(ns(2)))

!     Get square roots of the eigenvalues of the embedding matrix
      ifail = 0
      Call g05zrf(ns,xmin,xmax,ymin,ymax,maxm,var,icov2,norm,np,params,pad,    &
        icorr,lam,xx,yy,m,approx,rho,icount,eig,ifail)

      Call display_embedding_results(approx,m,rho,eig,icount)

!     Initialize state array
      Call initialize_state(state)

      Allocate (z(ns(1)*ns(2),s))

!     Compute s random field realizations
      ifail = 0
      Call g05zsf(ns,s,m,lam,rho,state,z,ifail)

      Call display_realizations(ns,s,xx,yy,z)

    Contains
      Subroutine read_input_data(icov2,np,params,norm,var,xmin,xmax,ymin,ymax, &
        ns,maxm,icorr,pad,s)

!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: var, xmax, xmin, ymax, ymin
        Integer, Intent (Out)          :: icorr, icov2, norm, np, pad, s
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: params(npmax)
        Integer, Intent (Out)          :: maxm(2), ns(2)
!       .. Executable Statements ..
!       Skip heading in data file
        Read (nin,*)

!       Read in covariance function number
        Read (nin,*) icov2

!       Read in number of parameters
        Read (nin,*) np

!       Read in parameters
        If (np>0) Then
          Read (nin,*) params(1:np)
        End If

!       Read in choice of norm to use
        Read (nin,*) norm

!       Read in variance of random field
        Read (nin,*) var

!       Read in domain endpoints
        Read (nin,*) xmin, xmax
        Read (nin,*) ymin, ymax

!       Read in number of sample points
        Read (nin,*) ns(1:2)

!       Read in maximum size of embedding matrix
        Read (nin,*) maxm(1:2)

!       Read in choice of scaling in case of approximation
        Read (nin,*) icorr

!       Read in choice of padding
        Read (nin,*) pad

!       Read in number of realization samples to be generated
        Read (nin,*) s

        Return

      End Subroutine read_input_data

      Subroutine display_embedding_results(approx,m,rho,eig,icount)

!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: rho
        Integer, Intent (In)           :: approx, icount
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: eig(3)
        Integer, Intent (In)           :: m(2)
!       .. Executable Statements ..
!       Display size of embedding matrix
        Write (nout,*)
        Write (nout,99999) 'Size of embedding matrix = ', m(1)*m(2)

!       Display approximation information if approximation used
        Write (nout,*)
        If (approx==1) Then
          Write (nout,*) 'Approximation required'
          Write (nout,*)
          Write (nout,99998) 'RHO = ', rho
          Write (nout,99997) 'EIG = ', eig(1:3)
          Write (nout,99999) 'ICOUNT = ', icount
        Else
          Write (nout,*) 'Approximation not required'
        End If

        Return

99999   Format (1X,A,I7)
99998   Format (1X,A,F10.5)
99997   Format (1X,A,3(F10.5,1X))

      End Subroutine display_embedding_results

      Subroutine initialize_state(state)

!       .. Use Statements ..
        Use nag_library, Only: g05kff
!       .. Implicit None Statement ..
        Implicit None
!       .. Parameters ..
        Integer, Parameter             :: genid = 1, inseed = 14965,           &
                                          lseed = 1, subid = 1
!       .. Array Arguments ..
        Integer, Intent (Out)          :: state(lenst)
!       .. Local Scalars ..
        Integer                        :: ifail, lstate
!       .. Local Arrays ..
        Integer                        :: seed(lseed)
!       .. Executable Statements ..
!       Initialize the generator to a repeatable sequence
        lstate = lenst
        seed(1) = inseed
        ifail = 0
        Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

      End Subroutine initialize_state

      Subroutine display_realizations(ns,s,xx,yy,z)

!       .. Use Statements ..
        Use nag_library, Only: x04cbf
!       .. Implicit None Statement ..
        Implicit None
!       .. Parameters ..
        Integer, Parameter             :: indent = 0, ncols = 80
        Character (1), Parameter       :: charlab = 'C', intlab = 'I',         &
                                          matrix = 'G', unit = 'n'
        Character (5), Parameter       :: form = 'F10.5'
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: s
!       .. Array Arguments ..
        Integer, Intent (In)           :: ns(2)
        Real (Kind=nag_wp), Intent (In) :: xx(ns(1)), yy(ns(2)),               &
                                          z(ns(1)*ns(2),s)
!       .. Local Scalars ..
        Integer                        :: i, ifail, j, nn
        Character (61)                 :: title
!       .. Local Arrays ..
        Character (1)                  :: clabs(0)
        Character (12), Allocatable    :: rlabs(:)
!       .. Executable Statements ..
        nn = ns(1)*ns(2)
        Allocate (rlabs(nn))

!       Set row labels to grid points (column label is realization number).
        Do j = 1, ns(2)
          Do i = 1, ns(1)
            If (i==1) Then
              Write (rlabs((j-1)*ns(1)+i),99999) xx(i), yy(j)
            Else
              Write (rlabs((j-1)*ns(1)+i),99998) xx(i)
            End If
          End Do
        End Do

!       Display random field results
        title = 'Random field realizations (x,y coordinates first):'
        Write (nout,*)
        Flush (nout)
        ifail = 0
        Call x04cbf(matrix,unit,nn,s,z,nn,form,title,charlab,rlabs,intlab,     &
          clabs,ncols,indent,ifail)

99999   Format (2F6.1)
99998   Format (F6.1,5X,'.')

      End Subroutine display_realizations

    End Program g05zsfe