Developer Reference

  • 0.10
  • 10/21/2020
  • Public Content
Contents

?gesvj

Computes the singular value decomposition of a real matrix using Jacobi plane rotations.

Syntax

lapack_int
LAPACKE_sgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
float
*
a
,
lapack_int
lda
,
float
*
sva
,
lapack_int
mv
,
float
*
v
,
lapack_int
ldv
,
float
*
stat
);
lapack_int
LAPACKE_dgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
double
*
a
,
lapack_int
lda
,
double
*
sva
,
lapack_int
mv
,
double
*
v
,
lapack_int
ldv
,
double
*
stat
);
lapack_int
LAPACKE_cgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
lapack_complex_float
*
a
,
lapack_int
lda
,
float
*
sva
,
lapack_int
mv
,
lapack_complex_float
*
v
,
lapack_int
ldv
,
float
*
stat
);
lapack_int
LAPACKE_zgesvj
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
lapack_int
m
,
lapack_int
n
,
lapack_complex_double
*
a
,
lapack_int
lda
,
double
*
sva
,
lapack_int
mv
,
lapack_complex_double
*
v
,
lapack_int
ldv
,
double
*
stat
);
Include Files
  • mkl.h
Description
The routine computes the singular value decomposition (SVD) of a real or complex
m
-by-
n
matrix
A
, where
m
n
.
The SVD of
A
is written as
A
=
U
*
Σ
*
V
T
for real flavors, or
A
=
U
*
Σ
*
V
H
for complex flavors,
where
Σ
is an
m
-by-
n
diagonal matrix,
U
is an
m
-by-
n
orthonormal matrix, and
V
is an
n
-by-
n
orthogonal/unitary matrix. The diagonal elements of
Σ
are the singular values of
A
; the columns of
U
and
V
are the left and right singular vectors of
A
, respectively. The matrices
U
and
V
are computed and stored in the arrays
u
and
v
, respectively. The diagonal of
Σ
is computed and stored in the array
sva
.
The
?gesvj
routine can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines.
The
n
-by-
n
orthogonal matrix
V
is obtained as a product of Jacobi plane rotations. The rotations are implemented as fast scaled rotations of Anda and Park [AndaPark94]. In the case of underflow of the Jacobi angle, a modified Jacobi transformation of Drmac ([Drmac08-4]) is used. Pivot strategy uses column interchanges of de Rijk ([deRijk98]). The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [Demmel92]. The condition number that determines the accuracy in the full rank case is essentially
where
κ
(.) is the spectral condition number. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [Drmac08-1], [Drmac08-2].
The computational range for the nonzero singular values is the machine number interval (
UNDERFLOW
,
OVERFLOW
). In extreme cases, even denormalized singular values can be computed with the corresponding gradual loss of accurate digit.
Input Parameters
matrix_layout
Specifies whether matrix storage layout is row major (
LAPACK_ROW_MAJOR
) or column major (
LAPACK_COL_MAJOR
).
joba
Must be
'L'
,
'U'
or
'G'
.
Specifies the structure of
A
:
If
joba
=
'L'
, the input matrix
A
is lower triangular.
If
joba
=
'U'
, the input matrix
A
is upper triangular.
If
joba
=
'G'
, the input matrix
A
is a general
m
-by-
n
,
m
n
.
jobu
Must be
'U'
,
'C'
or
'N'
.
Specifies whether to compute the left singular vectors (columns of
U
):
If
jobu
=
'U'
, the left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of
A
. See more details in the description of
a
. The default numerical orthogonality threshold is set to approximately
TOL=CTOL*EPS, CTOL=sqrt(
m
), EPS =
?lamch
(
'E'
)
If
jobu
=
'C'
, analogous to
jobu
=
'U'
, except that you can control the level of numerical orthogonality of the computed left singular vectors.
TOL
can be set to
TOL=CTOL*EPS
, where
CTOL
is given on input in the array
stat
. No
CTOL
smaller than ONE is allowed.
CTOL
greater than
1 / EPS
is meaningless. The option
'C'
can be used if
m
*EPS
is satisfactory orthogonality of the computed left singular vectors, so
CTOL=
m
could save a few sweeps of Jacobi rotations. See the descriptions of
a
and
stat[0]
.
If
jobu
=
'N'
,
u
is not computed. However, see the description of
a
.
jobv
Must be
'V'
,
'A'
or
'N'
.
Specifies whether to compute the right singular vectors, that is, the matrix
V
:
If
jobv
=
'V'
, the matrix
V
is computed and returned in the array
v
.
If
jobv
=
'A'
, the Jacobi rotations are applied to the
mv
-by
n
array
v
. In other words, the right singular vector matrix
V
is not computed explicitly, instead it is applied to an
mv
-by
n
matrix initially stored in the first
mv
rows of
V
.
If
jobv
=
'N'
, the matrix
V
is not computed and the array
v
is not referenced.
m
The number of rows of the input matrix
A
.
1/slamch('E')>
m
0
for
sgesvj
.
1/dlamch('E')>
m
0
for
dgesvj
.
n
The number of columns in the input matrix
A
;
m
n
 0.
a
,
v
Array
a
(size at least
lda
*
n
for column major layout and
lda
*
m
for row major layout)
is an array containing the
m
-by-
n
matrix
A
.
Array
v
(size at least max(1,
ldv
*
n
)) contains, if
jobv
= 'A'
the
mv
-by-
n
matrix to be post-multiplied by Jacobi rotations
.
lda
The leading dimension of the array
a
. Must be at least
max(1,
m
)
for column major layout and at least max(1,
n
) for row major layout
.
mv
If
jobv
=
'A'
, the product of Jacobi rotations in
?gesvj
is applied to the first
mv
rows of
v
. See the description of
jobv
.
0
mv
ldv
.
ldv
The leading dimension of the array
v
;
ldv
1
.
jobv
=
'V'
,
ldv
max(1,
n
)
.
jobv
=
'A'
,
ldv
max(1,
mv
)
for column major layout and
ldv
max(1,
n
)
for row major layout
.
stat
Array size 6. If
jobu
= 'C'
,
stat
[0] =
CTOL
, where
CTOL
defines the threshold for convergence. The process stops if all columns of
A
are mutually orthogonal up to
CTOL
*
EPS
, where
EPS
=
?lamch
('E')
. It is required that
CTOL
1 - that is, it is not allowed to force the routine to obtain orthogonality below
ε
.
Output Parameters
a
On exit:
If
jobu
=
'U'
or
jobu
=
'C'
:
  • if
    info
    = 0
    , the leading columns of
    A
    contain left singular vectors corresponding to the computed singular values of
    a
    that are above the underflow threshold
    ?lamch
    (
    'S'
    ), that is, non-zero singular values. The number of the computed non-zero singular values is returned in
    stat[1]
    . Also see the descriptions of
    sva
    and
    stat
    . The computed columns of
    u
    are mutually numerically orthogonal up to approximately
    TOL=sqrt(
    m
    )*EPS
    (default); or
    TOL=CTOL*EPS
    jobu
    =
    'C'
    , see the description of
    jobu
    .
  • if
    info
    > 0
    , the procedure
    ?gesvj
    did not converge in the given number of iterations (sweeps). In that case, the computed columns of
    u
    may not be orthogonal up to
    TOL
    . The output
    u
    (stored in
    a
    ),
    sigma
    (given by the computed singular values in
    sva
    (1:n)
    ) and
    v
    is still a decomposition of the input matrix
    A
    in the sense that the residual
    ||
    A
    -scale*
    U
    *
    sigma
    *
    V
    T
    ||
    2
    / ||
    A
    ||
    2
    for real flavors or
    ||
    A
    -scale*
    U
    *
    sigma
    *
    V
    H
    ||
    2
    / ||
    A
    ||
    2
    for complex flavors (where
    scale =
    stat
    [0]
    ) is small.
If
jobu
=
'N'
:
  • if
    info
    = 0
    , note that the left singular vectors are 'for free' in the one-sided Jacobi SVD algorithm. However, if only the singular values are needed, the level of numerical orthogonality of
    u
    is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately
    m
    *EPS
    . Thus, on exit,
    a
    contains the columns of
    u
    scaled with the corresponding singular values.
  • if
    info
    > 0
    , the procedure
    ?gesvj
    did not converge in the given number of iterations (sweeps).
sva
Array size
n
.
If
info
= 0
, depending on the value
scale
=
stat[0]
, where
scale
is the scaling factor:
  • if
    scale
    = 1
    ,
    sva
    [0:
    n
    - 1]
    contains the computed singular values of
    a
    .
  • if
    scale
    1
    , the singular values of
    a
    are
    scale
    *
    sva
    (1:n)
    , and this factored representation is due to the fact that some of the singular values of
    a
    might underflow or overflow.
If
info
> 0
, the procedure
?gesvj
did not converge in the given number of iterations (sweeps) and
scale
*
sva
(1:n)
may not be accurate.
v
On exit:
If
jobv
=
'V'
, contains the
n
-by-
n
matrix of the right singular vectors.
If
jobv
=
'A'
, then
v
contains the product of the computed right singular vector matrix and the initial matrix in the array
v
.
If
jobv
=
'N'
,
v
is not referenced.