!   F02FJF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module f02fjfe_mod

!     F02FJF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter        :: half = 0.5_nag_wp
      Real (Kind=nag_wp), Parameter        :: zero = 0.0_nag_wp
      Integer, Parameter                   :: nin = 5, nout = 6
    Contains
      Function dot(iflag,n,z,w,ruser,lruser,iuser,liuser)
!       This function implements the dot product - transpose(W)*B*Z.
!       DOT assumes that N is at least 3.

!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: dot
!       .. Scalar Arguments ..
        Integer, Intent (Inout)              :: iflag
        Integer, Intent (In)                 :: liuser, lruser, n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(lruser)
        Real (Kind=nag_wp), Intent (In)      :: w(n), z(n)
        Integer, Intent (Inout)              :: iuser(liuser)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: s
        Integer                              :: i
!       .. Executable Statements ..
        s = zero
        s = s + (z(1)-half*z(2))*w(1)
        s = s + (-half*z(n-1)+z(n))*w(n)
        Do i = 2, n - 1
          s = s + (-half*z(i-1)+z(i)-half*z(i+1))*w(i)
        End Do
        dot = s
        Return
      End Function dot
      Subroutine image(iflag,n,z,w,ruser,lruser,iuser,liuser)
!       This routine solves  A*W = B*Z  for W.
!       The routine assumes that N is at least 3.

!       The data A, NNZ, LA, IROW, ICOL, IPIV and ISTR on exit from
!       F11JAF have been packed into the xUSER communication arrays in
!       the following way:
!       IUSER(1:2) = (/NNZ, LA/)
!       RUSER(1:LA) = A
!       IUSER(3:(2*LA+2*N+3)) = (/IROW, ICOL, IPIV, ISTR/)
!       We'll also use RUSER((LA+1):(LA+N)) as space for F11JCF's dummy
!       arg. B, and RUSER((LA+N+1):(LA+7*N+120)) as space for F11JCF's
!       dummy arg. WORK

!       .. Use Statements ..
        Use nag_library, Only: f11jcf, x02ajf
!       .. Scalar Arguments ..
        Integer, Intent (Inout)              :: iflag
        Integer, Intent (In)                 :: liuser, lruser, n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(lruser)
        Real (Kind=nag_wp), Intent (Out)     :: w(n)
        Real (Kind=nag_wp), Intent (In)      :: z(n)
        Integer, Intent (Inout)              :: iuser(liuser)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: rnorm, tol
        Integer                              :: ifail, itn, j, la, lwork,      &
                                                maxitn, nnz
        Character (2)                        :: method
!       .. Executable Statements ..
        nnz = iuser(1)
        la = iuser(2)

!       Form B*Z in RUSER((LA+1):(LA+N)) and initialize W to
!       zero.
        w(1:n) = zero
        ruser(la+1) = z(1) - half*z(2)
        Do j = 2, n - 1
          ruser(la+j) = -half*z(j-1) + z(j) - half*z(j+1)
        End Do
        ruser(la+n) = -half*z(n-1) + z(n)

!       Call F11JCF to solve the equations  A*W = B*Z.
        method = 'CG'
        tol = x02ajf()
        maxitn = 100
        lwork = 6*n + 120

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 1
        Call f11jcf(method,n,nnz,ruser,la,iuser(3),iuser(la+3),iuser(2*la+3), &
          iuser(2*la+n+3),ruser(la+1),tol,maxitn,w,rnorm,itn,ruser(la+n+1), &
          lwork,ifail)

        If (ifail>0) iflag = -ifail
        Return
      End Subroutine image
      Subroutine monit(istate,nextit,nevals,nevecs,k,f,d)
!       Monitoring routine for F02FJF.

!       .. Parameters ..
        Integer, Parameter                   :: nout = 6
!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: istate, k, nevals, nevecs,     &
                                                nextit
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: d(k), f(k)
!       .. Local Scalars ..
        Integer                              :: i
!       .. Executable Statements ..
        If (istate/=0) Then
          Write (nout,*)
          Write (nout,99999) '  ISTATE = ', istate, ' NEXTIT = ', nextit
          Write (nout,99999) '  NEVALS = ', nevals, ' NEVECS = ', nevecs
          Write (nout,*) '       F           D'
          Write (nout,99998)(f(i),d(i),i=1,k)
        End If
        Return

99999   Format (1X,A,I4,A,I4)
99998   Format (1X,1P,E11.3,3X,E11.3)
      End Subroutine monit
    End Module f02fjfe_mod
    Program f02fjfe

!     F02FJF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: f02fjf, f02fjz, f11jaf, nag_wp, x04cbf
      Use f02fjfe_mod, Only: dot, image, nin, nout, zero
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: dscale, dtol, tol
      Integer                              :: i, ifail, k, l, la, ldx, lfill,  &
                                              liuser, lruser, lwork, m, n,     &
                                              nnz, nnzc, noits, novecs, npivm
      Character (1)                        :: mic, pstrat
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: a(:), d(:), ruser(:), work(:),   &
                                              x(:,:)
      Integer, Allocatable                 :: icol(:), ipiv(:), irow(:),       &
                                              istr(:), iuser(:)
      Character (1)                        :: clabs(1), rlabs(1)
!     .. Executable Statements ..
      Write (nout,*) 'F02FJF Example Program Results'
      Write (nout,*)
      Flush (nout)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, m, k, tol
      la = 10*n
      ldx = n
      liuser = 2*la + 2*n + 3
      lruser = la + 7*n + 120
      lwork = 5*k + 2*n
      Allocate (a(la),d(n),ruser(lruser),work(lwork),x(ldx,k),icol(la), &
        ipiv(n),irow(la),istr(n+1),iuser(liuser))

!     Set up the sparse symmetric coefficient matrix A.
      l = 0
      Do i = 1, n
        If (i>=5) Then
          l = l + 1
          a(l) = -0.25_nag_wp
          irow(l) = i
          icol(l) = i - 4
        End If
        If (i>=2) Then
          l = l + 1
          a(l) = -0.25_nag_wp
          irow(l) = i
          icol(l) = i - 1
        End If
        l = l + 1
        a(l) = 1.0_nag_wp
        irow(l) = i
        icol(l) = i
      End Do
      nnz = l

!     Call F11JAF to find an incomplete Cholesky factorisation of A.
      lfill = 2
      dtol = zero
      mic = 'M'
      dscale = zero
      pstrat = 'M'

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call f11jaf(n,nnz,a,la,irow,icol,lfill,dtol,mic,dscale,pstrat,ipiv,istr, &
        nnzc,npivm,iuser,liuser,ifail)

!     Call F02FJF to find eigenvalues and eigenvectors.

      noits = 1000
      novecs = 0

!     Communicate A, NNZ, LA, IROW, ICOL, IPIV and ISTR to IMAGE 
!     thread-safely using RUSER and IUSER.
!     In addition to using RUSER for storing A, we'll also use
!     7*N+120 elements of RUSER in place of local arrays in IMAGE.

!     Initialized A goes into ruser.
      ruser(1:nnz+nnzc) = a(1:nnz+nnzc)

!     NNZ, LA, IROW, ICOL, IPIV and ISTR go into IUSER, in that order.
!     Only the first NNZ+NNZC elements of IROW and ICOL will have been
!     initialized:
      iuser(1) = nnz
      iuser(2) = la
      iuser(3:2+nnz+nnzc) = irow(1:nnz+nnzc)
      iuser(la+3:la+2+nnz+nnzc) = icol(1:nnz+nnzc)
      iuser(2*la+3:2*la+2+n) = ipiv(1:n)
      iuser(liuser-n:liuser) = istr(1:n+1)

!     * To obtain monitoring information from the supplied subroutine MONIT,
!     replace the name F02FJZ by MONIT in the call to F02FJF, and USE MONIT
!     from f02fjfe_mod *
      ifail = -1
      Call f02fjf(n,m,k,noits,tol,dot,image,f02fjz,novecs,x,ldx,d,work,lwork, &
        ruser,lruser,iuser,liuser,ifail)

      If (ifail>=0) Then
        If (ifail/=1 .And. ifail<=4 .And. m>=1) Then
          Do i = 1, m
            d(i) = 1.0_nag_wp/d(i)
          End Do
          Write (nout,*) 'Final results'
          Write (nout,*)
          Write (nout,*) '  Eigenvalues'
          Write (nout,99999) d(1:m)
          Write (nout,*)
          Flush (nout)
!         Normalize eigenvectors
          Do i = 1, m
            x(1:n,i) = x(1:n,i)/x(1,i)
          End Do
          Call x04cbf('General',' ',n,m,x,ldx,'1P,E12.3','  Eigenvectors','N', &
            rlabs,'N',clabs,80,0,ifail)
        End If
      End If

99999 Format (1X,1P,4E12.3)
    End Program f02fjfe
