    PROGRAM d03ubfe

!      D03UBF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : d03ubf, 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)              :: adel, aparam, ares, delmax, delmn,   &
                                          resmax, resmn, root2, x1, x2, y1,    &
                                          y2, yy, z1, z2
       INTEGER                         :: i, ifail, it, j, k, lda, n1, n2, n3, &
                                          nits, sda
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:,:), b(:,:,:), c(:,:,:),        &
                                          d(:,:,:), e(:,:,:), f(:,:,:),        &
                                          g(:,:,:), q(:,:,:), r(:,:,:),        &
                                          t(:,:,:), wrksp1(:,:,:),             &
                                          wrksp2(:,:,:), wrksp3(:,:,:), x(:),  &
                                          y(:), z(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          abs, cos, exp, max, real, sqrt
!      .. Executable Statements ..
       WRITE (nout,*) 'D03UBF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) n1, n2, n3, nits
       lda = n1
       sda = n2
       ALLOCATE (a(lda,sda,n3),b(lda,sda,n3),c(lda,sda,n3),d(lda,sda,n3), &
          e(lda,sda,n3),f(lda,sda,n3),g(lda,sda,n3),q(lda,sda,n3), &
          r(lda,sda,n3),t(lda,sda,n3),wrksp1(lda,sda,n3),wrksp2(lda,sda,n3), &
          wrksp3(lda,sda,n3),x(n1),y(n2),z(n3))

       READ (nin,*) x(1:n1)
       READ (nin,*) y(1:n2)
       READ (nin,*) z(1:n3)
       root2 = sqrt(two)
       aparam = one
!      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
       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
!      Specification for internal nodes
       DO k = 2, n3 - 1
          a(2:n1-1,2:n2-1,k) = two/((z(k)-z(k-1))*(z(k+1)-z(k-1)))
          g(2:n1-1,2:n2-1,k) = two/((z(k+1)-z(k))*(z(k+1)-z(k-1)))
       END DO
       DO j = 2, n2 - 1
          b(2:n1-1,j,2:n3-1) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1)))
          f(2:n1-1,j,2:n3-1) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1)))
       END DO
       DO i = 2, n1 - 1
          c(i,2:n2-1,2:n3-1) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1)))
          e(i,2:n2-1,2:n3-1) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1)))
       END DO
       d(1:n1,1:n2,1:n3) = -a(1:n1,1:n2,1:n3) - b(1:n1,1:n2,1:n3) - &
          c(1:n1,1:n2,1:n3) - e(1:n1,1:n2,1:n3) - f(1:n1,1:n2,1:n3) - &
          g(1:n1,1:n2,1:n3)
!      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
!      Iterative loop
       DO it = 1, nits
          resmax = zero
          resmn = zero
          DO k = 1, n3
             DO j = 1, n2
                DO i = 1, n1
                   IF (d(i,j,k)/=zero) THEN
!                     Seven point molecule formula
                      r(i,j,k) = q(i,j,k) - a(i,j,k)*t(i,j,k-1) - &
                         b(i,j,k)*t(i,j-1,k) - c(i,j,k)*t(i-1,j,k) - &
                         d(i,j,k)*t(i,j,k) - e(i,j,k)*t(i+1,j,k) - &
                         f(i,j,k)*t(i,j+1,k) - g(i,j,k)*t(i,j,k+1)
                   ELSE
!                     Explicit equation
                      r(i,j,k) = q(i,j,k) - t(i,j,k)
                   END IF
                   ares = abs(r(i,j,k))
                   resmax = max(resmax,ares)
                   resmn = resmn + ares
                END DO
             END DO
          END DO
          resmn = resmn/(real(n1*n2*n3,kind=nag_wp))

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          CALL d03ubf(n1,n2,n3,lda,sda,a,b,c,d,e,f,g,aparam,it,r,wrksp1, &
             wrksp2,wrksp3,ifail)

          IF (it==1) THEN
             WRITE (nout,99997) 'Iteration', 'Residual', 'Change'
             WRITE (nout,99996) 'No', 'Max.', 'Mean', 'Max.', 'Mean'
          END IF

!         Update the dependent variable
          delmax = zero
          delmn = zero
          DO k = 1, n3
             DO j = 1, n2
                DO i = 1, n1
                   t(i,j,k) = t(i,j,k) + r(i,j,k)
                   adel = abs(r(i,j,k))
                   delmax = max(delmax,adel)
                   delmn = delmn + adel
                END DO
             END DO
          END DO
          delmn = delmn/(real(n1*n2*n3,kind=nag_wp))
          WRITE (nout,99999) it, resmax, resmn, delmax, delmn
!      Convergence tests here if required
       END DO
!      End of iterative loop
       WRITE (nout,*)
       WRITE (nout,*) 'Table of calculated function values'
       WRITE (nout,99995)
       WRITE (nout,*)
       WRITE (nout,99998) ((k,j,(i,t(i,j,k),i=1,n1),j=1,n2),k=1,n3)

99999  FORMAT (1X,I5,4(2X,E11.4))
99998  FORMAT ((1X,I1,I3,1X,4(1X,I3,2X,F8.3)))
99997  FORMAT (1X,A,6X,A,19X,A)
99996  FORMAT (2X,A,7X,A,8X,A,11X,A,6X,A/)
99995  FORMAT (1X,'K  J',2X,4(1X,'(I      T   )'))
    END PROGRAM d03ubfe
