Developer Reference for Intel® oneAPI Math Kernel Library - C

Developer Reference

Contents

?gejsv

Computes the singular value decomposition using a preconditioned Jacobi SVD method.

Syntax

lapack_int
LAPACKE_sgejsv
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
char
jobr
,
char
jobt
,
char
jobp
,
lapack_int
m
,
lapack_int
n
,
float
*
a
,
lapack_int
lda
,
float
*
sva
,
float
*
u
,
lapack_int
ldu
,
float
*
v
,
lapack_int
ldv
,
float
*
stat
,
lapack_int
*
istat
);
lapack_int
LAPACKE_dgejsv
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
char
jobr
,
char
jobt
,
char
jobp
,
lapack_int
m
,
lapack_int
n
,
double
*
a
,
lapack_int
lda
,
double
*
sva
,
double
*
u
,
lapack_int
ldu
,
double
*
v
,
lapack_int
ldv
,
double
*
stat
,
lapack_int
*
istat
);
lapack_int
LAPACKE_cgejsv
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
char
jobr
,
char
jobt
,
char
jobp
,
lapack_int
m
,
lapack_int
n
,
lapack_complex_float
*
a
,
lapack_int
lda
,
float
*
sva
,
lapack_complex_float
*
u
,
lapack_int
ldu
,
lapack_complex_float
*
v
,
lapack_int
ldv
,
float
*
stat
,
lapack_int
*
istat
);
lapack_int
LAPACKE_zgejsv
(
int
matrix_layout
,
char
joba
,
char
jobu
,
char
jobv
,
char
jobr
,
char
jobt
,
char
jobp
,
lapack_int
m
,
lapack_int
n
,
lapack_complex_double
*
a
,
lapack_int
lda
,
double
*
sva
,
lapack_complex_double
*
u
,
lapack_int
ldu
,
lapack_complex_double
*
v
,
lapack_int
ldv
,
double
*
stat
,
lapack_int
*
istat
);
Include Files
  • mkl.h
Description
The routine computes the singular value decomposition (SVD) of a real/complex
m
-by-
n
 matrix
A
, where
m
≥
n
.
The SVD is written as
A
 =
U
*
Σ
*
V
T
, for real routines
A
 =
U
*
Σ
*
V
H
, for complex routines
where
Σ
 is an
m
-by-
n
 matrix which is zero except for its
n
 diagonal elements,
U
 is an
m
-by-
n
 (or
m
-by-
m
) orthonormal matrix, and
V
 is an
n
-by-
n
 orthogonal 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
?gejsv
 routine can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines.
The routine implements a preconditioned Jacobi SVD algorithm. It uses
?geqp3
,
?geqrf
, and
?gelqf
 as preprocessors and preconditioners. Optionally, an additional row pivoting can be used as a preprocessor, which in some cases results in much higher accuracy. An example is matrix
A
 with the structure
A = D1 * C * D2
, where
D1
,
D2
 are arbitrarily ill-conditioned diagonal matrices and
C
 is a well-conditioned matrix. In that case, complete pivoting in the first QR factorizations provides accuracy dependent on the condition number of
C
, and independent of
D1
,
D2
. Such higher accuracy is not completely understood theoretically, but it works well in practice.
If
A
 can be written as
A = B*D
, with well-conditioned
B
 and some diagonal
D
, then the high accuracy is guaranteed, both theoretically and in software, independent of
D
. For more details see [Drmac08-1], [Drmac08-2].
The computational range for the singular values can be the full range (
UNDERFLOW
,
OVERFLOW
 ), provided that the machine arithmetic and the BLAS and LAPACK routines called by
?gejsv
 are implemented to work in that range. If that is not the case, the restriction for safe computation with the singular values in the range of normalized IEEE numbers is that the spectral condition number
kappa(A)=sigma_max(A)/sigma_min(A)
 does not overflow. This code (
?gejsv
) is best used in this restricted range, meaning that singular values of magnitude below
||A||_2 / slamch('O')
 (for single precision) or
||A||_2 / dlamch('O')
 (for double precision) are returned as zeros. See
jobr
 for details on this.
This implementation is slower than the one described in [Drmac08-1], [Drmac08-2] due to replacement of some non-LAPACK components, and because the choice of some tuning parameters in the iterative part (
?gesvj
) is left to the implementer on a particular machine.
The rank revealing QR factorization (in this code:
?geqp3
) should be implemented as in [Drmac08-3].
If
m
 is much larger than
n
, it is obvious that the inital QRF with column pivoting can be preprocessed by the QRF without pivoting. That well known trick is not used in
?gejsv
 because in some cases heavy row weighting can be treated with complete pivoting. The overhead in cases
m
 much larger than
n
 is then only due to pivoting, but the benefits in accuracy have prevailed. You can incorporate this extra QRF step easily and also improve data movement (matrix transpose, matrix copy, matrix transposed copy) - this implementation of
?gejsv
 uses only the simplest, naive data movement.
Product and Performance Information
Performance varies by use, configuration and other factors. Learn more at www.Intel.com/PerformanceIndex.
Notice revision #20201201
Input Parameters
matrix_layout
Specifies whether matrix storage layout is row major (
LAPACK_ROW_MAJOR
) or column major (
LAPACK_COL_MAJOR
).
joba
Must be
'C'
,
'E'
,
'F'
,
'G'
,
'A'
, or
'R'
.
Specifies the level of accuracy:
If
joba
 =
'C'
, high relative accuracy is achieved if
A
 =
B
*
D
 with well-conditioned
B
 and arbitrary diagonal matrix
D
. The accuracy cannot be spoiled by column scaling. The accuracy of the computed output depends on the condition of
B
, and the procedure aims at the best theoretical accuracy. The relative error
max_{i=1:N}|d sigma_i| / sigma_i
 is bounded by
f(M,N)*epsilon* cond(B)
, independent of
D
. The input matrix is preprocessed with the QRF with column pivoting. This initial preprocessing and preconditioning by a rank revealing QR factorization is common for all values of
joba
. Additional actions are specified as follows:
If
joba
 =
'E'
, computation as with
'C'
 with an additional estimate of the condition number of
B
. It provides a realistic error bound.
If
joba
 =
'F'
, accuracy higher than in the
'C'
 option is achieved, if
A
 =
D1
*
C
*
D2
 with ill-conditioned diagonal scalings
D1
,
D2
, and a well-conditioned matrix
C
. This option is advisable, if the structure of the input matrix is not known and relative accuracy is desirable. The input matrix
A
 is preprocessed with QR factorization with full (row and column) pivoting.
If
joba
 =
'G'
, computation as with
'F'
 with an additional estimate of the condition number of
B
, where
A
 =
B
*
D
. If
A
 has heavily weighted rows, using this condition number gives too pessimistic error bound.
If
joba
 =
'A'
, small singular values are the noise and the matrix is treated as numerically rank defficient. The error in the computed singular values is bounded by
f(m,n)*epsilon*||A||
. The computed SVD
A = U*S*V**t
 (for real flavors) or
A = U*S*V**H
 (for complex flavors) restores
A
 up to
f(m,n)*epsilon*||A||
. This enables the procedure to set all singular values below
n*epsilon*||A||
 to zero.
If
joba
 =
'R'
, the procedure is similar to the
'A'
 option. Rank revealing property of the initial QR factorization is used to reveal (using triangular factor) a gap
sigma_{r+1}
<
 epsilon * sigma_r
, in which case the numerical rank is declared to be
r
. The SVD is computed with absolute error bounds, but more accurately than with
'A'
.
jobu
Must be
'U'
,
'F'
,
'W'
, or
'N'
.
Specifies whether to compute the columns of the matrix
U
:
If
jobu
 =
'U'
,
n
 columns of
U
 are returned in the array
u
If
jobu
 =
'F'
, a full set of
m
 left singular vectors is returned in the array
u
.
If
jobu
 =
'W'
,
u
 may be used as workspace of length
m
*
n
. See the description of
u
.
If
jobu
 =
'N'
,
u
 is not computed.
jobv
Must be
'V'
,
'J'
,
'W'
, or
'N'
.
Specifies whether to compute the matrix
V
:
If
jobv
 =
'V'
,
n
 columns of
V
 are returned in the array
v
; Jacobi rotations are not explicitly accumulated.
If
jobv
 =
'J'
,
n
 columns of
V
 are returned in the array
v
 but they are computed as the product of Jacobi rotations. This option is allowed only if
jobu
≠
'N'
If
jobv
 =
'W'
,
v
 may be used as workspace of length
n
*
n
. See the description of
v
.
If
jobv
 =
'N'
,
v
 is not computed.
jobr
Must be
'N'
 or
'R'
.
Specifies the range for the singular values. If small positive singular values are outside the specified range, they may be set to zero. If
A
 is scaled so that the largest singular value of the scaled matrix is around
sqrt(big), big =
?lamch
(
'O'
)
, the function can remove columns of
A
 whose norm in the scaled matrix is less than
sqrt(
?lamch
(
'S'
))
 (for
jobr
 =
'R'
), or less than
small =
?lamch
(
'S'
)/
?lamch
(
'E'
)
.
If
jobr
 =
'N'
, the function does not remove small columns of the scaled matrix. This option assumes that BLAS and QR factorizations and triangular solvers are implemented to work in that range. If the condition of
A
 if greater that
big
, use
?gesvj
.
If
jobr
 =
'R'
, restricted range for singular values of the scaled matrix
A
 is
[sqrt(
?lamch
(
'S'
), sqrt(big)]
, roughly as described above. This option is recommended.
For computing the singular values in the full range
[
?lamch
(
'S'
),big]
, use
?gesvj
.
jobt
Must be
'T'
 or
'N'
.
If the matrix is square, the procedure may determine to use a transposed
A
 if
A
T
 (for real flavors) or
A
H
 (for complex flavors) seems to be better with respect to convergence. If the matrix is not square,
jobt
 is ignored.
The decision is based on two values of entropy over the adjoint orbit of
A
T
 *
A
 (for real flavors) or
A
H
 *
A
 (for complex flavors). See the descriptions of
stat
[5]
 and
stat
[6]
.
If
jobt
 =
'T'
, the function performs transposition if the entropy test indicates possibly faster convergence of the Jacobi process, if
A
 is taken as input. If
A
 is replaced with
A
T
 or
A
H
, the row pivoting is included automatically.
If
jobt
 =
'N'
, the functions attempts no speculations. This option can be used to compute only the singular values, or the full SVD (
u
,
sigma
, and
v
). For only one set of singular vectors (
u
 or
v
), the caller should provide both
u
 and
v
, as one of the arrays is used as workspace if the matrix
A
 is transposed. The implementer can easily remove this constraint and make the code more complicated. See the descriptions of
u
 and
v
.
The
jobt
 =
'T'
 option is experimental and its effect might not be the same in subsequent releases. Consider using the
jobt
 =
'N'
 instead.
jobp
Must be
'P'
 or
'N'
.
Enables structured perturbations of denormalized numbers. This option should be active if the denormals are poorly implemented, causing slow computation, especially in cases of fast convergence. For details, see [Drmac08-1], [Drmac08-2] . For simplicity, such perturbations are included only when the full SVD or only the singular values are requested. You can add the perturbation for the cases of computing one set of singular vectors.
If
jobp
 =
'P'
, the function introduces perturbation.
If
jobp
 =
'N'
, the function introduces no perturbation.
m
The number of rows of the input matrix
A
;
m
≥
 0
.
n
The number of columns in the input matrix
A
;
m
≥
n
≥
 0.
a
,
u
,
v
Array
a
(size
lda
*
n
 for column major layout and
lda
*
m
 for row major layout)
 is an array containing the
m
-by-
n
 matrix
A
.
u
 is a workspace array,
its size for column major layout is
ldu
*
n
 for
jobu
='U' or 'W' and
ldu
*
m
 for
jobu
='F'; for row major layout its size is at least
ldu
*
m
. When
jobt
 = 'T' and
m
 =
n
,
u
 must be provided even though
jobu
 = 'N'.
v
 is