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