ARPACK_SVD

Computes some singular values and left and right singular vectors of a real rectangular matrix AMxN = USVT. There is no restriction on the relative sizes, M and N. 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 operations y = Ax and y = ATx ≡ xTA. 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	y = Ax
ARPACK_xt_A	y = xTA

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 M ≥ N, ARPACK_Normal specifies computing singular values based on eigenvalues and eigenvectors of the normal symmetric matrix ATA; for values M < N this will be the alternate symmetric matrix AAT.

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

Default: ITERATION_TYPE = ARPACK_Normal.

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 and ITERATION_TYPE 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,x,y, switch between the values n, m. The values n, m 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 A. This approximation is B = USVT. The matrices U, S and V are available in these respective routine arguments. The terms UMxNSV and VMxNSV have orthogonal columns, UTU = INSV = VTV. The diagonal matrix SNSVxNSV has its entries in SVALUES(:), ordered from largest to smallest. Use the value NSV = min(size(SVALUES), NACC), where NACC is 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 ATA, U is computed from the relation AV = US. The result is then processed with the modified Gram-Schmidt algorithm to assure that U is orthogonal. When iterating with AAT we first compute the left singular vectors U and then obtain V by the Gram-Schmidt algorithm. If we use Iteration_Type=ARPACK_Expanded, U and V are computed simultaneously, and both are orthogonal.

Example

We define the M ×N matrix A = {ai,j} with entries ai,j = i + j. This matrix has two non-zero singular values. With the pair of values (M,N) = (512,265) and (M,N) = (265,512) 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 matrix A has its storage requirements changed from M,N to the value 2(M + N + 1). The resulting product B = USVT, when rounded to the nearest integer, satisfies B = A.

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

(Example 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

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

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

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, &

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

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

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

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

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

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

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

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

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

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

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