! F11DTF Example Program Text
! Mark 28.6 Release. NAG Copyright 2022.
Program f11dtfe
! .. Use Statements ..
Use nag_library, Only: f11dtf, 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 ..
Complex (Kind=nag_wp), Allocatable :: a(:), b(:)
Real (Kind=nag_wp), Allocatable :: 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 ..
! Print example header
Write (nout,*) 'F11DTF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
! Get the square matrix size
Read (nin,*) n
! Allocate arrays with lengths based on mesh.
liwork = 9*n + 3
Allocate (b(n),iwork(liwork))
! Get the number of non zero (nnz) matrix entries
Read (nin,*) nnz
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 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 f11dtfe_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,99998) ' N =', n
Write (nout,99998) ' NNZ =', nnz
Write (nout,99998) ' NB =', nb
Do k = 1, nb
Write (nout,99994) ' 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 = 0
Call f11dtf(n,nnz,a,la,irow,icol,nb,istb,indb,lindb,lfill,dtol,pstrat, &
milu,ipivp,ipivq,istr,idiag,nnzc,npivm,iwork,liwork,ifail)
! Output details of the factorization
Write (nout,99997) ' Factorization'
Write (nout,99998) ' NNZC =', nnzc
Write (nout,99997) ' Elements of factorization'
Write (nout,99997) &
' I J C(I,J) Index'
Do k = 1, nb
Write (nout,99995) ' C_', k, &
' -----------------------------------------------'
! Elements of the k-th block
Do i = istr(istb(k)), istr(istb(k+1)) - 1
Write (nout,99993) irow(i), icol(i), a(i), i
End Do
End Do
Write (nout,99997) ' Details of factorized blocks'
If (maxval(npivm(1:nb))>0) Then
! Including pivoting details.
Write (nout,99997) &
' K I ISTR(I) IDIAG(I) INDB(I) IPIVP(I) IPIVQ(I)'
Do k = 1, nb
i = istb(k)
Write (nout,99996) k, i, istr(i), idiag(i), indb(i), ipivp(i), &
ipivq(i)
Do i = istb(k) + 1, istb(k+1) - 1
Write (nout,99999) 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,99997) ' K I ISTR(I) IDIAG(I) INDB(I)'
Do k = 1, nb
i = istb(k)
Write (nout,99996) k, i, istr(i), idiag(i), indb(i)
Do i = istb(k) + 1, istb(k+1) - 1
Write (nout,99999) i, istr(i), idiag(i), indb(i)
End Do
Write (nout,*) ' ------------------------------------'
End Do
End If
100 Continue
99999 Format (1X,I7,5(I10))
99998 Format (1X,A,I4)
99997 Format (1X,/,1X,A)
99996 Format (1X,I3,1X,I3,5(I10))
99995 Format (1X,A3,I1,A)
99994 Format (1X,A,I3,A,I3,A,I3)
99993 Format (4X,2I4,4X,'(',E13.5,',',E13.5,')',I8)
Contains
Subroutine f11dtfe_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 F11DTF 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 F11DTF 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 ..
! Find the number of nonzero 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 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.
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 f11dtfe_overlap
End Program f11dtfe