Developer Reference for Intel® oneAPI Math Kernel Library - C

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.