Program d03fafe

!     D03FAF Example Program Text

!     Mark 28.3 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: d03faf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: dx, dy, dz, error, lambda, pertrb,   &
                                          t, xf, xs, yf, ys, yy, zf, zs, zz
      Integer                          :: i, ifail, j, k, l, lbdcnd, ldf,      &
                                          ldf2, lwrk, m, mbdcnd, n, nbdcnd
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: bdxf(:,:), bdxs(:,:), bdyf(:,:),     &
                                          bdys(:,:), bdzf(:,:), bdzs(:,:),     &
                                          f(:,:,:), x(:), y(:), z(:)
      Real (Kind=nag_wp)               :: w(0)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, cos, real, sin
!     .. Executable Statements ..
      Write (nout,*) 'D03FAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) l, m, n
      ldf = l + 1
      ldf2 = m + 1
      lwrk = 0
      Allocate (bdxf(ldf2,n+1),bdxs(ldf2,n+1),bdyf(ldf,n+1),bdys(ldf,n+1),     &
        bdzf(ldf,m+1),bdzs(ldf,m+1),f(ldf,ldf2,n+1),x(l+1),y(m+1),z(n+1))
      Read (nin,*) lambda
      Read (nin,*) xs, xf
      Read (nin,*) ys, yf
      Read (nin,*) zs, zf
      Read (nin,*) lbdcnd, mbdcnd, nbdcnd

!     Define the grid points for later use.

      dx = (xf-xs)/real(l,kind=nag_wp)
      Do i = 1, l + 1
        x(i) = xs + real(i-1,kind=nag_wp)*dx
      End Do
      dy = (yf-ys)/real(m,kind=nag_wp)
      Do j = 1, m + 1
        y(j) = ys + real(j-1,kind=nag_wp)*dy
      End Do
      dz = (zf-zs)/real(n,kind=nag_wp)
      Do k = 1, n + 1
        z(k) = zs + real(k-1,kind=nag_wp)*dz
      End Do

!     Define the array of derivative boundary values.

      Do j = 1, m + 1
        yy = sin(y(j))
        bdzf(1:l+1,j) = -yy*x(1:l+1)**4
      End Do

!     Note that for this example all other boundary arrays are
!     dummy variables.

!     We define the function boundary values in the F array.

      f(1,1:m+1,1:n+1) = 0.0_nag_wp
      Do k = 1, n + 1
        zz = cos(z(k))
        Do j = 1, m + 1
          f(l+1,j,k) = zz*sin(y(j))
        End Do
      End Do
      f(1:l+1,1:m+1,1) = -bdzf(1:l+1,1:m+1)

!     Define the values of the right hand side of the Helmholtz
!     equation.

      Do k = 2, n + 1
        zz = 4.0_nag_wp*cos(z(k))
        Do j = 1, m + 1
          yy = sin(y(j))*zz
          Do i = 2, l
            f(i,j,k) = yy*x(i)**2*(3.0_nag_wp-x(i)**2)
          End Do
        End Do
      End Do

!     Call D03FAF to generate and solve the finite difference equation.
!     ifail: behaviour on error exit
!            =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call d03faf(xs,xf,l,lbdcnd,bdxs,bdxf,ys,yf,m,mbdcnd,bdys,bdyf,zs,zf,n,   &
        nbdcnd,bdzs,bdzf,lambda,ldf,ldf2,f,pertrb,w,lwrk,ifail)

!     Compute discretization error.  The exact solution to the
!     problem is

!     U(X,Y,Z) = X**4*SIN(Y)*COS(Z)

      error = 0.0_nag_wp
      Do k = 1, n + 1
        zz = cos(z(k))
        Do j = 1, m + 1
          yy = sin(y(j))*zz
          Do i = 1, l + 1
            t = abs(f(i,j,k)-yy*x(i)**4)
            If (t>error) Then
              error = t
            End If
          End Do
        End Do
      End Do
      Write (nout,*)
      Write (nout,99999) error

99999 Format (1X,'Maximum component of discretization error =',1P,E13.6)
    End Program d03fafe