Chapter 2: Eigensystem Analysis

ARPACK_SVD

Computes some singular values and left and right singular vectors of a real rectangular matrix.  There is no restriction on the relative sizes,  and.  This routine calls ARPACK_SYMMETRIC.

Required Arguments

M — The number of matrix rows.   (Input)

N — The number of matrix columns.   (Input)

F — User-supplied FUNCTION to return matrix-vector operations. This user function is written corresponding to the abstract interface for the function SVDMV(…).The operations provided as code in the function F will be made based on the two matrix operationsand. The usage is F (X, TASK, EXTYPE), where

Function Return Value

F — An array of length N containing matrix-vector operations or linear equations solutions. Operations provided as code in the function F will be made depending upon the value of argument TASK.

Required Arguments

X — An array of length N containing the vector to which the operator will be applied.   (Input)

TASK — An enumerated type which specifies the operation to be performed.   (Input)
TASK
is an enumerated integer value, use-associated from the module ARPACK_INT. It will be one of the following:

 

Value

Description

ARPACK_Prepare

Take initial steps to prepare for the operations that follow. These steps can include defining the data for the matrices, factorizations for upcoming linear system solves, or recording the vectors used in the operations.

ARPACK_A_x

ARPACK_xt_A

 

EXTYPE — A derived type of the extensible class ARPACKBASE, which may be used to pass additional information to/from the user-supplied function.   (Input/Output)
The user must include a USE ARPACK_INT statement in the calling program to define this derived type. If EXTYPE is not included as an argument to ARPACK_SVD it should be ignored in the user-function, F.

            The function F must be written according to the abstract interface for SVDMV. If F is not use-associated nor contained in the calling program, declare it with PROCEDURE(SVDMV) F.

SVALUES — A rank-1 array of singular values.   (Output)
The value NEV = size(SVALUES) defines the number of singular values to be computed.  The calling program declares or allocates the array SVALUES(1:NEV).

Optional Arguments

PLACE — Indicates the placement in the spectrum for required singular values.   (Input)
PLACE can be one of the following enumerated integers as defined in ARPACK_INT:

 

Value

ARPACK_Largest_Algebraic

ARPACK_Smallest_Magnitude

 

            Default: PLACE = ARPACK_Largest_Algebraic.

ITERATION_TYPE — Indicates the method for obtaining the required singular values.   (Input)
ITERATION_TYPE can be one of the following enumerated integers as defined in ARPACK_INT:

 

Value

ARPACK_Normal

ARPACK_Expanded

 

            For values  , ARPACK_Normal specifies computing singular values based on eigenvalues and eigenvectors of the normal symmetric matrix; for values M  N this will be the alternate symmetric matrix.

            For all values of , ARPACK_Expanded specifies computing singular values based on the symmetric matrix eigenvalue problem for the matrices

  or 

            Default: ITERATION_TYPE = ARPACK_Normal.

CATEGORY — An integer from a packaged enumeration with values that are passed to ARPACK_SYMMETRIC.   (Input)
CATEGORY can be one of the following enumerated integers as defined in ARPACK_INT:

 

Value

ARPACK_Regular

ARPACK_Regular_Inverse

 

            Default: CATEGORY = Regular.

EXTYPE — A derived type of the extensible class ARPACKBASE, which may be used to pass additional information to/from the user-supplied function.   (Input/Output)
The user must include a USE ARPACK_INT statement in the calling program to define this derived type. If EXTYPE is not included as an argument to ARPACK_SVD it must still be supplied as an argument to user-function, F, but is not used.

U_VECTORS — An allocatable array of orthogonal left singular vectors.   (Output)
It is not necessary to allocate U_VECTORS(:,:).  If this argument is present, the allocation occurs within the routine ARPACK_SVD  The output sizes are UVECTORS(1:M,1:NCV).  The second dimension value is NCV=min(M, max(FACTOR_MAXNCV*NEV,NEV+1)), where the value FACTOR_MAXNCV is a component of the base class, ARPACKBASE.  The first NEV columns of U_VECTORS(:,:) are the left singular vectors.

V_VECTORS — An allocatable array of orthogonal right singular vectors.   (Output)
It is not necessary to allocate V_VECTORS(:,:).  If this argument is present, the allocation occurs within the routine ARPACK_SVD.  The output sizes are V_VECTORS(1:N,1:NCV).  The second dimension value is NCV=min(M, max(FACTOR_MAXNCV*NEV,NEV+1)), where the value FACTOR_MAXNCV is a component of the base class, ARPACKBASE.  The first NEV columns of V_VECTORS(:,:) are the right singular vectors.

FORTRAN 2003 Interface

Generic:          ARPACK_SVD (M, N, F,SVALUES [,…])

Specific:         The specific interface name is D_ARPACK_SVD.

FORTRAN 90 Interface

A Fortran 90 compiler version is not available for this routine.

Description

Routine ARPACK_SVD calls ARPACK_SYMMETRIC to compute partial singular value decompositions for rectangular real matrices.  There is no restriction on the relative sizes of the number of rows and columns.  A function internal to ARPACK_SVD is used in the call to ARPACK_SYMMETRIC.  The internal function calls the user function, F, which provides matrix-vector products of the matrix and an internally generated vector.  By including the class object EXTYPE as an argument to this function, user data and procedure pointers are available for the evaluations.  A user code must extend the base class EXTYPE to include the extra data and procedure pointers.

Comments

The user function supplies requests for the matrix operations.  Those cases that follow from the settings of PLACE, Iteration_Type and Category need to be provided in the user code. The enumerated TASK constants, ARPACK_A_x and ARPACK_xt_A are available by use-association from the module ARPACK_INT.  The sizes of the inputs and outputs, , switch between the values.  The values are alternated in the base class components EXTYPE%NCOLS and EXTYPE%MROWS.

The value of Iteration_Type may impact the number of iterations required.  Generally one expects Iteration_Type=ARPACK_Normal (the default) to result in the fewest iterations, and Iteration_Type=ARPACK_Expanded to result in singular values with the greatest accuracy.

The output arrays U_VECTORS(:,:), SVALUES(:), and V_VECTORS(:,:) allow for reconstruction of an approximation to the matrix.  This approximation is.  The matrices  and  are available in these respective routine arguments. The terms and  have orthogonal columns, .  The diagonal matrix  has its entries in SVALUES(:), ordered from largest to smallest.  Use the value, whereis the number of accurate singular values computed by ARPACK_SYMMETRIC.  This is the component EXTYPE%NACC of the base class EXTYPE.

After computing the singular values and right singular vectors by iteration with the normal matrix , is computed from the relation .  The result is then processed with the modified Gram-Schmidt algorithm to assure thatis orthogonal.  When iterating withwe first compute the left singular vectorsand then obtainby the Gram-Schmidt algorithm.  If we use Iteration_Type=ARPACK_Expanded, and are computed simultaneously, and both are orthogonal.

Example 1

We define the matrixwith entries.  This matrix has two non-zero singular values.  With the pair of values  and  we obtain the singular decomposition for these rectangular matrices.   With each pair we compute the decomposition using the input Iteration_Type=ARPACK_Normal and Iteration_Type=ARPACK_Expanded.  The latter value requires the larger number of iterations.  The matrixhas its storage requirements changed fromto the value.  The resulting product, when rounded to the nearest integer, satisfies.

The base class ARPACKBASE is extended to include an allocatable array, EXTYPE%A(:,:).  This is allocated and defined and stores the matrix.  The matrix operations and are computed with DGEMV.

Link to example source (arpack_svd_ex1.f90)

 

      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

Output

 Largest Singular values, Normal Method

              1   523955.7

              2    36644.2

 Number of matrix-vector products

 --------------------------------

       12

 

Number of singular values, S and columns of U,V =     2

 

Integer units of error with U,V and S =  0.00

 Largest Singular values, Expanded Method

               1   523955.7

               2    36644.2

 Number of matrix-vector products

 --------------------------------

       22

Number of singular values, S and columns of U,V =     2

Integer units of error with U,V and S =  0.00

 Largest Singular values, Normal Method

              1   523955.7

              2    36644.2

 Number of matrix-vector products

 --------------------------------

       12

 

Number of singular values, S and columns of U,V =     2

 

Integer units of error with U,V and S =  0.00

 Largest Singular values, Expanded Method

               1   523955.7

               2    36644.2

 Number of matrix-vector products

 --------------------------------

       18

Number of singular values, S and columns of U,V =     2

Integer units of error with U,V and S =  0.00



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260