XLPack 7.0
XLPack Numerical Library (C API) Reference Manual
Loading...
Searching...
No Matches

◆ zgesvdq()

void zgesvdq ( char  joba,
char  jobp,
char  jobr,
char  jobu,
char  jobv,
int  m,
int  n,
int  lda,
doublecomplex  a[],
double  s[],
int  ldu,
doublecomplex  u[],
int  ldv,
doublecomplex  v[],
int *  numrank,
doublecomplex  work[],
int  lwork,
double  rwork[],
int  lrwork,
int  iwork[],
int  liwork,
int *  info 
)

Singular value decomposition (SVD) of a complex matrix (Preconditioned QR method)

Purpose
This routine computes the singular value decomposition (SVD) of a complex m x n matrix A, where m >= n. The SVD of A is written as
[++] [xx] [x0] [xx]
A = U * SIGMA * V^H, [++] = [xx] * [ox] * [xx]
[++] [xx]
where SIGMA is an n x n diagonal matrix, U is an m x n orthonormal matrix, and V is an n x n unitary matrix. The diagonal elements of SIGMA are the singular values of A. The columns of U and V are the left and the right singular vectors of A, respectively.
Parameters
[in]jobaSpecifies the level of accuracy in the computed SVD.
= 'A': The requested accuracy corresponds to having the backward error bounded by || delta A ||_F <= f(m, n) * EPS * || A ||_F, where EPS = dlamch('Epsilon'). This authorises this routine to truncate the computed triangular factor in a rank revealing QR factorization whenever the truncated part is below the threshold of the order of EPS * ||A||_F. This is aggressive truncation level.
= 'M': Similarly as with 'A', but the truncation is more gentle: it is allowed only when there is a drop on the diagonal of the triangular factor in the QR factorization. This is medium truncation level.
= 'H': High accuracy requested. No numerical rank determination based on the rank revealing QR factorization is attempted.
= 'E': Same as 'H', and in addition the condition number of column scaled A is estimated and returned in rwork[0]. n^(-1/4)*rwork[0] <= ||pinv(A_scaled)||_2 <= n^(1/4)*rwork[0].
[in]jobp= 'P': The rows of A are ordered in decreasing order with respect to ||A(i,:)||_infinity. This enhances numerical accuracy at the cost of extra data movement. Recommended for numerical robustness.
= 'N': No row pivoting.
[in]jobr= 'T': After the initial pivoted QR factorization, zgesvd is applied to the transposed R^H of the computed triangular factor R. This involves some extra data movement (matrix transpositions). Useful for experiments, research and development.
= 'N': The triangular factor R is given as input to zgesvd. This may be preferred as it involves less data movement.
[in]jobu= 'A': All m left singular vectors are computed and returned in the array u[][]. See the description of u[][].
= 'S' or 'U': min(m, n) (= n) left singular vectors are computed and returned in the array u[][]. See the description of u[][].
= 'R': Numerical rank numrank is determined and only numrank left singular vectors are computed and returned in the array u[][].
= 'F': The n left singular vectors are returned in factored form as the product of the Q factor from the initial QR factorization and the n left singular vectors of (R^H, 0)^H. If row pivoting is used, then the necessary information on the row pivoting is stored in iwork[n:n+m-2].
= 'N': The left singular vectors are not computed.
[in]jobv= 'A' or 'V': All n right singular vectors are computed and returned in the array v[][].

= 'R': Numerical rank numrank is determined and only numrank right singular vectors are computed and returned in the array v[][]. This option is allowed only if jobu = 'R' or jobu = 'N'. Otherwise it is illegal.
= 'N': The right singular vectors are not computed.
[in]mNumber of rows of the input matrix A. (m >= 0) (If m = 0, returns without computation)
[in]nNumber of columns of the input matrix A. (m >= n >= 0) (If n = 0, returns without computation)
[in]ldaLeading dimension of the two dimensional array a[][]. (lda >= max(1, m))
[in,out]a[][]Array a[la][lda] (la >= n)
[in] The input matrix A.
[out] If jobu != 'N' or jobv != 'N', the lower triangle of a[][] contains the Householder vectors as stored by zgeqp3. If jobu = 'F', these Householder vectors together with work[0:n-1] can be used to restore the Q factors from the initial pivoted QR factorization of A. See the description of u[][].
[out]s[]Array s[ls] (ls >= n)
The singular values of A, ordered so that s[i] >= s[i+1].
[in]lduLeading dimension of the two dimensional array u[][]. (ldu >= max(1, m) if jobu = 'A', 'S', 'U' or 'R'. ldu >= max(1, n) if jobu = 'F'. ldu >= 1 otherwise.)
[out]u[][]Array u[lu][ldu] (lu >= m if jobu = 'A'. lu >= n if jobu = 'S', 'U', 'R' or 'F'.)
jobu = 'A': u[][] contains the m left singular vectors.
jobu = 'S', 'U' or 'R': u[][] contains the leading n or the leading numrank left singular vectors.
jobu = 'F': u[][] contains the n x n unitary matrix that can be used to form the left singular vectors.
jobu = 'N': u[][] is not referenced.
[in]ldvLeading dimension of the two dimensional array v[][]. (ldv >= max(1, n) if jobv = 'A', 'V' or 'R', or joba = 'E'. ldv >= 1 otherwise.)
[out]v[][]Array v[lv][ldv] (lv >= n if jobv = 'A', 'V' or 'R', or if joba = 'E')
jobv = 'A' or 'V': v[][] contains the n x n unitary matrix V^H.
jobv = 'R': v[][] contains the first numrank rows of V^H (the right singular vectors, stored rowwise, of the numrank largest singular values).
jobv = 'N' and joba = 'E': v[][] is used as a workspace.
jobv = 'N' and joba != 'E': v[][] is not referenced.
[out]numrankThe numerical rank first determined after the rank revealing QR factorization, following the strategy specified by the value of joba. if jobv = 'R' and jobu = 'R', only numrank leading singular values and vectors are then requested in the call of zgesvd. The final value of numrank might be further reduced if some singular values are computed as zeros.
[out]work[]Array work[lwork]
Work array.
On exit, if, on entry, lwork != -1, work[0:n-1] contains parameters needed to recover the Q factor from the QR factorization computed by zgeqp3.

If lwork, lrwork, or liwork = -1, then on exit, if info = 0, work[0] returns the optimal lwork, and work[1] returns the minimal lwork.
[in]lworkThe dimension of the array work[]. (lwork must be at least 2: lwork = max(2, lwork))

Let
lwqp3 = n + 1;
lwcon = 2*n;
lwunq = max(1, n) if jobu = 'R', 'S', or 'U', lwunq = max(1, m) if jobu = 'A';
lwsvd = max(1, 3*n);
lwlqf = max(1, n/2);
lwsvd2 = max(1, 3*(n/2));
lwunlq = max(1, n);
lwqrf = max(1, n/2); and
lwunq2 = max(1, n).

Then the minimal value of lwork
= max(n + lwqp3, lwsvd) if only the singular values are needed;
= max(n + lwqp3, lwcon, lwsvd) if only the singular values are needed and a scaled condition estimate requested;
= n + max(lwqp3, lwsvd, lwunq) if the singular values and the left singular vectors are requested;
= n + max(lwqp3, lwcon, lwsvd, lwunq) if the singular values and the left singular vectors are requested, and also a scaled condition estimate requested;
= n + max(lwqp3, lwsvd) if the singular values and the right singular vectors are requested;
= n + max(lwqp3, lwcon, lwsvd) if the singular values and the right singular vectors are requested, and also a scaled condition etimate requested;
= n + max(lwqp3, lwsvd, lwunq) if the full SVD is requested with jobv = 'R' independent of jobr;
= n + max(lwqp3, lwcon, lwsvd, lwunq) if the full SVD is requested, jobv = 'R' and, also a scaled condition estimate requested independent of jobr;
= max(n + max(lwqp3, lwsvd, lwunq), n + max(lwqp3, n/2+lwlqf, n/2+lwsvd2, n/2+lwunlq, lwunq)) if the full SVD is requested with jobv = 'A' or 'V', and jobr ='N';
= max(n + max(lwqp3, lwcon, lwsvd, lwunq), n + max(lwqp3, lwcon, n/2+lwlqf, n/2+lwsvd2, n/2+lwunlq, lwunq)) if the full SVD is requested with jobv = 'A' or 'V', and jobr ='N', and also a scaled condition number estimate requested;
= max(n + max(lwqp3, lwsvd, lwunq), n + max(lwqp3, n/2+lwqrf, n/2+lwsvd2, n/2+lwunq2, lwunq)) if the full SVD is requested with jobv = 'A' or 'V', and jobr ='T';
= max(n + max(lwqp3, lwcon, lwsvd, lwunq), n + max(lwqp3, lwcon, n/2+lwqrf, n/2+lwsvd2, n/2+lwunq2, lwunq)) if the full SVD is requested with jobv = 'A' or 'V', and jobr ='T', and also a scaled condition number estimate requested.

If lwork = -1, then a workspace query is assumed. The routine only calculates and returns the optimal and minimal sizes for the work, rwork, and iwork arrays.
[out]rwork[]Array rwork[lrwork]
On exit,
  1. If joba = 'E', rwork[0] contains an estimate of the condition number of column scaled A. If A = C * D where D is diagonal and C has unit columns in the Euclidean norm, then, assuming full column rank, n^(-1/4) * rwork[0] <= ||pinv(c)||_2 <= n^(1/4) * rwork[0]. Otherwise, rwork[0] = -1.
  2. rwork[1] contains the number of singular values computed as exact zeros in zgesvd applied to the upper triangular or trapezoidal R (from the initial QR factorization). In case of early exit (no call to zgesvd, such as in the case of zero matrix) rwork[1] = -1.

    If lwork, lrwork, or liwork = -1, then on exit, if info = 0, rwork[0] returns the minimal lrwork.
[in]lrworkThe dimension of the array rwork[]. (lrwork >= max(2, m, 5*n) if jobp = 'P'. lrwork >= max(2, 5*n) otherwise.)

If lrwork = -1, then a workspace query is assumed. The routine only calculates and returns the optimal and minimal sizes for the work, rwork, and iwork arrays.
[out]iwork[]Array iwork[liwork]
On exit, iwork[0:n-1] contains column pivoting permutation of the rank revealing QR factorization.
If jobp = 'P', iwork[n:n+m-2] contains the indices of the sequence of row swaps used in row pivoting. These can be used to restore the left singular vectors in the case jobu = 'F'.

If lwork, lrwork, or liwork = -1, then on exit, if info = 0, iwork[0] returns the minimal liwork.
[in]liworkThe dimension of the array iwork[]. (liwork must be at least 1: liwork = max(1, liwork))
liwork >= n + m - 1 if jobp = 'P'; and
liwork >= n if jobp = 'N'.

If lwork = -1, then a workspace query is assumed. The routine only calculates and returns the optimal and minimal sizes for the work, rwork, and iwork arrays.
[out]info= 0: Successful exit
= -1: The argument joba had an illegal value (joba != 'A', 'M', 'H' nor 'E')
= -2: The argument jobp had an illegal value (jobp != 'P' nor 'N')
= -3: The argument jobr had an illegal value (jobr != 'T' nor 'N')
= -4: The argument jobu had an illegal value (jobu != 'A', 'S', 'R', 'F' nor 'N')
= -5: The argument jobv had an illegal value (jobv != 'A', 'V', 'R' nor 'N')
= -6: The argument m had an illegal value (m < 0)
= -7: The argument n had an illegal value (n < 0 or n > m)
= -8: The argument lda had an illegal value (lda < max(1, m))
= -9: The argument a had an illegal value (a[][] had a NAN entry)
= -11: The argument ldu had an illegal value (ldu too small)
= -13: The argument ldv had an illegal value (ldv too small)
= -17: The argument lwork had an illegal value (lwork too small)
= -19: The argument lrwork had an illegal value (lrwork too small)
= -21: The argument liwork had an illegal value (liwork too small)
> 0: zbdsqr did not converge. info specifies how many superdiagonals of an intermediate bidiagonal form B (computed in zgesvd) did not converge to zero.
Reference
LAPACK