
      MODULE ARPACK_NONSYMMETRIC_EX1_INT
      USE ARPACK_INT
      USE LSLCQ_INT
      USE N1RTY_INT
      IMPLICIT NONE

      TYPE, EXTENDS(ARPACKBASE) :: ARPACKBASE_EXT
         INTEGER ::     NX=0
         INTEGER ::     NV=0
! This example extends the base type to 
! information for solving complex tridiagonal systems.
         COMPLEX(DKIND), ALLOCATABLE :: A(:), B(:), C(:)
         INTEGER, ALLOCATABLE ::  IR(:), IS(:) 
! This controls the type of shifting.  When
! the value is 1, use real part of inv(A-*M)*x.
! If value is 2, use imaginary part of same.
         INTEGER :: SHIFT_STRATEGY=1
      END TYPE ARPACKBASE_EXT 
      
      CONTAINS
      FUNCTION FCN(X, TASK, EX) RESULT(Y)
         USE UMACH_INT
         
         CLASS (ARPACKBASE), INTENT(INOUT) :: EX
         REAL (DKIND), INTENT(INOUT) :: X(:)
         INTEGER, INTENT(IN) :: TASK
         INTEGER, PARAMETER :: NSIZE=12
         REAL(DKIND) Y(size(X))
         REAL(DKIND) DL, DD, DU
         COMPLEX(DKIND) Z(2*size(X))
         REAL(DKIND) U(2*size(X))
         INTEGER J, IERR, NOUT, IJOB
         CALL UMACH(2, NOUT)
         SELECT TYPE(EX)
            TYPE IS(ARPACKBASE_EXT)
            ASSOCIATE(N   =>    EX % NCOLS,&
                      NV  =>    EX % NV, &
                      SHIFT =>  EX % SHIFT)

            SELECT CASE(TASK)

               CASE(ARPACK_A_x)
                  DL =  -2._DKIND
                  DD =   2._DKIND 
                  DU =   3._DKIND
                  Y(1) =  DD*X(1) + DU*X(2) 
                  DO J = 2,N-1 
                     Y(J) = DL*X(J-1) + DD*X(J) + DU*X(J+1)
                  END DO 
                  Y(N) =  DL*X(N-1) + DD*X(N)
                  NV=NV+1 
                            
               CASE(ARPACK_B_x)
                  Y(1) = 4._DKIND*X(1) + X(2)
                  DO J = 2,N-1 
                     Y(J) = X(J-1) + 4._DKIND*X(J) + X(J+1)
                  END DO   
                  Y(N) = X(N-1) + 4._DKIND*X(N) 
                  NV=NV+1
                  
               CASE(ARPACK_inv_A_minus_Shift_B_x)
! Compute Y=REAL/AIMAG(inv(A-*B)*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.
                  Z=CMPLX(X,0._DKIND,DKIND)
                  IJOB = 2
                  CALL LSLCQ (EX%C, EX%A, EX%B, Z, U, &
                              EX%IR, EX%IS, N=N, IJOB=IJOB)
                  IERR= N1RTY(1)
                  IF (IERR==4 .OR. IERR==5)  &
                      STOP 'IMSl_FATAL_ERROR_SOLVING'


                  IF(EX % SHIFT_STRATEGY == 1) THEN
                     Y(1:N)=REAL( Z(1:N),DKIND)
                  ELSE IF (EX % SHIFT_STRATEGY == 2)THEN
                     Y(1:N)=AIMAG(Z(1:N))
                  END IF

! 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%IR)) DEALLOCATE(EX%IR)
                  IF (ALLOCATED(EX%IS)) DEALLOCATE(EX%IS)
                  ALLOCATE(EX%A(2*N), EX%B(2*N), EX%C(2*N), &
                           EX%IR(NSIZE), EX%IS(NSIZE), STAT=IERR)
                  IF (IERR /= 0) STOP 'IMSL_ERROR_ALLOCATING_WORKSPACE'
 
! Define matrix, A-SHIFT*B.
                  EX % B(1:N) =  -2._DKIND-SHIFT
                  EX % A(1:N) =  2._DKIND-4._DKIND*SHIFT
                  EX % C(1:N) =  3._DKIND-SHIFT 
! Factor the matrix with LU and partial pivoting.
                  IJOB = 3
                  CALL LSLCQ (EX%C, EX%A, EX%B, Z, 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 ZVALUES to satisfy compiler.
                  Y=0._DKIND
                  NV=0
                  
               CASE DEFAULT
                  WRITE(NOUT,*) TASK, ': INVALID OPERATION REQUESTED'
                  STOP 'IMSL_ERROR_WRONG_OPERATION'
               END SELECT
            END ASSOCIATE 
            END SELECT
      END FUNCTION
      END MODULE

      
!         Suppose we want to solve A*x = lambda*B*x in -invert mode
!         The matrix A is the tridiagonal matrix with 2 on the diagonal,
!         -2 on the subdiagonal and 3 on the superdiagonal.  The matrix 
!         is the tridiagonal matrix with 4 on the diagonal and 1 on the 
!         off-diagonals.
!         The  sigma is a complex number (sigmar, sigmai).
!         OP = Real_Part{invA-(SIGMAR,SIGMAI)*B*B.  
      
      USE ARPACK_NONSYMMETRIC_EX1_INT
      USE UMACH_INT
      USE WRCRN_INT
      
      INTEGER, PARAMETER ::  NEV=6, N=100
      COMPLEX(DKIND)  :: ZVALUES(NEV), RES(N),U(N),V(N),W(N)
      REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)
      REAL(DKIND) NORM
      LOGICAL SKIP, SMALL, SOLVED
      INTEGER J, STRATEGY, NOUT
      CHARACTER(LEN=12) TAG    
      CHARACTER(LEN=60) TITLE
      TYPE(ARPACKBASE_EXT) EX
    
      ASSOCIATE(NX  => EX % NX, &
                NV  => EX % NV, &
                SHIFT  => EX % SHIFT,&
                FACTOR => EX % FACTOR_MAXNCV,& 
                NACC   => EX % NACC)
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program.  That happens within the
! routine ARPACK_NONSYMMETRIC().  It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_NONSYMMETRIC.

      CALL UMACH(2, NOUT)
      SOLVED=.TRUE. 
      DO STRATEGY=1,2
         SHIFT=CMPLX(0.4_DKIND,0.6_DKIND,DKIND)
         FACTOR=5._DKIND
         EX % SHIFT_STRATEGY=STRATEGY
    
         CALL ARPACK_NONSYMMETRIC(N, FCN, ZVALUES,         &
                 TYPE=ARPACK_Generalized,                   &
                 CATEGORY=ARPACK_Complex_Part_Shift_Invert, &
                 EXTYPE=EX, VECTORS=VECTORS)
                 
         WRITE(NOUT, *) &
             'Number of Matrix-Vector Products Required, NS Ex-1'
         WRITE(NOUT, *) &
             '--------------------------------------------------'
         WRITE(NOUT, '(5X,I4)') NV
         WRITE(NOUT, *) 'Number of accurate values determined'
         WRITE(NOUT, *) '------------------------------------'
         WRITE(NOUT, '(5X, I4)') NACC
! Check residuals, A*vectors = ZVALUES*M*vectors:

         SKIP=.FALSE.         
         DO J=1,NACC
            IF(SKIP) THEN
               SKIP=.FALSE.
               CYCLE
            END IF
! The eigenvalue is complex and the pair of vectors
! for the complex eigenvector is returned.         
            IF(AIMAG(ZVALUES(J)) /= 0._DKIND)THEN
! Make calls for real and imaginary parts of eigenvectors
! applied to the operators A, M                     
               U=CMPLX(FCN(VECTORS(:,J),ARPACK_A_x,EX),&
                       FCN(VECTORS(:,J+1),ARPACK_A_x,EX),DKIND)
               V=CMPLX(FCN(VECTORS(:,J),ARPACK_B_x,EX),&
                       FCN(VECTORS(:,J+1),ARPACK_B_x,EX),DKIND)
! Since the matrix is real, there is an additional conjugate:
               RES=U-ZVALUES(J)*V
               SKIP=.TRUE.
            ELSE
! The eigenvalue is real and the real eigenvector is returned.
        
               RES=FCN(VECTORS(:,J),ARPACK_A_x,EX)-ZVALUES(J)*&
                   FCN(VECTORS(:,J),ARPACK_B_x,EX)
            END IF
            NORM=maxval(abs(RES))
            SMALL=(NORM <= ABS(ZVALUES(J))*SQRT(EPSILON(NORM)))
            SOLVED=SOLVED .and. SMALL
         END DO
         IF(STRATEGY==1) TAG='REAL SHIFT'
         IF(STRATEGY==2) TAG='IMAG SHIFT'
  
         TITLE = 'Largest Raleigh Quotient Eigenvalues,'//TAG
         CALL WRCRN(TITLE, ZVALUES)

         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 DO ! Shift strategy
      END ASSOCIATE  

      END 
