?gejsv
?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
A
= U
*Σ
*V
H
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.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.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.Optimization Notice
|
---|
Intel's compilers may or may not optimize to the same degree for non-Intel microprocessors for optimizations that are not unique to Intel microprocessors. These optimizations include SSE2, SSE3, and SSSE3 instruction sets and other optimizations. Intel does not guarantee the availability, functionality, or effectiveness of any optimization on microprocessors not manufactured by Intel. Microprocessor-dependent optimizations in this product are intended for use with Intel microprocessors. Certain optimizations not specific to Intel microarchitecture are reserved for Intel microprocessors. Please refer to the applicable product User and Reference Guides for more information regarding the specific instruction sets covered by this notice.
Notice revision #20110804
|
This notice covers the following instruction sets: SSE2, SSE4.2, AVX2, AVX-512.
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, high relative accuracy is achieved ifjoba='C'with well-conditionedA=B*DBand arbitrary diagonal matrixD. The accuracy cannot be spoiled by column scaling. The accuracy of the computed output depends on the condition ofB, and the procedure aims at the best theoretical accuracy. The relative errormax_{i=1:N}|d sigma_i| / sigma_iis bounded byf(M,N)*epsilon* cond(B), independent ofD. 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 ofjoba. Additional actions are specified as follows:If, computation as withjoba='E''C'with an additional estimate of the condition number ofB. It provides a realistic error bound.If, accuracy higher than in thejoba='F''C'option is achieved, ifwith ill-conditioned diagonal scalingsA=D1*C*D2D1,D2, and a well-conditioned matrixC. This option is advisable, if the structure of the input matrix is not known and relative accuracy is desirable. The input matrixAis preprocessed with QR factorization with full (row and column) pivoting.If, computation as withjoba='G''F'with an additional estimate of the condition number ofB, where. IfA=B*DAhas heavily weighted rows, using this condition number gives too pessimistic error bound.If, small singular values are the noise and the matrix is treated as numerically rank defficient. The error in the computed singular values is bounded byjoba='A'f(m,n)*epsilon*||A||. The computed SVDA = U*S*V**t(for real flavors) orA = U*S*V**H(for complex flavors) restoresAup tof(m,n)*epsilon*||A||. This enables the procedure to set all singular values belown*epsilon*||A||to zero.If, the procedure is similar to thejoba='R''A'option. Rank revealing property of the initial QR factorization is used to reveal (using triangular factor) a gapsigma_{r+1}, in which case the numerical rank is declared to be<epsilon * sigma_rr. 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 matrixU:If,jobu='U'ncolumns ofUare returned in the arrayuIf, a full set ofjobu='F'mleft singular vectors is returned in the arrayu.If,jobu='W'umay be used as workspace of lengthm*n. See the description ofu.If,jobu='N'uis not computed.
- jobv
- Must be'V','J','W', or'N'.Specifies whether to compute the matrixV:If,jobv='V'ncolumns ofVare returned in the arrayv; Jacobi rotations are not explicitly accumulated.If,jobv='J'ncolumns ofVare returned in the arrayvbut they are computed as the product of Jacobi rotations. This option is allowed only ifjobu≠'N'If,jobv='W'vmay be used as workspace of lengthn*n. See the description ofv.If,jobv='N'vis 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. IfAis scaled so that the largest singular value of the scaled matrix is aroundsqrt(big), big =, the function can remove columns of?lamch('O')Awhose norm in the scaled matrix is less thansqrt((for?lamch('S'))jobr='R'), or less thansmall =.?lamch('S')/?lamch('E')If, 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 ofjobr='N'Aif greater thatbig, use?gesvj.If, restricted range for singular values of the scaled matrixjobr='R'Ais[sqrt(, roughly as described above. This option is recommended.?lamch('S'), sqrt(big)]For computing the singular values in the full range[, use?lamch('S'),big]?gesvj.
- jobt
- Must be'T'or'N'.If the matrix is square, the procedure may determine to use a transposedAif(for real flavors) orAT(for complex flavors) seems to be better with respect to convergence. If the matrix is not square,AHjobtis ignored.The decision is based on two values of entropy over the adjoint orbit of*ATA(for real flavors) or*AHA(for complex flavors). See the descriptions of.andstat[5]stat[6]If, the function performs transposition if the entropy test indicates possibly faster convergence of the Jacobi process, ifjobt='T'Ais taken as input. IfAis replaced withorAT, the row pivoting is included automatically.AHIf, the functions attempts no speculations. This option can be used to compute only the singular values, or the full SVD (jobt='N'u,sigma, andv). For only one set of singular vectors (uorv), the caller should provide bothuandv, as one of the arrays is used as workspace if the matrixAis transposed. The implementer can easily remove this constraint and make the code more complicated. See the descriptions ofuandv.Theoption is experimental and its effect might not be the same in subsequent releases. Consider using thejobt='T'instead.jobt='N'
- 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, the function introduces perturbation.jobp='P'If, the function introduces no perturbation.jobp='N'
- m
- The number of rows of the input matrixA;.m≥0
- n
- The number of columns in the input matrixA;m≥n≥0.
- a,u,v
- Arraya(sizeis an array containing thelda*nfor column major layout andlda*mfor row major layout)m-by-nmatrixA.uis a workspace array,its size for column major layout is. Whenldu*nforjobu='U' or 'W' andldu*mforjobu='F'; for row major layout its size is at leastldu*mjobt= 'T' andm=n,umust be provided even thoughjobu= 'N'.vis a workspace array, its size is. Whenldv*njobt= 'T' andm=n,vmust be provided even thoughjobv= 'N'.
- lda
- The leading dimension of the arraya. Must be at leastmax(1,m)for column major layout and at least max(1,.n) for row major layout
- sva
- svais a workspace array, its size isn.
- ldu
- The leading dimension of the arrayu;.ldu≥1jobu='U'or'F'or'W',ldu≥mfor column major layout; for row major layout ifjobu='U'or