!   F02FJF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    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
!dec$ attributes offload : mic :: dot
       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.

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: dot
       END FUNCTION dot
!dec$ attributes offload : mic :: image
       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
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: image
       END SUBROUTINE image
!dec$ attributes offload : mic :: monit
       SUBROUTINE monit(istate,nextit,nevals,nevecs,k,f,d)
!         Monitoring routine for F02FJF.

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)
!dec$ attributes offload : mic :: monit
       END SUBROUTINE monit
    END MODULE f02fjfe_mod
    PROGRAM f02fjfe

!      F02FJF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : f02fjf, f02fjz, f11jaf, x04cbf
       USE f02fjfe_mod, ONLY : dot, image, nag_wp, 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
