      PROGRAM ARPACK_SYMMETRIC_EX1
      USE ARPACK_INT
      USE UMACH_INT
      USE WRRRN_INT
      IMPLICIT NONE
!     Compute the smallest eigenvalues of a discrete Laplacian,
!     based on second order divided differences.
!     The matrix used is the 2 dimensional discrete Laplacian on 
!     the unit square with zero Dirichlet boundary condition.                    
      INTEGER                  :: J, NOUT
      INTEGER, parameter       :: NEV=5   !number of Eigenvalues required
      INTEGER, parameter       :: NV=0, NX=10 
      INTEGER, parameter       :: N=NX**2 !size of matrix problem
      REAL(DKIND)              :: VALUES(NEV), RES(N), EF(NX, NX)
      REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)
      REAL(DKIND)              :: NORM
      LOGICAL                  :: SMALL, SOLVED
      TYPE(ARPACKBASE)         :: Q
      PROCEDURE(DMV)           :: FCN

      CALL UMACH(2, NOUT)
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program.  That happens within the
! routine ARPACK_SYMMETRIC().  It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_SYMMETRIC.  
      CALL ARPACK_SYMMETRIC(N, FCN, VALUES, &
           PLACE=ARPACK_Smallest_Magnitude,  VECTORS=VECTORS)
      WRITE(NOUT, *) 'Number of eigenvalues requested, and declared accurate'
      WRITE(NOUT, *) '------------------------------------------------------'
      WRITE(NOUT, '(5X, I4, 5X, I4)') NEV, Q%NACC
      WRITE(NOUT, *) 'Number of Matrix-Vector Products Recorded, EX-11'
      WRITE(NOUT, *) '------------------------------------------------'
      WRITE(NOUT, '(5X, I4)') NV
      CALL WRRRN('Smallest Laplacian Eigenvalues', VALUES)

! Check residuals, A*vectors = values*vectors:
      DO J=1,NEV
! Necessary to have an unused TYPE(ARPACKBASE) :: Q as an argument:         
          RES=FCN(VECTORS(:,J), ARPACK_A_x, Q)-VALUES(J)*VECTORS(:,J)
          NORM=maxval(abs(RES)) 
          SMALL=(NORM <= ABS(VALUES(J))*SQRT(EPSILON(NORM)))
          IF(J==1) SOLVED=SMALL
          SOLVED=SOLVED .and. SMALL
      END DO
      IF(SOLVED) THEN
         WRITE(nout,'(A///)') &
              'All Ritz Values and Vectors have small residuals.'
      ELSE
         WRITE(nout,'(A///)') &
               'Some Ritz Values and Vectors have large residuals.'
      ENDIF

! The first eigenvector is scaled to be positive.
! It defines the eigenfunction for the smallest
! eigenvalue at the grid defined by the differencing.         
      EF=reshape(VECTORS(:,1),(/NX,NX/))
      CALL WRRRN('First 2D Laplacian Eigenfunction at Grid Points', EF)
    
      END
      
      FUNCTION FCN(X, TASK, EX)RESULT(Y)
         USE ARPACK_INT

         CLASS(ARPACKBASE),INTENT(INOUT) :: EX
         REAL(DKIND), INTENT(INOUT) :: X(:)
         INTEGER, INTENT(IN) :: TASK
         REAL(DKIND) Y(SIZE(X))
        
! Local variables:                
         INTEGER J
         INTEGER, SAVE :: NX
         REAL(DKIND), SAVE :: HSQ

         SELECT CASE(TASK)
            CASE(ARPACK_A_x)
!     Computes y <-- A*x, where A is the N**2 by N**2 block        
!     tridiagonal matrix                                                
!                                                                       
!                  | T -I          |                                    
!                  |-I  T -I       |                                    
!              A = |   -I  T       |                                    
!                  |        ...  -I|                                    
!                  |           -I T|            
               Y(1:NX)=T(NX,X(1:NX)) - X(NX+1:2*NX)
               DO J=NX+1,NX**2-NX,NX
                  Y(J:J+NX-1)=T(NX,X(J:J+NX-1))&
                           - X(J-NX:J-1)-X(J+NX:J+2*NX-1)		   
               END DO 
               Y((NX-1)*NX+1:NX**2)= - X((NX-1)*NX-NX+1:(NX-1)*NX)&
                                  + T(NX,X((NX-1)*NX+1:NX**2))
! Note that HSQ is passed as a component of the extended type.		                                   
               Y=(1._DKIND/HSQ)*Y
            
            CASE(ARPACK_Prepare)
! Define NX, 1/H**2 so they are later available in the evaluator. 
               NX=10 ! This value is fixed in the evaluator.           
               HSQ = 1._DKIND/REAL(NX+1,DKIND)**2
               Y=0._DKIND
            CASE DEFAULT
               WRITE(NOUT,*) TASK, ': INVALID TASK REQUESTED'
               STOP 'IMSL_ERROR_WRONG_OPERATION'
            END SELECT
      CONTAINS
         FUNCTION T(NX, X)RESULT(V)
            INTEGER, INTENT(IN) :: NX
            REAL(DKIND), INTENT(IN) :: X(:)
            REAL(DKIND) :: V(NX)
            REAL(DKIND) :: MONE=-1._DKIND, FOUR=4._DKIND
            INTEGER J
         
            V(1)=FOUR*X(1)+MONE*X(2)
            DO J=2,NX-1
               V(J)=MONE*X(J-1)+FOUR*X(J)+MONE*X(J+1)
            END DO
            V(NX)=MONE*X(NX-1)+FOUR*X(NX)
         END FUNCTION    
      END FUNCTION
