      MODULE ARPACK_SVD_EX1_INT
      USE ARPACK_INT
      IMPLICIT NONE

      TYPE, EXTENDS(ARPACKBASE) :: ARPACKBASE_EXT
          REAL(DKIND), ALLOCATABLE :: A(:,:)
          INTEGER :: NV = 0
          INTEGER :: MM = 0
          INTEGER :: NN = 0
      END TYPE ARPACKBASE_EXT 
      CONTAINS
    
      FUNCTION FCN(X, TASK, EXTYPE) RESULT(Y)
         USE UMACH_INT
         CLASS (ARPACKBASE), INTENT(INOUT) :: EXTYPE
         REAL (DKIND), INTENT(INOUT) :: X(EXTYPE % NCOLS)
         INTEGER, INTENT(IN) :: TASK
         REAL (DKIND) Y(EXTYPE % MROWS)
         INTEGER I, J, NOUT
         
         CALL UMACH(2, NOUT)
         SELECT TYPE(EXTYPE)
            TYPE IS(ARPACKBASE_EXT)
            ASSOCIATE(M => EXTYPE % MM,&
                      N => EXTYPE % NN,&
                      A => EXTYPE % A)
            SELECT CASE(TASK)
            CASE(ARPACK_A_x)
!     Computes y <-- A*x		    
               CALL DGEMV('N',M,N,1._DKIND,A,M,X,1,0._DKIND,Y,1)
               EXTYPE % NV = EXTYPE % NV + 1
            CASE(ARPACK_xt_A)

!     Computes y <-- A^T*x = x^T * A
               CALL DGEMV('T',M,N,1._DKIND,A,M,X,1,0._DKIND,Y,1)
               EXTYPE % NV = EXTYPE % NV + 1
 
            CASE(ARPACK_Prepare)
               EXTYPE % NV = 0
               IF(ALLOCATED(EXTYPE % A)) DEALLOCATE(EXTYPE % A)
               ALLOCATE(EXTYPE % A(M,N))
               
               DO J=1,N
                  DO I=1,M
                     EXTYPE % A (I,J) = I + J
                  END DO
               END DO
            CASE DEFAULT
               WRITE(NOUT,*) TASK, ': INVALID TASK REQUESTED'
               STOP 'IMSL_ERROR_WRONG_OPERATION'
            END SELECT
            
            END ASSOCIATE
         END SELECT
      END FUNCTION
      END MODULE

      USE ARPACK_SVD_EX1_INT
      USE UMACH_INT
      USE WRRRN_INT
! Compute the largest and smallest singular values of a 
! patterned matrix.
      INTEGER, PARAMETER :: NSV=2
      INTEGER ::  COUNT, I, J, N, M, nout
      REAL(DKIND)  :: SVALUESMax(NSV)
      REAL(DKIND), ALLOCATABLE :: SVALUEsMin(:)
      REAL(DKIND), ALLOCATABLE :: VECTORS(:,:), B(:,:)
      REAL(DKIND), ALLOCATABLE :: U_VECTORS(:,:), V_VECTORS(:,:)
      REAL(DKIND) NORM
      LOGICAL SMALL, SOLVED
      TYPE(ARPACKBASE_EXT) EX

      ASSOCIATE(M=>EX % MM,&
                N=>EX % NN,&
                NACC=>EX % NACC,&
                TOL =>EX % TOL,&
                MAXMV => EX % MAXMV)
      SOLVED = .true.          
      CALL UMACH(2, NOUT)       
! Define size of matrix problem.
      N=800
      M=600
      DO COUNT =1,2
! Some values will not be accurately determined for rank
! deficient problems.  This next value drops the number
! requested after every sequence of iterations of this size.    
         MAXMV=500
         CALL ARPACK_SVD(M, N, FCN, SVALUESMax, &
                PLACE=ARPACK_Largest_Algebraic, & !Default
                Iteration_TYPE=ARPACK_Normal,   & !Default
                CATEGORY=ARPACK_Regular,        & !Default
                EXTYPE=EX, U_VECTORS=U_VECTORS,   &
                V_VECTORS=V_VECTORS)
         CALL WRRRN( 'Largest Singular values, Normal Method', &
                     SVALUESMax)
         WRITE(NOUT, *) 'Number of matrix-vector products'
         WRITE(NOUT, *) '--------------------------------'
         WRITE(NOUT, '(5X, I4)') EX % NV
         IF(ALLOCATED(B))DEALLOCATE(B)
         ALLOCATE(B(M,N))
! Reconstruct an approximation to A, B = U * S * V ^T.  
! Use only the singular values accurately determined. 
     
         DO I=1,NACC
            U_VECTORS(:,I)=U_VECTORS(:,I)*SVALUESMax(I)
         END DO
         B=matmul(U_VECTORS(:,1:NACC),transpose(V_VECTORS(:,1:NACC)))
! Truncate the approximation to nearest integers.       
! Subtract known integer matrix and check agreement with
! the approximation.
         DO I=1,M
            DO J=1,N
               B(I,J)=REAL(NINT(B(I,J)),DKIND)
               B(I,J)=B(I,J)-EX % A(I,J)
            END DO
         END DO    
         WRITE(NOUT,'(/A,I6)')&
         'Number of singular values, S and columns of U,V =', NACC
         WRITE(NOUT,'(/A,F6.2)')&
               'Integer units of error with U,V and S =', maxval(B)
         if (maxval(B) > 0.0d0) then
            solved = .false.
         else 
            solved = solved .and. .true.
         end if
         SVALUESMax=0._DKIND
! Do same SVD with the Expanded form of the symmetric matrix. 
         CALL ARPACK_SVD(M, N, FCN, SVALUESMax,&
                         PLACE=ARPACK_Largest_Algebraic, & !Default
                         Iteration_TYPE=ARPACK_Expanded, &    
                         CATEGORY=ARPACK_Regular,        & !Default
                         EXTYPE=EX, U_VECTORS=U_VECTORS,   &
                         V_VECTORS=V_VECTORS)
                         
         CALL WRRRN('Largest Singular values, Expanded Method', SVALUESMax)
         WRITE(NOUT, *) 'Number of matrix-vector products'
         WRITE(NOUT, *) '--------------------------------'
         WRITE(NOUT, '(5X,I4)') EX % NV
! Reconstruct an approximation to A, B = U * S * V ^T.  
! Use only the singular values accurately determined.  
         DO I=1,NACC
            U_VECTORS(:,I)=U_VECTORS(:,I)*SVALUESMax(I)
         END DO
         B=matmul(U_VECTORS(:,1:NACC),transpose(V_VECTORS(:,1:NACC)))
! Truncate the approximation to nearest integers.       
! Subtract known integer matrix and check agreement with
! the approximation. 
         DO I=1,M
            DO J=1,N
               B(I,J)=REAL(NINT(B(I,J)),DKIND)
               B(I,J)=B(I,J)-EX % A(I,J)
            END DO
         END DO      
         WRITE(NOUT,'(A,I6)')&
              'Number of singular values, S and columns of U,V =', NACC
         WRITE(NOUT,'(A,F6.2)')&
               'Integer units of error with U,V and S =', maxval(B)
         if (maxval(B) > 0.0d0) then
            solved = .false.
         else 
            solved = solved .and. .true.
         end if

         M=800
         N=600
         DEALLOCATE(U_VECTORS, V_VECTORS)
      END DO         
      
      END ASSOCIATE
      END
