Program d03uafe

!     D03UAF Example Program Text

!     Mark 28.3 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: d03uaf, 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
      Integer                          :: i, ifail, it, j, lda, n1, n2, nits
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), b(:,:), c(:,:), d(:,:),      &
                                          e(:,:), q(:,:), r(:,:), t(:,:),      &
                                          wrksp1(:,:), wrksp2(:,:), x(:), y(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, cos, exp, max, real
!     .. Executable Statements ..
      Write (nout,*) 'D03UAF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n1, n2, nits
      lda = n1
      Allocate (a(lda,n2),b(lda,n2),c(lda,n2),d(lda,n2),e(lda,n2),q(lda,n2),   &
        r(lda,n2),t(lda,n2),wrksp1(lda,n2),wrksp2(lda,n2),x(n1),y(n2))
      Read (nin,*) x(1:n1)
      Read (nin,*) y(1:n2)
      aparam = one

!     Set up difference equation coefficients, source terms and
!     initial S
      a(1:n1,1:n2) = zero
      b(1:n1,1:n2) = zero
      d(1:n1,1:n2) = zero
      e(1:n1,1:n2) = zero
      q(1:n1,1:n2) = zero
      t(1:n1,1:n2) = zero
!     Specification for internal nodes
      Do j = 2, n2 - 1
        a(2:n1-1,j) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1)))
        e(2:n1-1,j) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1)))
      End Do
      Do i = 2, n1 - 1
        b(i,2:n2-1) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1)))
        d(i,2:n2-1) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1)))
      End Do
      c(1:n1,1:n2) = -a(1:n1,1:n2) - b(1:n1,1:n2) - d(1:n1,1:n2) -             &
        e(1:n1,1:n2)
!     Specification for boundary nodes
      Do j = 1, n2
        q(1,j) = exp((x(1)+one)/y(n2))*cos(y(j)/y(n2))
        q(n1,j) = exp((x(n1)+one)/y(n2))*cos(y(j)/y(n2))
      End Do
      Do i = 1, n1
        q(i,1) = exp((x(i)+one)/y(n2))*cos(y(1)/y(n2))
        q(i,n2) = exp((x(i)+one)/y(n2))*cos(y(n2)/y(n2))
      End Do

!     Iterative loop
      Do it = 1, nits
!       Calculate the residuals
        resmax = zero
        resmn = zero
        Do j = 1, n2
          Do i = 1, n1
            If (c(i,j)/=zero) Then
!             Five point molecule formula
              r(i,j) = q(i,j) - a(i,j)*t(i,j-1) - b(i,j)*t(i-1,j) -            &
                c(i,j)*t(i,j) - d(i,j)*t(i+1,j) - e(i,j)*t(i,j+1)
            Else
!             Explicit equation
              r(i,j) = q(i,j) - t(i,j)
            End If
            ares = abs(r(i,j))
            resmax = max(resmax,ares)
            resmn = resmn + ares
          End Do
        End Do
        resmn = resmn/(real(n1*n2,kind=nag_wp))

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d03uaf(n1,n2,lda,a,b,c,d,e,aparam,it,r,wrksp1,wrksp2,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 j = 1, n2
          Do i = 1, n1
            t(i,j) = t(i,j) + r(i,j)
            adel = abs(r(i,j))
            delmax = max(delmax,adel)
            delmn = delmn + adel
          End Do
        End Do
        delmn = delmn/(real(n1*n2,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,*)
      Write (nout,99995) 'I', 1, (i,i=2,6)
      Write (nout,*) ' J'
      Do j = 1, n2
        Write (nout,99998) j, (t(i,j),i=1,n1)
      End Do

99999 Format (1X,I3,4(2X,E11.4))
99998 Format (1X,I2,1X,6(F9.3,2X))
99997 Format (1X,A,6X,A,19X,A)
99996 Format (3X,A,7X,A,8X,A,11X,A,6X,A,/)
99995 Format (4X,A,4X,I1,5I11)
    End Program d03uafe