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.
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).
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.
Generic: ARPACK_SVD (M, N, F,SVALUES [,…])
Specific: The specific interface name is D_ARPACK_SVD.
A Fortran 90 compiler version is not available for this routine.
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.
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
, where
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 ,
is computed from the
relation
. The result is then
processed with the modified Gram-Schmidt algorithm to assure that
is orthogonal. When iterating with
we first compute the left singular vectors
and then obtain
by the Gram-Schmidt
algorithm. If we use Iteration_Type=ARPACK_Expanded,
and
are computed simultaneously, and both are orthogonal.
We define the matrix
with 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 matrix
has its storage requirements changed from
to 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
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
PHONE: 713.784.3131 FAX:713.781.9260 |