NAG Library Routine Document
F02XUF
1 Purpose
F02XUF returns all, or part, of the singular value decomposition of a complex upper triangular matrix.
2 Specification
SUBROUTINE F02XUF ( |
N, A, LDA, NCOLB, B, LDB, WANTQ, Q, LDQ, SV, WANTP, RWORK, CWORK, IFAIL) |
INTEGER |
N, LDA, NCOLB, LDB, LDQ, IFAIL |
REAL (KIND=nag_wp) |
SV(N), RWORK(*) |
COMPLEX (KIND=nag_wp) |
A(LDA,*), B(LDB,*), Q(LDQ,*), CWORK(max(1,N-1)) |
LOGICAL |
WANTQ, WANTP |
|
3 Description
The
by
upper triangular matrix
is factorized as
where
and
are
by
unitary matrices and
is an
by
diagonal matrix with real 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.
F02XUF 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.
Good background descriptions to the singular value decomposition are given in
Dongarra et al. (1979),
Hammarling (1985) and
Wilkinson (1978).
Note that if
is any unitary diagonal matrix so that
then
is also a singular value decomposition of
.
4 References
Dongarra J J, Moler C B, Bunch J R and Stewart G W (1979) LINPACK Users' Guide SIAM, Philadelphia
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,) – COMPLEX (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
unitary 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 F02XUF 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,) – COMPLEX (KIND=nag_wp) arrayInput/Output
-
Note: the second dimension of the array
B
must be at least
.
On entry: if
, the leading
by
part of the array
B must contain the matrix to be transformed.
On exit: is overwritten by the by matrix .
- 6: LDB – INTEGERInput
On entry: the first dimension of the array
B as declared in the (sub)program from which F02XUF is called.
Constraints:
- if , ;
- otherwise .
- 7: WANTQ – LOGICALInput
On entry: must be .TRUE. if the matrix
is required.
If
then the array
Q is not referenced.
- 8: Q(LDQ,) – COMPLEX (KIND=nag_wp) arrayOutput
-
Note: the second dimension of the array
Q
must be at least
if
, and at least
otherwise.
On exit: if
, the leading
by
part of the array
Q will contain the unitary 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 F02XUF is called.
Constraints:
- if , ;
- otherwise .
- 10: SV(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the diagonal elements of the matrix .
- 11: WANTP – LOGICALInput
On entry: must be .TRUE. if the matrix
is required, in which case
is returned in the array
A, otherwise
WANTP must be .FALSE..
- 12: RWORK() – REAL (KIND=nag_wp) arrayOutput
-
Note: the dimension of the array
RWORK
must be at least
if
and
and
,
if
and
and
or
and
or
and
, and at least
otherwise.
On exit:
RWORK(
N) contains the total number of iterations taken by the
algorithm.
The rest of the array is used as workspace.
- 13: CWORK() – COMPLEX (KIND=nag_wp) arrayWorkspace
- 14: 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
.
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
, where
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 (f02xufe.f90)
9.2 Program Data
Program Data (f02xufe.d)
9.3 Program Results
Program Results (f02xufe.r)