      MODULE ARPACK_COMPLEX_EX1_INT
      USE ARPACK_INT
      IMPLICIT NONE

      TYPE, EXTENDS(ARPACKBASE), PUBLIC :: ARPACKBASE_EXT
         REAL(DKIND) :: H  =0._DKIND
         REAL(DKIND) :: HSQ=0._DKIND
         COMPLEX(DKIND) :: RHO=(0._DKIND,0._DKIND)
         COMPLEX(DKIND) :: DL=(0._DKIND,0._DKIND)
         COMPLEX(DKIND) :: DD=(0._DKIND,0._DKIND)
         COMPLEX(DKIND) :: DU=(0._DKIND,0._DKIND)
         INTEGER :: NX=0
         INTEGER :: NV=0
      END TYPE ARPACKBASE_EXT 
      CONTAINS
      FUNCTION FZ1(X, TASK, EXTYPE) RESULT(Y)
         USE UMACH_INT
         
         CLASS (ARPACKBASE), INTENT(INOUT) :: EXTYPE
         COMPLEX (DKIND), INTENT(INOUT) :: X(:)	
         INTEGER, INTENT(IN) :: TASK
         COMPLEX (DKIND) Y(size(X))
         COMPLEX (DKIND) DT(3)
         REAL(DKIND) :: ONE=1._DKIND
         INTEGER J, NOUT

         CALL UMACH(2, NOUT)
         SELECT TYPE(EXTYPE)
            TYPE IS(ARPACKBASE_EXT)
            ASSOCIATE(NX  => EXTYPE % NX,&
                      H   => EXTYPE % H,&
                      HSQ => EXTYPE % HSQ,&
                      RHO => EXTYPE % RHO,&
                      DL  => EXTYPE % DL,&
                      DD  => EXTYPE % DD,&
                      DU  => EXTYPE % DU,&
                      NV  => EXTYPE % NV)
                      
            SELECT CASE(TASK)
               CASE(ARPACK_A_x)
!     Computes y <-- A*x, where A is the N**2 by N**2 block        
!     tridiagonal matrix deriving from (Laplacian u) + rho*(du/dx).                                                    

                  DT=(/DL,DD,DU/)            
                  Y(1:NX)=T(NX,X(1:NX),DT) - X(NX+1:2*NX)/HSQ
                  DO J=NX+1,NX**2-NX,NX
                     Y(J:J+NX-1)=T(NX,X(J:J+NX-1),DT) &
                                 - (X(J-NX:J-1)+ X(J+NX:J+2*NX-1))/HSQ		   
                  END DO 
                  Y((NX-1)*NX+1:NX**2)= - X((NX-1)*NX-NX+1:(NX-1)*NX) &
                                  / HSQ + T(NX,X((NX-1)*NX+1:NX**2),DT)
! Total the number of matrix-vector products.		    
                  NV=NV+1
            
               CASE(ARPACK_Prepare)
! Define 1/H**2, etc. so they are available in the evaluator.
                  H   = ONE/REAL(NX+1,DKIND)          
                  HSQ = H**2
                  DD  = (4.0D+0, 0.0D+0)  / HSQ
                  DL  = -ONE/HSQ - (5.0D-1, 0.0D+0) *RHO/H 
                  DU  = -ONE/HSQ + (5.0D-1, 0.0D+0) *RHO/H                
                  NV = 0    	       
               CASE DEFAULT
                  WRITE(nout,*) TASK, ': INVALID TASK REQUESTED'
                  STOP 'IMSL_ERROR_WRONG_OPERATION'
            END SELECT
            END ASSOCIATE
         END SELECT
      END FUNCTION

      FUNCTION T(NX, X, DT)RESULT(V)
         INTEGER, INTENT(IN) :: NX
         COMPLEX(DKIND), INTENT(IN) :: X(:), DT(3)
         COMPLEX(DKIND) :: V(NX)
         INTEGER J
         ASSOCIATE(DL => DT(1),&
                   DD => DT(2),&
                   DU => DT(3))
         V(1) =  DD*X(1) + DU*X(2) 
         DO J = 2,NX-1 
            V(J) = DL*X(J-1) + DD*X(J) + DU*X(J+1) 
         END DO
         V(NX) =  DL*X(NX-1) + DD*X(NX)
         END ASSOCIATE
      END FUNCTION    
      END MODULE

! Compute the largest magnitude eigenvalues of a discrete Laplacian,
! based on second order divided differences.

!     The matrix used is obtained from the standard central difference      
!         discretization of the convection-diffusion operator           
!                 (Laplacian u) + rho*(du / dx)                         
!         on the unit squre 0,1x0,1 with zero Dirichlet boundary    
!         conditions.                                                                                   
      USE ARPACK_COMPLEX_EX1_INT
      USE UMACH_INT
      USE WRCRN_INT
  
      INTEGER, PARAMETER :: NEV=6
      INTEGER ::  J, N, NOUT
      COMPLEX(DKIND) :: VALUES(NEV)
      COMPLEX(DKIND), ALLOCATABLE :: RES(:), EF(:,:)
      COMPLEX(DKIND), ALLOCATABLE :: VECTORS(:,:)
      REAL(DKIND) NORM
      LOGICAL SMALL, SOLVED
      TYPE(ARPACKBASE_EXT) EX

      ASSOCIATE(NX   => EX % NX, &
                NV   => EX % NV, &
                RHO  => EX % RHO,&
                NACC => EX % NACC)
              
      CALL UMACH(2, NOUT)        
      NX=10
      RHO=(100._DKIND,1._DKIND)

! Define size of matrix problem.    
      N=NX**2
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program.  That happens within the
! routine ARPACK_COMPLEX().  It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_COMPLEX.  
      CALL ARPACK_COMPLEX(N, FZ1, VALUES, EXTYPE=EX, VECTORS=VECTORS)
      
      WRITE(NOUT, *) 'Number of eigenvalues requested, and accurate'
      WRITE(NOUT, *) '---------------------------------------------'
      WRITE(NOUT, '(5X, I4, 5X, I4)') NEV, NACC
      WRITE(NOUT, *) 'Number of Matrix-Vector Products Required, ZEX-1'
      WRITE(NOUT, *) '------------------------------------------------'
      WRITE(NOUT,  '(5X, I4)') NV
      CALL WRCRN ('Largest Magnitude Operator Eigenvalues', VALUES)
! Check residuals, A*vectors = values*vectors:
      ALLOCATE(RES(N))
      DO J=1,NACC
         RES=FZ1(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.'
      END IF

      END ASSOCIATE
      END 
