gesvdq(3)

📦
version 3.12.1 · 9 versions
, in arch rolling

💅 gesvdq(3) — Library Functions Manual

okay so bestie, you're reading gesvdq(3) out of the Library Functions Manual, and honestly? this one's a whole meal. it's the LAPACK squad that does your SVD and your QR with pivoting, and there are four flavors of it depending on whether your matrix is real or complex, single or double precision. no cap, it's giving "linear algebra final boss" energy 👑

✨ NAME ✨

gesvdq — gesvdq: SVD, QR with pivoting 🎀

🌈 SYNOPSIS 🌈

🫶 Functions

so there's literally four of them, one per data type, and they all do the same iconic thing. here's the lineup:

  • subroutine cgesvdq (joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork, rwork, lrwork, info) — CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices 💖
  • subroutine dgesvdq (joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work, lwork, rwork, lrwork, info) — DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices 💖
  • subroutine sgesvdq (joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work, lwork, rwork, lrwork, info) — SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices 💖
  • subroutine zgesvdq (joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork, rwork, lrwork, info) — ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices 💖

🦋 Detailed Description 🦋

🧚 Function Documentation

⭐ subroutine cgesvdq ⭐

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

CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices 🎀

🌸 Purpose 🌸

okay so, what does it actually do? here's the tea:

CGESVDQ 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.

periodt. that's the whole vibe 💅

👑 Parameters 👑

💜 JOBA

JOBA is CHARACTER*1 Specifies 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 = SLAMCH('Epsilon'). This authorises CGESVDQ 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(1). N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1) ⭐

💜 JOBP

JOBP is 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 🛡️
  • 'N' — no row pivoting 🙅

💜 JOBR

JOBR is CHARACTER*1
  • 'T' — after the initial pivoted QR factorization, CGESVD 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

JOBU is 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

JOBV is 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

M is INTEGER The number of rows of the input matrix A. M >= 0. 📏

💜 N

N is INTEGER The number of columns of the input matrix A. M >= N >= 0. 📏

💜 A

A is COMPLEX 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 CGEQP3. 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

LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M). 📐

💜 S

S is REAL array of dimension N. The singular values of A, ordered so that S(i) >= S(i+1). 🌈

💜 U

U is COMPLEX 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

LDU is 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

V is COMPLEX 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

LDV is INTEGER The leading dimension of the array V.
  • If JOBV = 'A', 'V', 'R', or JOBA = 'E' → LDV >= max(1,N).
  • Otherwise → LDV >= 1.

💜 NUMRANK

NUMRANK is 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 👀

💜 IWORK

IWORK is 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

LIWORK is 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

CWORK is a COMPLEX array, dimension (max(2, LCWORK)), and it's here to be your workspace bestie 💅

On exit, if — on entry — LCWORK.NE.-1, then CWORK(1:N) contains the parameters you need to recover the Q factor from the QR factorization that CGEQP3 computed. (iykyk 🎀)

And if LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0:

  • CWORK(1) returns the optimal LCWORK
  • CWORK(2) returns the minimal LCWORK

LCWORK — the size math, bestie 🔥

LCWORK is INTEGER. It's the dimension of the array CWORK, and it's determined as follows. First, 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 the minimal value of LCWORK is:

  • = 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.

Finally — and this is the floor, no negotiating — LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ). 🫶

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. (So don't go looking for a drag, there isn't one 💅)

RWORK 🌸

RWORK is REAL array, dimension (max(1, LRWORK)).

On exit:

  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.
  2. RWORK(2) contains the number of singular values computed as exact zeros in CGESVD applied to the upper triangular or trapezoidal R (from the initial QR factorization). In case of early exit (no call to CGESVD, such as in the case of zero matrix) RWORK(2) = -1.

If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, RWORK(1) returns the minimal LRWORK. ⭐

LRWORK 🦋

LRWORK is INTEGER. It's the dimension of the array RWORK.

  • If JOBP ='P', then LRWORK >= MAX(2, M, 5*N)
  • Otherwise, LRWORK >= MAX(2, 5*N)
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 here — not a typo of mine, that's literally what the source says 😭)

INFO 👑

INFO is INTEGER:

  • = 0: successful exit. periodt. ✨
  • < 0: if INFO = -i, the i-th argument had an illegal value. 🙈
  • > 0: if CBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B (computed in CGESVD) did not converge to zero. 💀

Further Details (the behind-the-scenes tea ☕)

  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.
  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 (yes, that's the spelling in the source, love that for them 🫶) 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. 🫶

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 cgesvdq.f.

💜 and now, the DGESVDQ girlies 💜

   subroutine dgesvdq (character joba, character jobp, character jobr,
       character jobu, character jobv, integer m, integer n, double precision,
       dimension( lda, * ) a, integer lda, double precision, dimension( * ) s,
       double precision, dimension( ldu, * ) u, integer ldu, double precision,
       dimension( ldv, * ) v, integer ldv, integer numrank, integer, dimension(
       * ) iwork, integer liwork, double precision, dimension( * ) work,
       integer lwork, double precision, dimension( * ) rwork, integer lrwork,
       integer info)

DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices. iconic behaviour, honestly.

Purpose: ✨

DGESVDQ computes the singular value decomposition (SVD) of a real 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 orthogonal 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

JOBA 🎯

JOBA is CHARACTER*1. It specifies 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 DGESVDQ 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(1). N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)

JOBP 🎀

JOBP is 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. 👑
  • = 'N' — No row pivoting.

JOBR

JOBR

  • = 'T' — so after the initial pivoted QR factorization, DGESVD gets applied to the transposed R**T of the computed triangular factor R. Heads up: that means some extra data movement (matrix transpositions, ugh). Lowkey only really useful for experiments, research and development. 🧪
  • = 'N' — the triangular factor R goes straight into DGESVD as input. This one might be preferred fr, since it involves less data movement. ✨

JOBU

JOBU is CHARACTER*1:

  • = 'A' — all M left singular vectors get computed and returned in the matrix U. Peep the description of U. 💅
  • = 'S' or 'U' — N = min(M,N) left singular vectors get computed and returned in the matrix U. Again, see the description of U.
  • = 'R' — the 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 come back in factored form, as the product of the Q factor from the initial QR factorization and the N left singular vectors of (R**T , 0)**T. Not me explaining this twice: if row pivoting is used, then the necessary info on the row pivoting is stored in IWORK(N+1:N+M-1). 📌
  • = 'N' — the left singular vectors are not computed. We're just not doing that today. 🙅‍♀️

JOBV

JOBV is CHARACTER*1:

  • = 'A', 'V' — all N right singular vectors are computed and returned in the matrix V.
  • = 'R' — the 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. Like, genuinely illegal, don't.
  • = 'N' — the right singular vectors are not computed.

M

M is INTEGER

The number of rows of the input matrix A. M >= 0. 📏

N

N is INTEGER

The number of columns of the input matrix A. M >= N >= 0. (Yes bestie, M >= N — that's the vibe check and A has to pass it.) ✅

A

A is DOUBLE PRECISION 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 DGEQP3. And if JOBU = 'F', those Householder vectors, together with WORK(1:N), can be used to restore the Q factors from the initial pivoted QR factorization of A. See the description of U. 👀

LDA

LDA is INTEGER.

The leading dimension of the array A. LDA >= max(1,M). 📐

S

S is DOUBLE PRECISION array of dimension N.

The singular values of A, ordered so that S(i) >= S(i+1). Descending order, the way we like it. 📉

U

U is DOUBLE PRECISION 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 an N x N orthogonal matrix that can be used to form the left singular vectors.

If JOBU = 'N', U is not referenced. Iconic behaviour honestly, it's just not invited. 🙈

LDU

LDU is INTEGER.

The leading dimension of the array U.

  • If JOBU = 'A', 'S', 'U', 'R', then LDU >= max(1,M).
  • If JOBU = 'F', then LDU >= max(1,N).
  • Otherwise, LDU >= 1.

V

V is DOUBLE PRECISION 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 orthogonal matrix V**T.
  • If JOBV = 'R', V contains the first NUMRANK rows of V**T (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

LDV is INTEGER

The leading dimension of the array V.

  • If JOBV = 'A', 'V', 'R', or JOBA = 'E', then LDV >= max(1,N).
  • Otherwise, LDV >= 1.

NUMRANK

NUMRANK is 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 DGESVD. And heads up: the final value of NUMRANK might be further reduced if some singular values are computed as zeros. It can shrink, basically. 😭

IWORK

IWORK is INTEGER array, dimension (max(1, LIWORK)).

On exit, IWORK(1:N) contains the 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'. Keep them, you'll want them.

If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, IWORK(1) returns the minimal LIWORK. ✨

LIWORK

LIWORK is INTEGER

The dimension of the array IWORK. So:

LIWORK >= N + M - 1,     if JOBP = 'P' and JOBA .NE. 'E';
LIWORK >= N              if JOBP = 'N' and JOBA .NE. 'E';
LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E';
LIWORK >= N + N          if JOBP = 'N' and JOBA = 'E'.

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

WORK

WORK is DOUBLE PRECISION array, dimension (max(2, LWORK)), used as a workspace.

On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters needed to recover the Q factor from the QR factorization computed by DGEQP3. 🛠️

If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0:

  • WORK(1) returns the optimal LWORK, and
  • WORK(2) returns the minimal LWORK.

LWORK

LWORK is INTEGER

The dimension of the array WORK. It is determined as follows. Prepare yourself bestie, this is the whole formula sheet: 📝

Let  LWQP3 = 3*N+1,  LWCON = 3*N, and let
LWORQ = { MAX( N, 1 ),  if JOBU = 'R', 'S', or 'U'
        { MAX( M, 1 ),  if JOBU = 'A'
LWSVD = MAX( 5*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
Then the minimal value of LWORK is:
= 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, LWORQ ) if the singular values and the left
                         singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) 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, LWORQ ) if the full SVD is requested with JOBV = 'R';
                         independent of JOBR;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
                         JOBV = 'R' and, also a scaled condition
                         estimate requested; independent of JOBR;
= MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
               full SVD is requested with JOBV = 'A' or 'V', and
               JOBR ='N'
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
               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, LWORQ ),
N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
               full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
               if the full SVD is requested with JOBV = 'A' or 'V', and
               JOBR ='T', and also a scaled condition number estimate
               requested.
Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).

If LWORK = -1, then a workspace query is assumed; the routine only calculates and returns the optimal and minimal sizes for the WORK, IWORK, and RWORK arrays, and no error message related to LWORK is issued by XERBLA. No cap, that's a free pass on the error message. 🦋

RWORK

(coming right up, keep scrolling bestie 💜)

  • RWORK — DOUBLE PRECISION array, dimension (max(1, LRWORK)). On exit, bestie:
    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.
    2. RWORK(2) contains the number of singular values computed as exact zeros in DGESVD applied to the upper triangular or trapezoidal R (from the initial QR factorization). In case of early exit (no call to DGESVD, such as in the case of zero matrix) RWORK(2) = -1.

If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, RWORK(1) returns the minimal LRWORK.

  • LRWORK — INTEGER. The dimension of the array RWORK. If JOBP ='P', then LRWORK >= MAX(2, M). Otherwise, LRWORK >= 2.

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

  • INFO — INTEGER.
    • = 0: successful exit.
    • < 0: if INFO = -i, the i-th argument had an illegal value.
    • > 0: if DBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B (computed in DGESVD) 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.
  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 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.

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 412 of file dgesvdq.f.

subroutine sgesvdq (character joba, character jobp, character jobr,
       character jobu, character jobv, integer m, integer n, real, dimension(
       lda, * ) a, integer lda, real, dimension( * ) s, real, dimension( ldu, *
       ) u, integer ldu, real, dimension( ldv, * ) v, integer ldv, integer
       numrank, integer, dimension( * ) iwork, integer liwork, real, dimension(
       * ) work, integer lwork, real, dimension( * ) rwork, integer lrwork,
       integer info)

SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices

Purpose 🎯

SGESVDQ computes the singular value decomposition (SVD) of a real 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 orthogonal 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 ✨

  • JOBA — CHARACTER*1. Specifies 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 = SLAMCH('Epsilon'). This authorises CGESVDQ 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(1). N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
  • 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.
    • = 'N': No row pivoting.
  • JOBR — CHARACTER*1.
    • = 'T': After the initial pivoted QR factorization, SGESVD is applied to the transposed R**T 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 SGESVD. 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**T , 0)**T. 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 — REAL 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 SGEQP3. If JOBU = 'F', these Householder vectors together with WORK(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 — REAL array of dimension N. The singular values of A, ordered so that S(i) >= S(i+1).
  • U — REAL 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 orthogonal 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 —
  • V — this one is a REAL 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 orthogonal matrix V**T; if JOBV = 'R', V contains the first NUMRANK rows of V**T (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, bestie. It's the leading dimension of the array V. If JOBV = 'A', 'V', 'R', or JOBA = 'E', then 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 SGESVD. The final value of NUMRANK might be further reduced if some singular values are computed as zeros.
  • 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, LWORK, 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' and JOBA .NE. 'E'; LIWORK >= N if JOBP = 'N' and JOBA .NE. 'E'; LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E'; LIWORK >= N + N if JOBP = 'N' and JOBA = 'E'. If LIWORK = -1, then a workspace query is assumed; the routine only calculates and returns the optimal and minimal sizes for the WORK, IWORK, and RWORK arrays, and no error message related to LWORK is issued by XERBLA.
  • WORK — REAL array, dimension (max(2, LWORK)), used as a workspace. On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters needed to recover the Q factor from the QR factorization computed by SGEQP3. If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, WORK(1) returns the optimal LWORK, and WORK(2) returns the minimal LWORK.
  • LWORK — INTEGER. It's the dimension of the array WORK, and it is determined as follows — buckle up, because this is the formula drop:
Let  LWQP3 = 3*N+1,  LWCON = 3*N, and let
LWORQ = { MAX( N, 1 ),  if JOBU = 'R', 'S', or 'U'
        { MAX( M, 1 ),  if JOBU = 'A'
LWSVD = MAX( 5*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
Then the minimal value of LWORK is:
= 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, LWORQ ) if the singular values and the left
                         singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) 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, LWORQ ) if the full SVD is requested with JOBV = 'R';
                         independent of JOBR;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
                         JOBV = 'R' and, also a scaled condition
                         estimate requested; independent of JOBR;
= MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
               full SVD is requested with JOBV = 'A' or 'V', and
               JOBR ='N'
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
               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, LWORQ ),
N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
               full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
               if the full SVD is requested with JOBV = 'A' or 'V', and
               JOBR ='T', and also a scaled condition number estimate
               requested.
Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).

If LWORK = -1, then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
  • RWORK — REAL array, dimension (max(1, LRWORK)). On exit:
    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.
    2. RWORK(2) contains the number of singular values computed as exact zeros in SGESVD applied to the upper triangular or trapezoidal R (from the initial QR factorization). In case of early exit (no call to SGESVD, such as in the case of zero matrix) RWORK(2) = -1. If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, RWORK(1) returns the minimal LRWORK.
  • LRWORK — INTEGER. The dimension of the array RWORK. If JOBP ='P', then LRWORK >= MAX(2, M). Otherwise, LRWORK >= 2. If LRWORK = -1, then a workspace query is assumed; the routine only calculates and returns the optimal and minimal sizes for the WORK, IWORK, and RWORK arrays, and no error message related to LWORK is issued by XERBLA.
  • INFO — INTEGER. = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if SBDSQR did not converge, INFO specifies how many superdiagonals of an intermediate bidiagonal form B (computed in SGESVD) 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.
  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 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. 💖

📚 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 412 of file sgesvdq.f.

   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

Purpose

okay bestie, the doc block literally says 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. iconic math behavior 💅

Parameters

  • JOBA — JOBA is CHARACTER*1. It specifies 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 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 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(1). N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
  • JOBP — JOBP is 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.
    • = 'N' No row pivoting.
  • JOBR — JOBR is 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 — JOBU is 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 — JOBV is 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 — M is INTEGER. The number of rows of the input matrix A. M >= 0.
  • N — N is INTEGER. The number of columns of the input matrix A. M >= N >= 0.
  • A — A is 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 — LDA is INTEGER. The leading dimension of the array A. LDA >= max(1,M).
  • S — S is DOUBLE PRECISION array of dimension N. The singular values of A, ordered so that S(i) >= S(i+1).
  • U — U is 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 — LDU is 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 — V is 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 — LDV is INTEGER. The leading dimension of the array V.
    • If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N).
    • Otherwise, LDV >= 1.
  • NUMRANK — NUMRANK is 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.
  • IWORK — IWORK is 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 — LIWORK is 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 — CWORK is 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.
  • LCWORK

LCWORK is INTEGER, and it’s literally the dimension of the array CWORK 💅. It gets determined like this:

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 the minimal value of LCWORK is:

= 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 estimate 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.

Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ).

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.

RWORK 💖

RWORK is a DOUBLE PRECISION array, dimension (max(1, LRWORK)).

On exit:

  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.
  2. 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.

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

LRWORK 🦋

LRWORK is INTEGER. It’s the dimension of the array RWORK.

If JOBP = 'P', then LRWORK >= MAX(2, M, 5*N); Otherwise, LRWORK >= MAX(2, 5*N).

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.

INFO 🎀

INFO is INTEGER.

  • = 0: successful exit.
  • < 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.
  2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause column norm overflow. This is the minimal precaution 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.

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 🌸

Generated automatically by Doxygen for LAPACK from the source code.

Final sign-off, bestie: this page is from LAPACK, version 3.12.0, manual page gesvdq(3) — and the footer doesn’t give a date, so I’m not making one up 😭💖

📦 where this page comes from