/home/abui.../zgesvdq.f(3)

๐Ÿ“ฆ version 3.12.0, in opensuse-leap 16.0

๐ŸŒธ NAME ๐ŸŒธ

okay so bestie, this is the zgesvdq.f(3) page from the Library Functions Manual, and it's literally the page for this file:

/home/abuild/rpmbuild/BUILD/lapack-static-3.12.0-build/lapack-3.12.0/SRC/zgesvdq.f

โœจ SYNOPSIS โœจ

๐ŸŽ€ Functions/Subroutines ๐ŸŽ€

  • subroutine ZGESVDQ (...) โ€” ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices. iconic behavior, honestly.

Here's the exact call shape, bestie:

subroutine ZGESVDQ (joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u,
    ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork, rwork, lrwork,
    info)

๐Ÿ‘‘ Function/Subroutine Documentation ๐Ÿ‘‘

๐Ÿฆ‹ subroutine ZGESVDQ ๐Ÿฆ‹

subroutine ZGESVDQ (character joba, character jobp, character jobr,
    character jobu, character jobv, integer m, integer n, complex*16,
    dimension( lda, * ) a, integer lda, double precision, dimension( * ) s,
    complex*16, dimension( ldu, * ) u, integer ldu, complex*16, dimension(
    ldv, * ) v, integer ldv, integer numrank, integer, dimension( * ) iwork,
    integer liwork, complex*16, dimension( * ) cwork, integer lcwork, double
    precision, dimension( * ) rwork, integer lrwork, integer info)

ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices. periodt. ๐Ÿ’…

๐Ÿ’– Purpose ๐Ÿ’–

ZCGESVDQ computes the singular value decomposition (SVD) of a complex M-by-N matrix A, where M >= N. The SVD of A is written as:

                                   [++]   [xx]   [x0]   [xx]
             A = U * SIGMA * V^*,  [++] = [xx] * [ox] * [xx]
                                   [++]   [xx]

where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal matrix, and V is an N-by-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. slay. ๐ŸŒˆ

๐Ÿงš Parameters ๐Ÿงš

  • JOBA โ€” CHARACTER*1. Specifies the level of accuracy in the computed SVD, bestie.
    • '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 ZGESVDQ 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 the 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 the 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(1). N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1). big brain accuracy vibes. ๐Ÿง 
  • JOBP โ€” CHARACTER*1.
    • 'P' โ€” The rows of A are ordered in decreasing order with respect to ||A(i,:)||_\infty. This enhances numerical accuracy at the cost of extra data movement. Recommended for numerical robustness, no cap.
    • 'N' โ€” No row pivoting.
  • JOBR โ€” CHARACTER*1.
    • 'T' โ€” After the initial pivoted QR factorization, ZGESVD is applied to the adjoint 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 CGESVD. This may be preferred as it involves less data movement.
  • JOBU โ€” CHARACTER*1.
    • 'A' โ€” All M left singular vectors are computed and returned in the matrix U. See the description of U.
    • 'S' or 'U' โ€” N = min(M,N) left singular vectors are computed and returned in the matrix U. See the description of U.
    • 'R' โ€” Numerical rank NUMRANK is determined and only NUMRANK left singular vectors are computed and returned in the matrix 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+1:N+M-1).
    • 'N' โ€” The left singular vectors are not computed.
  • JOBV โ€” CHARACTER*1.
    • 'A', 'V' โ€” All N right singular vectors are computed and returned in the matrix V.
    • 'R' โ€” Numerical rank NUMRANK is determined and only NUMRANK right singular vectors are computed and returned in the matrix V. This option is allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal. ๐Ÿšจ
    • 'N' โ€” The right singular vectors are not computed.
  • M โ€” INTEGER. The number of rows of the input matrix A. M >= 0.
  • N โ€” INTEGER. The number of columns of the input matrix A. M >= N >= 0.
  • A โ€” COMPLEX*16 array of dimensions LDA x N. On entry, the input matrix A. On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains the Householder vectors as stored by ZGEQP3. If JOBU = 'F', these Householder vectors together with CWORK(1:N) can be used to restore the Q factors from the initial pivoted QR factorization of A. See the description of U.
  • LDA โ€” INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • S โ€” DOUBLE PRECISION array of dimension N. The singular values of A, ordered so that S(i) >= S(i+1). descending order, as she should. ๐Ÿ’ซ
  • U โ€” COMPLEX*16 array, dimension:
    • LDU x M if JOBU = 'A'; see the description of LDU. In this case, on exit, U contains the M left singular vectors.
    • LDU x N if JOBU = 'S', 'U', 'R'; see the description of LDU. In this case, U contains the leading N or the leading NUMRANK left singular vectors.
    • LDU x N if JOBU = 'F'; see the description of LDU. In this case U contains N x N unitary matrix that can be used to form the left singular vectors.
    • If JOBU = 'N', U is not referenced.
  • LDU โ€” INTEGER. The leading dimension of the array U.
    • If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M).
    • If JOBU = 'F', LDU >= max(1,N).
    • Otherwise, LDU >= 1.
  • V โ€” COMPLEX*16 array, dimension:
    • LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E'.
    • If JOBV = 'A', or 'V', V contains the N-by-N unitary matrix V**H.
    • If JOBV = 'R', V contains the first NUMRANK rows of V**H (the right singular vectors, stored rowwise, of the NUMRANK largest singular values).
    • If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
    • If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
  • LDV โ€” INTEGER. The leading dimension of the array V.
    • If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N).
    • Otherwise, LDV >= 1.
  • NUMRANK โ€” INTEGER. NUMRANK is the 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 CGESVD. The final value of NUMRANK might be further reduced if some singular values are computed as zeros. rank reveal tea. โ˜•
  • IWORK โ€” INTEGER array, dimension (max(1, LIWORK)). On exit, IWORK(1:N) contains column pivoting permutation of the rank revealing QR factorization. If JOBP = 'P', IWORK(N+1:N+M-1) 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 LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, IWORK(1) returns the minimal LIWORK.

  • LIWORK โ€” INTEGER. The dimension of the array IWORK.
    • LIWORK >= N + M - 1, if JOBP = 'P';
    • LIWORK >= N if JOBP = 'N'.

If LIWORK = -1, then a workspace query is assumed; the routine only calculates and returns the optimal and minimal sizes for the CWORK, IWORK, and RWORK arrays, and no error message related to LCWORK is issued by XERBLA.

  • CWORK โ€” COMPLEX*12 array, dimension (max(2, LCWORK)), used as a workspace. On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters needed to recover the Q factor from the QR factorization computed by ZGEQP3.

If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, CWORK(1) returns the optimal LCWORK, and CWORK(2) returns the minimal LCWORK.

yes bestie, it says COMPLEX*12 right there โ€” I'm keeping it exactly as the page has it. ๐Ÿ’…
  • LCWORK โ€” this is exactly where part 1 cuts off, bestie; part 2 picks it up right here. ๐Ÿซถ

๐Ÿ’พ LCWORK

ok so bestie, LCWORK is an INTEGER and it's literally just how much room you're giving the CWORK array. and the way you figure out the minimum? buckle up ๐ŸŽ€ โ€” it's a whole little recipe:

Let  LWQP3 = N+1,  LWCON = 2*N, and let
LWUNQ = { MAX( N, 1 ),  if JOBU = 'R', 'S', or 'U'
        { MAX( M, 1 ),  if JOBU = 'A'
LWSVD = MAX( 3*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )

then, with all those little numbers in your pocket, the minimal LCWORK works out to:

= 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', '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', 'V' and
               JOBR ='T', and also a scaled condition number estimate
               requested.

and then โ€” the final plot twist ๐Ÿ˜ญ โ€” LCWORK has to be at least two: LCWORK = MAX( 2, LCWORK ).

๐Ÿšจ but wait, there's a special case: if LCWORK = -1, then a workspace query is assumed. the routine only calculates and returns the optimal and minimal sizes for the CWORK, IWORK, and RWORK arrays, and no error message related to LCWORK is issued by XERBLA. it's giving โœจwindow shopping, not buyingโœจ

๐Ÿงฎ RWORK

RWORK is a DOUBLE PRECISION array, dimension (max(1, LRWORK)). and on exit? two little gifts ๐Ÿ’œ

  1. If JOBA = 'E', RWORK(1) 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(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).

Otherwise, RWORK(1) = -1. (we love a graceful little dash-one, ngl ๐Ÿ™ˆ)

  1. RWORK(2) 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(2) = -1. ๐Ÿ’€

and, quote, if LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, RWORK(1) returns the minimal LRWORK. cute. iconic. ๐Ÿซถ

๐Ÿ“ LRWORK

LRWORK is an INTEGER. it's the dimension of the array RWORK, bestie. so:

  • If JOBP ='P', then LRWORK >= MAX(2, M, 5*N); ๐ŸŒธ
  • Otherwise, LRWORK >= MAX(2, 5*N).
๐Ÿ’… same drill as before: if LRWORK = -1, then a workspace query is assumed โ€” the routine only calculates and returns the optimal and minimal sizes for the CWORK, IWORK, and RWORK arrays, and no error message related to LCWORK is issued by XERBLA. (yes, it says LCWORK there, we're just reading the page as written, don't come for me ๐Ÿ˜ญ)

๐Ÿšฆ INFO

INFO is an INTEGER. the vibe check at the end:

  • = 0: successful exit. slay ๐Ÿ”ฅ
  • < 0: if INFO = -i, the i-th argument had an illegal value. ๐Ÿ™ˆ
  • > 0: if ZBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B (computed in ZGESVD) did not converge to zero. ๐Ÿซ 

๐Ÿ” Further Details

  1. The data movement (matrix transpose) is coded using simple nested DO-loops because BLAS and LAPACK do not provide corresponding subroutines. Those DO-loops are easily identified in this source code โ€” by the CONTINUE statements labeled with 11**. In an optimized version of this code, the nested DO loops should be replaced with calls to an optimized subroutine. iykyk ๐Ÿงš
  2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause column norm overflow. This is the minial precaution (their spelling, not mine, love u Drmaฤ ๐Ÿ’–) and it is left to the SVD routine (CGESVD) to do its own preemptive scaling if potential over- or underflows are detected. To avoid repeated scanning of the array A, an optimal implementation would do all necessary scaling before calling CGESVD and the scaling in CGESVD can be switched off.
  3. Other comments related to code optimization are given in comments in the code, enclosed in [[double brackets]]. ๐ŸŽ€

๐Ÿ› Bugs, examples and comments

Please report all bugs and send interesting examples and/or comments to drmac@math.hr. Thank you. โ€” literally a sweet little "ty" note, we stan ๐Ÿซถ

๐Ÿ“š References

  • [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for Computing the SVD with High Accuracy. ACM Trans. Math. Softw. 44(1): 11:1-11:30 (2017)
  • SIGMA library, xGESVDQ section updated February 2016. Developed and coded by Zlatko Drmac, Department of Mathematics, University of Zagreb, Croatia, drmac@math.hr

๐ŸŒŸ Contributors:

Developed and coded by Zlatko Drmac, Department of Mathematics, University of Zagreb, Croatia, drmac@math.hr ๐Ÿ‘‘

โœ๏ธ Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 410 of file zgesvdq.f.

๐ŸŽฌ Author (again, yes, there are two of them)

Generated automatically by Doxygen for LAPACK from the source code.

and that's the whole page, bestie โ€” go forth and query your workspaces responsibly ๐Ÿ’…โœจ

๐Ÿ“ source: zgesvdq.f ยท LAPACK Version 3.12.0 ยท from the LAPACK source tree ๐Ÿ“ ily, bye ๐Ÿฆ‹

๐Ÿ“ฆ where this page comes from