! F11DTF Example Program ! ! NAG FORTRAN Library. ! Mark 24 Release. NAG Copyright 2012. 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 .. Continue ! 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 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 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 .. 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 f11dtfe_overlap End Program f11dtfe