Program f11ddfe

!     F11DDF Example Program Text

!     Mark 29.1 Release. NAG Copyright 2023.

!     .. Use Statements ..
      Use nag_library, Only: f11bdf, f11bef, f11bff, f11ddf, f11xaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: anorm, omega, sigmax, stplhs,        &
                                          stprhs, tol
      Integer                          :: i, ifail, irevcm, iterm, itn, la,    &
                                          liwork, lwneed, lwork, m, maxitn,    &
                                          monit, n, nnz
      Character (1)                    :: ckddf, ckxaf, norm, precon, trans,   &
                                          weight
      Character (8)                    :: method
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), b(:), rdiag(:), wgt(:),        &
                                          work(:), x(:)
      Integer, Allocatable             :: icol(:), irow(:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'F11DDF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)

!     Read algorithmic parameters
      Read (nin,*) n, m
      Read (nin,*) nnz
      la = 3*nnz
      lwork = max(n*(m+3)+m*(m+5)+101,7*n+100,(2*n+m)*(m+2)+n+100,10*n+100)
      liwork = 2*n + 1
      Allocate (a(la),b(n),rdiag(n),wgt(n),work(lwork),x(n),icol(la),irow(la), &
        iwork(liwork))
      Read (nin,*) method
      Read (nin,*) precon, norm, iterm
      Read (nin,*) tol, maxitn
      Read (nin,*) anorm, sigmax
      Read (nin,*) omega

!     Read the matrix A

      Do i = 1, nnz
        Read (nin,*) a(i), irow(i), icol(i)
      End Do

!     Read right-hand side vector b and initial approximate solution x

      Read (nin,*) b(1:n)
      Read (nin,*) x(1:n)

!     Call F11BDF to initialize solver

      weight = 'N'
      monit = 0

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call f11bdf(method,precon,norm,weight,iterm,n,m,tol,maxitn,anorm,sigmax, &
        monit,lwneed,work,lwork,ifail)

!     Calculate reciprocal diagonal matrix elements if necessary

      If (precon=='P' .Or. precon=='p') Then

        iwork(1:n) = 0

        Do i = 1, nnz
          If (irow(i)==icol(i)) Then
            iwork(irow(i)) = iwork(irow(i)) + 1
            If (a(i)/=0.E0_nag_wp) Then
              rdiag(irow(i)) = 1.E0_nag_wp/a(i)
            Else
              Write (nout,*) 'Matrix has a zero diagonal element'
              Go To 100
            End If
          End If
        End Do

        Do i = 1, n
          If (iwork(i)==0) Then
            Write (nout,*) 'Matrix has a missing diagonal element'
            Go To 100
          End If
          If (iwork(i)>=2) Then
            Write (nout,*) 'Matrix has a multiple diagonal element'
            Go To 100
          End If
        End Do

      End If

!     Call F11BEF to solve the linear system

      irevcm = 0
      ckxaf = 'C'
      ckddf = 'C'

loop: Do
        ifail = 0
        Call f11bef(irevcm,x,b,wgt,work,lwork,ifail)

        Select Case (irevcm)
        Case (1)
!         Compute matrix-vector product
          trans = 'N'

          Call f11xaf(trans,n,nnz,a,irow,icol,ckxaf,x,b,ifail)

          ckxaf = 'N'
        Case (-1)
!         Compute transposed matrix-vector product
          trans = 'T'

          Call f11xaf(trans,n,nnz,a,irow,icol,ckxaf,x,b,ifail)

          ckxaf = 'N'
        Case (2)

!         SSOR preconditioning
          trans = 'N'

          Call f11ddf(trans,n,nnz,a,irow,icol,rdiag,omega,ckddf,x,b,iwork,     &
            ifail)

          ckddf = 'N'
        Case (4)
!         Termination

          ifail = 0
          Call f11bff(itn,stplhs,stprhs,anorm,sigmax,work,lwork,ifail)

          Write (nout,'(A,I10,A)') ' Converged in', itn, ' iterations'
          Write (nout,'(A,1P,E16.3)') ' Matrix norm =', anorm
          Write (nout,'(A,1P,E16.3)') ' Final residual norm =', stplhs
          Write (nout,*)

!         Output x
          Write (nout,*) '           X'
          Write (nout,'(1X,1P,E16.4)') x(1:n)
          Exit loop
        Case Default
          Exit loop
        End Select
      End Do loop

100   Continue

    End Program f11ddfe