Program d03ecfe

!     D03ECF Example Program Text

!     Mark 30.3 Release. nAG Copyright 2024.

!     .. Use Statements ..
      Use nag_library, Only: d03ecf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: aparam, conchn, conres, root2, x1,   &
                                          x2, y1, y2, yy, z1, z2
      Integer                          :: i, ifail, itcoun, itmax, itused,     &
                                          ixn, iyn, izn, j, k, lda, n1, n2,    &
                                          n3, ndir, sda
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:,:), b(:,:,:), c(:,:,:),        &
                                          chngs(:), d(:,:,:), e(:,:,:),        &
                                          f(:,:,:), g(:,:,:), q(:,:,:),        &
                                          resids(:), t(:,:,:), wrksp1(:,:,:),  &
                                          wrksp2(:,:,:), wrksp3(:,:,:),        &
                                          wrksp4(:,:,:), x(:), y(:), z(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cos, exp, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'D03ECF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n1, n2, n3, itmax
      lda = n1
      sda = n2
      Allocate (a(lda,sda,n3),b(lda,sda,n3),c(lda,sda,n3),chngs(itmax),        &
        d(lda,sda,n3),e(lda,sda,n3),f(lda,sda,n3),g(lda,sda,n3),q(lda,sda,n3), &
        resids(itmax),t(lda,sda,n3),wrksp1(lda,sda,n3),wrksp2(lda,sda,n3),     &
        wrksp3(lda,sda,n3),wrksp4(lda,sda,n3),x(n1),y(n2),z(n3))

      Read (nin,*) x(1:n1)
      Read (nin,*) y(1:n2)
      Read (nin,*) z(1:n3)
      Read (nin,*) conres, conchn
      Read (nin,*) ndir
      root2 = sqrt(two)
      aparam = one
      itcoun = 0
!     Set up difference equation coefficients, source terms and
!     initial approximation.
      a(1:n1,1:n2,1:n3) = zero
      b(1:n1,1:n2,1:n3) = zero
      c(1:n1,1:n2,1:n3) = zero
      d(1:n1,1:n2,1:n3) = zero
      e(1:n1,1:n2,1:n3) = zero
      f(1:n1,1:n2,1:n3) = zero
      g(1:n1,1:n2,1:n3) = zero
      q(1:n1,1:n2,1:n3) = zero
      t(1:n1,1:n2,1:n3) = zero
!     Non-zero Specification for internal nodes
      Do k = 2, n3 - 1
        Do j = 2, n2 - 1
          Do i = 2, n1 - 1
            a(i,j,k) = two/((z(k)-z(k-1))*(z(k+1)-z(k-1)))
            g(i,j,k) = two/((z(k+1)-z(k))*(z(k+1)-z(k-1)))
            b(i,j,k) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1)))
            f(i,j,k) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1)))
            c(i,j,k) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1)))
            e(i,j,k) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1)))
            d(i,j,k) = -a(i,j,k) - b(i,j,k) - c(i,j,k) - e(i,j,k) - f(i,j,k) - &
              g(i,j,k)
          End Do
        End Do
      End Do
!     Non-zero specification for boundary nodes
      yy = one/y(n2)
      x1 = (x(1)+one)*yy
      x2 = (x(n1)+one)*yy
      Do j = 1, n2
        y1 = root2*y(j)*yy
        q(1,j,1:n3) = exp(x1)*cos(y1)*exp((-z(1:n3)-one)*yy)
        q(n1,j,1:n3) = exp(x2)*cos(y1)*exp((-z(1:n3)-one)*yy)
      End Do
      y1 = root2*y(1)*yy
      y2 = root2*y(n2)*yy
      Do i = 1, n1
        x1 = (x(i)+one)*yy
        q(i,1,1:n3) = exp(x1)*cos(y1)*exp((-z(1:n3)-one)*yy)
        q(i,n2,1:n3) = exp(x1)*cos(y2)*exp((-z(1:n3)-one)*yy)
      End Do
      z1 = (-z(1)-one)*yy
      z2 = (-z(n3)-one)*yy
      Do i = 1, n1
        x1 = (x(i)+one)*yy
        q(i,1:n2,1) = exp(x1)*cos(root2*y(1:n2)*yy)*exp(z1)
        q(i,1:n2,n3) = exp(x1)*cos(root2*y(1:n2)*yy)*exp(z2)
      End Do

!     ifail: behaviour on error exit
!            =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call d03ecf(n1,n2,n3,lda,sda,a,b,c,d,e,f,g,q,t,aparam,itmax,itcoun,      &
        itused,ndir,ixn,iyn,izn,conres,conchn,resids,chngs,wrksp1,wrksp2,      &
        wrksp3,wrksp4,ifail)

      Write (nout,*) 'Iteration      Maximum        Maximum'
      Write (nout,*) ' number        residual        change'
      If (itused/=0) Then
        Write (nout,99999)(i,resids(i),chngs(i),i=1,itused)
      End If
      Write (nout,*)
      Write (nout,*) 'Table of calculated function values'
      Write (nout,*)
      Write (nout,99998)
      Do k = 1, n3
        Do j = 1, n2
          Write (nout,99997) k, j, (i,t(i,j,k),i=1,n1)
        End Do
      End Do

99999 Format (2X,I3,9X,E11.4,4X,E11.4)
99998 Format (1X,'K  J  ',4(' (I      T   )'))
99997 Format (1X,I1,2X,I1,1X,4(1X,I3,2X,F8.3))
    End Program d03ecfe