Program d03ebfe

!     D03EBF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. Use Statements ..
      Use nag_library, Only: d03ebf, 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
      Integer                          :: i, ifail, itcoun, itmax, itused,     &
                                          ixn, iyn, j, lda, n1, n2, ndir
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), b(:,:), c(:,:), chngs(:),    &
                                          d(:,:), e(:,:), q(:,:), resids(:),   &
                                          t(:,:), wrksp1(:,:), wrksp2(:,:),    &
                                          wrksp3(:,:), x(:), y(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cos, exp
!     .. Executable Statements ..
      Write (nout,*) 'D03EBF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n1, n2, itmax
      lda = n1
      Allocate (a(lda,n2),b(lda,n2),c(lda,n2),chngs(itmax),d(lda,n2),          &
        e(lda,n2),q(lda,n2),resids(itmax),t(lda,n2),wrksp1(lda,n2),            &
        wrksp2(lda,n2),wrksp3(lda,n2),x(n1),y(n2))

      Read (nin,*) x(1:n1)
      Read (nin,*) y(1:n2)
      Read (nin,*) conres, conchn
      Read (nin,*) ndir
      aparam = one
      itcoun = 0
!     Set up difference equation coefficients, source terms and
!     initial conditions.
      a(1:n1,1:n2) = zero
      b(1:n1,1:n2) = zero
      c(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
!     Non-zero specification for internal nodes
      Do j = 2, n2 - 1
        Do i = 2, n1 - 1
          a(i,j) = two/((y(j)-y(j-1))*(y(j+1)-y(j-1)))
          e(i,j) = two/((y(j+1)-y(j))*(y(j+1)-y(j-1)))
          b(i,j) = two/((x(i)-x(i-1))*(x(i+1)-x(i-1)))
          d(i,j) = two/((x(i+1)-x(i))*(x(i+1)-x(i-1)))
          c(i,j) = -a(i,j) - b(i,j) - d(i,j) - e(i,j)
        End Do
      End Do
!     Non-zero specification for boundary nodes
      q(1:n1,1) = exp((x(1:n1)+one)/y(n2))*cos(y(1)/y(n2))
      q(1:n1,n2) = exp((x(1:n1)+one)/y(n2))*cos(y(n2)/y(n2))
      q(1,1:n2) = exp((x(1)+one)/y(n2))*cos(y(1:n2)/y(n2))
      q(n1,1:n2) = exp((x(n1)+one)/y(n2))*cos(y(1:n2)/y(n2))

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

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

99999 Format (2X,I2,10X,E11.4,4X,E11.4)
99998 Format (1X,' I',1X,6(I4,7X))
99997 Format (1X,I2,1X,6(F9.3,2X))
    End Program d03ebfe