! F11DGF Example Program ! ! NAG FORTRAN Library. ! Mark 24 Release. NAG Copyright 2012. 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 .. Continue ! Print example header Write (nout,*) 'F11DGF 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),x(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 in matrix A Read (nin,*)(a(i),irow(i),icol(i),i=1,nnz) ! Read in RHS Read (nin,*) b(1:n) ! Read algorithmic parameters Read (nin,*) method Read (nin,*) lfillg, dtolg Read (nin,*) pstrag Read (nin,*) milug 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 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 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 .. Continue ! Find the number of non-zero 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 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 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