!   F11DFF Example Program
!
!   NAG FORTRAN Library.
!   Mark 25 Release. NAG Copyright 2014.

    Program f11dffe

!     .. Use Statements ..
      Use nag_library, Only: f11dff, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: dtolg
      Integer                          :: i, ifail, k, la, lfillg, lindb,      &
                                          liwork, mb, n, nb, nnz, nnzc, nover
      Character (1)                    :: milug, pstrag
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:), b(:), dtol(:)
      Integer, Allocatable             :: icol(:), idiag(:), indb(:),          &
                                          ipivp(:), ipivq(:), irow(:),         &
                                          istb(:), istr(:), iwork(:),          &
                                          lfill(:), npivm(:)
      Character (1), Allocatable       :: milu(:), pstrat(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: maxval, minval
!     .. Executable Statements ..
      Continue

!     Print example header      
      Write (nout,*) 'F11DFF Example Program Results'
      Write (nout,*)

!     Skip heading in data file
      Read (nin,*)

!     Get the matrix order and number of non-zero entries.
      Read (nin,*) n
      Read (nin,*) nnz

!     Allocate arrays with lengths based on n and nnz..
      liwork = 9*n + 3
      Allocate (b(n),iwork(liwork))

      la = 20*nnz
      Allocate (a(la),irow(la),icol(la))

      lindb = 3*n
      Allocate (idiag(lindb),indb(lindb),ipivp(lindb),ipivq(lindb), &
        istr(lindb+1))

!     Read the matrix A
      Read (nin,*)(a(i),irow(i),icol(i),i=1,nnz)

!     Read algorithmic parameters
      Read (nin,*) lfillg, dtolg
      Read (nin,*) pstrag
      Read (nin,*) milug
      Read (nin,*) nb, nover

!     Allocate arrays with length based on number of blocks.
      Allocate (dtol(nb),istb(nb+1),lfill(nb),npivm(nb),milu(nb),pstrat(nb))

!     Define diagonal block indices. 
!     In this example use blocks of MB consecutive rows and initialise 
!     assuming no overlap.
      mb = (n+nb-1)/nb
      Do k = 1, nb
        istb(k) = (k-1)*mb + 1
      End Do
      istb(nb+1) = n + 1
      Do i = 1, n
        indb(i) = i
      End Do

!     Modify INDB and ISTB to account for overlap.
      Call f11dffe_overlap(n,nnz,la,irow,icol,nb,istb,indb,lindb,nover,iwork)
      If (iwork(1)==-999) Then
        Write (nout,*) '** LINDB too small, LINDB = ', lindb, '.'
        Go To 100
      End If

!     Output matrix and blocking details
      Write (nout,*) ' Original Matrix'
      Write (nout,99997) ' N     =', n
      Write (nout,99997) ' NNZ   =', nnz
      Write (nout,99997) ' NB    =', nb
      Do k = 1, nb
        Write (nout,99993) ' Block ', k, ': order = ', istb(k+1) - istb(k), &
          ', start row = ', minval(indb(istb(k):istb(k+1)-1))
      End Do

!     Set algorithmic parameters for each block from global values
      lfill(1:nb) = lfillg
      dtol(1:nb) = dtolg
      pstrat(1:nb) = pstrag
      milu(1:nb) = milug

!     Calculate factorization
      ifail = 1
      Call f11dff(n,nnz,a,la,irow,icol,nb,istb,indb,lindb,lfill,dtol,pstrat, &
        milu,ipivp,ipivq,istr,idiag,nnzc,npivm,iwork,liwork,ifail)
      If (ifail/=0) Then
        Write (nout,99995) ifail
        Go To 100
      End If

!     Output details of the factorization
      Write (nout,99996) ' Factorization'
      Write (nout,99997) '  NNZC  =', nnzc

      Write (nout,99996) ' Elements of factorization'
      Write (nout,99996) '        I   J        C(I,J)     Index'
      Do k = 1, nb
        Write (nout,99994) ' C_', k, ' --------------------------------'
!       Elements of the k-th block
        Do i = istr(istb(k)), istr(istb(k+1)) - 1
          Write (nout,99992) irow(i), icol(i), a(i), i
        End Do
      End Do

      Write (nout,99996) ' Details of factorized blocks'
      If (maxval(npivm(1:nb))>0) Then
!       Including pivoting details.
        Write (nout,99996) &
          '  K   I      ISTR(I)  IDIAG(I)   INDB(I)  IPIVP(I)  IPIVQ(I)'
        Do k = 1, nb
          i = istb(k)
          Write (nout,99999) k, i, istr(i), idiag(i), indb(i), ipivp(i), &
            ipivq(i)
          Do i = istb(k) + 1, istb(k+1) - 1
            Write (nout,99998) i, istr(i), idiag(i), indb(i), ipivp(i), &
              ipivq(i)
          End Do
          Write (nout,*) &
            ' -----------------------------------------------------'
        End Do
      Else
!       No pivoting on any block.
        Write (nout,99996) '  K   I      ISTR(I)  IDIAG(I)   INDB(I)'
        Do k = 1, nb
          i = istb(k)
          Write (nout,99999) k, i, istr(i), idiag(i), indb(i)
          Do i = istb(k) + 1, istb(k+1) - 1
            Write (nout,99998) i, istr(i), idiag(i), indb(i)
          End Do
          Write (nout,*) ' ------------------------------------'
        End Do
      End If

100   Continue

99999 Format (1X,I3,1X,I3,5(I10))
99998 Format (1X,I7,5(I10))
99997 Format (1X,A,I4)
99996 Format (1X/1X,A)
99995 Format (1X/1X,' ** F11DFF returned with IFAIL = ',I5)
99994 Format (1X,A3,I1,A)
99993 Format (1X,A,I3,A,I3,A,I3)
99992 Format (6X,2I4,E16.5,I8)

    Contains

      Subroutine f11dffe_overlap(n,nnz,la,irow,icol,nb,istb,indb,lindb,nover, &
        iwork)

!       Purpose
!       =======
!       This routine takes a set of row indices INDB defining the diagonal 
!       blocks to be used in F11DFF to define a block Jacobi or additive Schwarz 
!       preconditioner, and expands them to allow for NOVER levels of overlap.
!       The pointer array ISTB is also updated accordingly, so that the returned
!       values of ISTB and INDB can be passed to F11DFF to define overlapping 
!       diagonal blocks.
!       ------------------------------------------------------------------------

!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In)             :: la, lindb, n, nb, nnz, nover
!       .. Array Arguments ..
        Integer, Intent (In)             :: icol(la), irow(la)
        Integer, Intent (Inout)          :: indb(lindb), istb(nb+1)
        Integer, Intent (Out)            :: iwork(3*n+1)
!       .. Local Scalars ..
        Integer                          :: i, ik, ind, iover, k, l, n21,      &
                                            nadd, row
!       .. Executable Statements ..
        Continue

!       Find the number of non-zero elements in each row of the matrix A, and
!       and start address of each row. Store the start addresses in 
!       IWORK(N+1,...,2*N+1).

        iwork(1:n) = 0
        Do k = 1, nnz
          iwork(irow(k)) = iwork(irow(k)) + 1
        End Do
        iwork(n+1) = 1
        Do i = 1, n
          iwork(n+i+1) = iwork(n+i) + iwork(i)
        End Do

!       Loop over blocks.
blocks: Do k = 1, nb

!         Initialize marker array.
          iwork(1:n) = 0

!         Mark the rows already in block K in the workspace array.
          Do l = istb(k), istb(k+1) - 1
            iwork(indb(l)) = 1
          End Do

!         Loop over levels of overlap.
          Do iover = 1, nover

!           Initialize counter of new row indices to be added.
            ind = 0

!           Loop over the rows currently in the diagonal block.
            Do l = istb(k), istb(k+1) - 1
              row = indb(l)

!             Loop over non-zero elements in row ROW.
              Do i = iwork(n+row), iwork(n+row+1) - 1

!               If the column index of the non-zero element is not in the
!               existing set for this block, store it to be added later, and
!               mark it in the marker array.
                If (iwork(icol(i))==0) Then
                  iwork(icol(i)) = 1
                  ind = ind + 1
                  iwork(2*n+1+ind) = icol(i)
                End If
              End Do
            End Do

!           Shift the indices in INDB and add the new entries for block K.
!           Change ISTB accordingly.
            nadd = ind
            If (istb(nb+1)+nadd-1>lindb) Then
              iwork(1) = -999
              Exit blocks
            End If

            Do i = istb(nb+1) - 1, istb(k+1), -1
              indb(i+nadd) = indb(i)
            End Do
            n21 = 2*n + 1
            ik = istb(k+1) - 1
            indb(ik+1:ik+nadd) = iwork(n21+1:n21+nadd)
            istb(k+1:nb+1) = istb(k+1:nb+1) + nadd
          End Do
        End Do blocks

        Return

      End Subroutine f11dffe_overlap
    End Program f11dffe