      MODULE ARPACK_SYMMETRIC_EX2_INT
      USE ARPACK_INT
      USE LSLCR_INT
      USE N1RTY_INT
      IMPLICIT NONE

      TYPE, EXTENDS(ARPACKBASE) :: ARPACKBASE_EXT
         REAL(DKIND) :: HSQ=0._DKIND
         INTEGER ::     NX=0, NV=0
! This example extends the base type to 
! information for solving tridiagonal systems.       
         REAL(DKIND), ALLOCATABLE :: A(:), B(:), C(:)
         REAL(DKIND), ALLOCATABLE :: Y1(:), U(:)
         INTEGER, ALLOCATABLE :: IR(:), IS(:)
      END TYPE ARPACKBASE_EXT
      
      CONTAINS
      
      FUNCTION FCN(X, TASK, EX) RESULT(Y)
         CLASS (ARPACKBASE), INTENT(INOUT) :: EX
         REAL (DKIND), INTENT(INOUT) :: X(:)	
         INTEGER, INTENT(IN) :: TASK
         REAL (DKIND) Y(size(X))
         INTEGER J, IERR, IJOB, NSIZE

         SELECT TYPE(EX)
            TYPE IS(ARPACKBASE_EXT)
            ASSOCIATE(N     =>    EX%NCOLS, &
                      NV    =>    EX%NV, &
                      HSQ   =>    EX%HSQ, &
                      SHIFT =>  EX%SHIFT)             
            SELECT CASE(TASK)
               CASE(ARPACK_A_x)
                  Y(1) =  2._DKIND*X(1) - X(2)
                   
                  DO J = 2,N-1 
                     Y(J) = - X(J-1) + 2._DKIND*X(J) - X(J+1) 
                  END DO
                  
                  Y(N) = - X(N-1) + 2._DKIND*X(N) 
                  Y=Y/HSQ		    
                  CASE(ARPACK_inv_A_minus_Shift_x)
! Compute Y=inv(A-*I)*x.  This is done with a solve
! step, using the LU factorization.  Note that the data
! for the factorization is stored in the user's extended
! data type.

                     EX%Y1(1:N) = X
                     IJOB = 2
                     CALL LSLCR (EX%C, EX%A, EX%B, EX%Y1, EX%U, &
                                  EX%IR, EX%IS, N=N, IJOB=IJOB)
                     Y(1:N) = EX%Y1(1:N)
                     IERR= N1RTY(1)
                     IF (IERR==4 .OR. IERR==5)  STOP 'IMSl_FATAL_ERROR_SOLVING'
! Total number of solve steps.                
                     NV=NV+1
                  CASE(ARPACK_Prepare)
! Set up storage areas for factored tridiagonal matrix.
                     IF (ALLOCATED(EX%A)) DEALLOCATE(EX%A)
                     IF (ALLOCATED(EX%B)) DEALLOCATE(EX%B)
                     IF (ALLOCATED(EX%C)) DEALLOCATE(EX%C)
                     IF (ALLOCATED(EX%Y1)) DEALLOCATE(EX%Y1)
                     IF (ALLOCATED(EX%U)) DEALLOCATE(EX%U)
                     IF (ALLOCATED(EX%IR)) DEALLOCATE(EX%IR)
                     IF (ALLOCATED(EX%IS)) DEALLOCATE(EX%IS)
                     NSIZE = (log(dble(N))/log(2.0)) + 5
                     ALLOCATE(EX%A(2*N), EX%B(2*N), EX%C(2*N), EX%Y1(2*N), &
                               EX%U(2*N), EX%IR(NSIZE),               &
                               EX%IS(NSIZE),   STAT=IERR)
                     IF (IERR /= 0) STOP 'IMSL_ERROR_ALLOCATING_WORKSPACE' 
! Define matrix values.
                     HSQ=1._DKIND/REAL((N+1)**2,DKIND)
                     EX%B(1:N) = -1._DKIND/HSQ
                     EX%A(1:N) = 2._DKIND/HSQ - SHIFT
                     EX%C(1:N) = EX%B(1:N)
                     EX%Y1(:) = 0.0_DKIND
! Factor the matrix with LU and partial pivoting.
                     IJOB = 3
                     CALL LSLCR (EX%C, EX%A, EX%B, EX%Y1, EX%U, &
                                  EX%IR, EX%IS, N=N, IJOB=IJOB)
                     IERR = N1RTY(1)   
                     IF(IERR == 4 .or. IERR == 5) STOP 'IMSL FATAL ERROR'
! Give output some values to satisfy compiler.                
                     Y=0._DKIND
                     NV=0
                  CASE DEFAULT
                     STOP 'IMSL_ERROR_WRONG_OPERATION'
                  END SELECT
            END ASSOCIATE 
         END SELECT
      END FUNCTION    
      END MODULE

      USE ARPACK_SYMMETRIC_EX2_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 is the 1 dimensional discrete Laplacian on       
! the interval 0,1 with zero Dirichlet boundary condition.                      
      
      INTEGER, PARAMETER ::  NEV=4, N=100
      REAL(DKIND) :: VALUES(NEV), RES(N)
      REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)
      REAL(DKIND) NORM
      LOGICAL SMALL, SOLVED
      INTEGER J, NOUT
      TYPE(ARPACKBASE_EXT) EX
    
      ASSOCIATE(NX  => EX%NX, &
                NV  => EX%NV, &
              SIGMA => EX%SHIFT)
              
      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.
      SIGMA=0._DKIND
      CALL ARPACK_SYMMETRIC(N, FCN, VALUES,&
         CATEGORY=ARPACK_Shift_Invert, EXTYPE=EX, VECTORS=VECTORS)
      
      WRITE(NOUT,*) 'Number of Matrix-Vector Products Required, EX-2'
      WRITE(NOUT,*) '-----------------------------------------------'
      WRITE(NOUT, '(5X, I4)') NV
      CALL WRRRN('Largest Laplacian Eigenvalues Near Zero Shift', &
                 VALUES)
! Check residuals, A*vectors = values*vectors:
      DO J=1,NEV
         RES=FCN(VECTORS(:,J),ARPACK_A_x,EX)-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
      END ASSOCIATE              
      END
      
