NAG Library Manual, Mark 28.5
Interfaces:  FL   CL   CPP   AD
```!   F11DGF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.

Program f11dgfe

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

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

!     Skip heading in data file

!     Get the square matrix size

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

!     Get the number of non zero (nnz) matrix entries
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 in matrix A

!     Read in RHS

!     Read algorithmic parameters
Read (nin,*) lfillg, dtolg
Read (nin,*) m, tol, maxitn
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))

!     Set up initial approximate solution x
x(1:n) = 0.0_nag_wp

!     Define diagonal block indices.
!     In this example use blocks of MB consecutive rows and initialize
!     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 f11dgfe_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

!     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

!     Set size of real workspace
lwork = 2*n*m + 8*n + m*(m+2) + 100
Allocate (work(lwork))

!     Calculate factorization
ifail = 0
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)

!     Solve Ax = b using F11DGF
ifail = 0
Call f11dgf(method,n,nnz,a,la,irow,icol,nb,istb,indb,lindb,ipivp,ipivq,  &
istr,idiag,b,m,tol,maxitn,x,rnorm,itn,work,lwork,ifail)

Write (nout,99999) itn
Write (nout,99998) rnorm
Write (nout,*)

!     Output x
Write (nout,*) ' Solution vector  X'
Write (nout,*) ' ------------------'
Write (nout,99997) x(1:n)

100   Continue

99999 Format (' Converged in',I10,' iterations')
99998 Format (' Final residual norm =',1P,D16.3)
99997 Format (F8.4)

Contains

Subroutine f11dgfe_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, ind, iover, k, l, m, nadd, row
!       .. Executable Statements ..

!       Find the number of nonzero elements in each row of the matrix A, and
!       the 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 nonzero elements in row ROW.
Do i = iwork(n+row), iwork(n+row+1) - 1

!               If the column index of the nonzero 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.
iwork(1) = -999
Exit blocks
End If

Do i = istb(nb+1) - 1, istb(k+1), -1
End Do

Do i = 1, nadd
l = istb(k+1) + i - 1
indb(l) = iwork(2*n+1+i)
End Do

Do m = k + 1, nb + 1
istb(m) = istb(m) + nadd
End Do
End Do
End Do blocks

Return

End Subroutine f11dgfe_overlap
End Program f11dgfe
```