NAG Library Routine Document
F02WUF
1 Purpose
F02WUF returns all, or part, of the singular value decomposition of a real upper triangular matrix.
2 Specification
SUBROUTINE F02WUF ( |
N, A, LDA, NCOLB, B, LDB, WANTQ, Q, LDQ, SV, WANTP, WORK, IFAIL) |
INTEGER |
N, LDA, NCOLB, LDB, LDQ, IFAIL |
REAL (KIND=nag_wp) |
A(LDA,*), B(LDB,*), Q(LDQ,*), SV(N), WORK(*) |
LOGICAL |
WANTQ, WANTP |
|
3 Description
The
by
upper triangular matrix
is factorized as
where
and
are
by
orthogonal matrices and
is an
by
diagonal matrix with non-negative diagonal elements,
, ordered such that
The columns of
are the left-hand singular vectors of
, the diagonal elements of
are the singular values of
and the columns of
are the right-hand singular vectors of
.
Either or both of
and
may be requested and the matrix
given by
where
is an
by
given matrix, may also be requested.
The routine obtains the singular value decomposition by first reducing to bidiagonal form by means of Givens plane rotations and then using the algorithm to obtain the singular value decomposition of the bidiagonal form.
Note that if
is any orthogonal diagonal matrix so that
(that is the diagonal elements of
are
or
) then
is also a singular value decomposition of
.
4 References
Chan T F (1982) An improved algorithm for computing the singular value decomposition ACM Trans. Math. Software 8 72–83
Dongarra J J, Moler C B, Bunch J R and Stewart G W (1979) LINPACK Users' Guide SIAM, Philadelphia
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Hammarling S (1985) The singular value decomposition in multivariate statistics SIGNUM Newsl. 20(3) 2–25
Wilkinson J H (1978) Singular Value Decomposition – Basic Aspects Numerical Software – Needs and Availability (ed D A H Jacobs) Academic Press
5 Parameters
- 1: N – INTEGERInput
On entry:
, the order of the matrix
.
If , an immediate return is effected.
Constraint:
.
- 2: A(LDA,) – REAL (KIND=nag_wp) arrayInput/Output
-
Note: the second dimension of the array
A
must be at least
.
On entry: the leading
by
upper triangular part of the array
A must contain the upper triangular matrix
.
On exit: if
, the
by
part of
A will contain the
by
orthogonal matrix
, otherwise the
by
upper triangular part of
A is used as internal workspace, but the strictly lower triangular part of
A is not referenced.
- 3: LDA – INTEGERInput
On entry: the first dimension of the array
A as declared in the (sub)program from which F02WUF is called.
Constraint:
.
- 4: NCOLB – INTEGERInput
On entry:
, the number of columns of the matrix
.
If
, the array
B is not referenced.
Constraint:
.
- 5: B(LDB,) – REAL (KIND=nag_wp) arrayInput/Output
-
Note: the second dimension of the array
B
must be at least
.
On entry: with
, the leading
by
part of the array
B must contain the matrix to be transformed.
On exit: the leading
by
part of the array
B is overwritten by the matrix
.
- 6: LDB – INTEGERInput
On entry: the first dimension of the array
B as declared in the (sub)program from which F02WUF is called.
Constraints:
- if , ;
- otherwise .
- 7: WANTQ – LOGICALInput
On entry: must be .TRUE. if the matrix
is required.
If
, the array
Q is not referenced.
- 8: Q(LDQ,) – REAL (KIND=nag_wp) arrayOutput
-
Note: the second dimension of the array
Q
must be at least
if
, and at least
otherwise.
On exit: with
, the leading
by
part of the array
Q will contain the orthogonal matrix
. Otherwise the array
Q is not referenced.
- 9: LDQ – INTEGERInput
On entry: the first dimension of the array
Q as declared in the (sub)program from which F02WUF is called.
Constraints:
- if , ;
- otherwise .
- 10: SV(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the array
SV will contain the
diagonal elements of the matrix
.
- 11: WANTP – LOGICALInput
On entry: must be .TRUE. if the matrix
is required, in which case
is overwritten on the array
A, otherwise
WANTP must be .FALSE..
- 12: WORK() – REAL (KIND=nag_wp) arrayOutput
-
Note: the dimension of the array
WORK
must be at least
if
and
and
,
if (
and
and
) or (
and (
or
)), and at least
otherwise.
On exit:
contains the total number of iterations taken by the
algorithm.
The rest of the array is used as internal workspace.
- 13: IFAIL – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
.
When the value is used it is essential to test the value of IFAIL on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
On entry, | , |
or | , |
or | , |
or | and , |
or | and . |
The algorithm has failed to converge in iterations. In this case may not have been found correctly and the remaining singular values may not be the smallest. The matrix will nevertheless have been factorized as , where is a bidiagonal matrix with as the diagonal elements and as the superdiagonal elements.
This failure is not likely to occur.
7 Accuracy
The computed factors
,
and
satisfy the relation
where
is the
machine precision,
is a modest function of
and
denotes the spectral (two) norm. Note that
.
A similar result holds for the computed matrix .
The computed matrix
satisfies the relation
where
is exactly orthogonal and
where
is a modest function of
. A similar result holds for
.
For given values of
NCOLB,
WANTQ and
WANTP, the number of floating point operations required is approximately proportional to
.
Following the use of this routine the rank of
may be estimated by a call to the INTEGER FUNCTION
F06KLF. The statement
IRANK = F06KLF(N,SV,1,TOL)
returns the value
in
IRANK, where
is the smallest integer for which
, and
is the tolerance supplied in
TOL, so that
IRANK is an estimate of the rank of
and thus also of
. If
TOL is supplied as negative then the
machine precision is used in place of
TOL .
9 Example
This example finds the singular value decomposition of the
by
upper triangular matrix
together with the vector
for the vector
9.1 Program Text
Program Text (f02wufe.f90)
9.2 Program Data
Program Data (f02wufe.d)
9.3 Program Results
Program Results (f02wufe.r)