NAG Library Manual, Mark 29
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
!   F02EKF Example Program Text
!   Mark 29.0 Release. NAG Copyright 2023.
    Module f02ekfe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: mymonit, myoption
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine myoption(icomm,comm,istat,iuser,ruser)

!       .. Use Statements ..
        Use nag_library, Only: f12adf
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (Inout)        :: istat
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: comm(*), ruser(*)
        Integer, Intent (Inout)        :: icomm(*), iuser(*)
!       .. Local Scalars ..
        Integer                        :: ifail1
        Character (25)                 :: rec
!       .. Intrinsic Procedures ..
        Intrinsic                      :: max
!       .. Executable Statements ..

        istat = 0

        If (iuser(1)>0) Then
          Write (rec,99999) 'Print Level=', iuser(1)
          ifail1 = 1
          Call f12adf(rec,icomm,comm,ifail1)
          istat = max(istat,ifail1)
        End If
        If (iuser(2)>100) Then
          Write (rec,99999) 'Iteration Limit=', iuser(2)
          ifail1 = 1
          Call f12adf(rec,icomm,comm,ifail1)
          istat = max(istat,ifail1)
        End If
        If (iuser(3)>0) Then
          ifail1 = 1
          Call f12adf('Shifted Inverse Real',icomm,comm,ifail1)
          istat = max(istat,ifail1)
        End If
99999   Format (A,I5)
      End Subroutine myoption
      Subroutine mymonit(ncv,niter,nconv,w,rzest,istat,iuser,ruser)

!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (Inout)        :: istat
        Integer, Intent (In)           :: nconv, ncv, niter
!       .. Array Arguments ..
        Complex (Kind=nag_wp), Intent (In) :: w(ncv)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: rzest(ncv)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Executable Statements ..

        If (iuser(4)>0) Then
          If (niter==1 .And. iuser(3)>0) Then
            Write (nout,99999) ' Arnoldi basis vectors used:', ncv
            Write (nout,*)                                                     &
              ' The following Ritz values (mu) are related to the'
            Write (nout,*)                                                     &
              ' true eigenvalues (lambda) by lambda = sigma + 1/mu'
          End If
          Write (nout,*)
          Write (nout,99999) ' Iteration number ', niter
          Write (nout,99998) ' Ritz values converged so far (', nconv,         &
            ') and their Ritz estimates:'
          Do i = 1, nconv
            Write (nout,99997) i, w(i), rzest(i)
          End Do
          Write (nout,*) ' Next (unconverged) Ritz value:'
          Write (nout,99996) nconv + 1, w(nconv+1)
        End If
        istat = 0
99999   Format (1X,A,I4)
99998   Format (1X,A,I4,A)
99997   Format (1X,1X,I4,1X,'(',E13.5,',',E13.5,')',1X,E13.5)
99996   Format (1X,1X,I4,1X,'(',E13.5,',',E13.5,')')
      End Subroutine mymonit
    End Module f02ekfe_mod
    Program f02ekfe

!     Example problem for F02EKF.

!     .. Use Statements ..
      Use f02ekfe_mod, Only: mymonit, myoption, nin, nout
      Use nag_library, Only: f02ekf, nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: three = 3.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: h, rho, s, sigma
      Integer                          :: i, ifail, imon, k, ldv, maxit, mode, &
                                          n, nconv, ncv, nev, nnz, nx, prtlvl
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: w(:)
      Real (Kind=nag_wp), Allocatable  :: a(:), resid(:), v(:,:)
      Real (Kind=nag_wp)               :: ruser(1)
      Integer, Allocatable             :: icolzp(:), irowix(:)
      Integer                          :: iuser(4)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'F02EKF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) nx
      Read (nin,*) nev
      Read (nin,*) ncv
      Read (nin,*) rho
      Read (nin,*) sigma

      n = nx*nx
      nnz = 3*n - 2
      ldv = n

      Allocate (resid(ncv),a(nnz),icolzp(n+1),irowix(nnz),w(ncv),v(ldv,ncv))

!     Construct  A in compressed column storage (CCS) format where:
!      A_{i,i}   = 2 + i
!      A_{i+1,i) = 3
!      A_{i,i+1} = rho/(2n+2) - 1

      h = one/real(n+1,kind=nag_wp)
      s = rho*h/two - one

      a(1) = two + one
      a(2) = three
      icolzp(1) = 1
      irowix(1) = 1
      irowix(2) = 2
      k = 3
      Do i = 2, n - 1
        icolzp(i) = k
        irowix(k) = i - 1
        irowix(k+1) = i
        irowix(k+2) = i + 1
        a(k) = s
        a(k+1) = two + real(i,kind=nag_wp)
        a(k+2) = three
        k = k + 3
      End Do
      icolzp(n) = k
      icolzp(n+1) = k + 2
      irowix(k) = n - 1
      irowix(k+1) = n
      a(k) = s
      a(k+1) = two + real(n,kind=nag_wp)

!     Set some options via iuser array and routine argument OPTION.
!     iuser(1) = print level, iuser(2) = iteration limit,
!     iuser(3)>0 means shifted-invert mode
!     iuser(4)>0 means print monitoring info

      Read (nin,*) prtlvl
      Read (nin,*) maxit
      Read (nin,*) mode
      Read (nin,*) imon

      If (prtlvl>0) Then
        imon = 0
      End If

      iuser(1) = prtlvl
      iuser(2) = maxit
      iuser(3) = mode
      iuser(4) = imon

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call f02ekf(n,nnz,a,icolzp,irowix,nev,ncv,sigma,mymonit,myoption,nconv,  &
        w,v,ldv,resid,iuser,ruser,ifail)

      Write (nout,99999) nconv, sigma
      Do i = 1, nconv
        If (resid(i)>real(100*n,kind=nag_wp)*x02ajf()) Then
          Write (nout,99998) i, w(i), resid(i)
        Else
          Write (nout,99998) i, w(i)
        End If
      End Do

99999 Format (1X,/,' The ',I4,' Ritz values of closest to ',E13.5,' are:',/)
99998 Format (1X,I8,5X,'( ',E13.5,' , ',E13.5,' )',5X,E13.5)
    End Program f02ekfe