gesvj(3)

📦

version 3.12.1 · 9 versions

3.12.1 · arch rolling (latest)
3.12.0 · opensuse-leap 16.0
3.12.1 · ubuntu stonking
3.12.1 · ubuntu questing
3.12.0 · ubuntu noble
3.12.0 · ubuntu noble
3.12.0 · centos-stream 10, fedora 43, fedora 44
3.12.1 · ubuntu resolute
3.12.1 · debian trixie

🌸 NAME 🌸

gesvj — gesvj: SVD, Jacobi, low-level. This is gesvj(3) from the Library Functions Manual, and it's the low-level Jacobi SVD page. 💅

✨ SYNOPSIS ✨

okay so here are the four function signatures, exactly as they appear, no edits bestie:

🦋 Functions 🦋

subroutine cgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv,
    cwork, lwork, rwork, lrwork, info)
     CGESVJ
subroutine dgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv,
    work, lwork, info)
    DGESVJ
subroutine sgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv,
    work, lwork, info)
    SGESVJ
subroutine zgesvj (joba, jobu, jobv, m, n, a, lda, sva, mv, v, ldv,
    cwork, lwork, rwork, lrwork, info)
     ZGESVJ

💫 DETAILED DESCRIPTION 💫

This page skips straight to the function docs, so let's go. 👑

👑 Function Documentation 👑

🌟 `CGESVJ` 🌟

subroutine cgesvj (character*1 joba, character*1 jobu, character*1 jobv,
    integer m, integer n, complex, dimension( lda, * ) a, integer lda, real,
    dimension( n ) sva, integer mv, complex, dimension( ldv, * ) v, integer
    ldv, complex, dimension( lwork ) cwork, integer lwork, real, dimension(
    lrwork ) rwork, integer lrwork, integer info)
     CGESVJ

Purpose:

CGESVJ 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. 💖

Parameters

🌸 `JOBA` 🌸

JOBA is CHARACTER*1. Specifies the structure of A.

= 'L': the input matrix A is lower triangular;
= 'U': the input matrix A is upper triangular;
= 'G': the input matrix A is general M-by-N matrix, M >= N.

🦋 `JOBU` 🦋

JOBU is CHARACTER*1. Specifies whether to compute the left singular vectors (columns of U):

= 'U' or 'F': the left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A. See more details in the description of A. The default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E').
= 'C': analogous to JOBU='U', except that the user can control the level of numerical orthogonality of the computed left singular vectors. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK. No CTOL smaller than ONE is allowed. CTOL greater than 1 / EPS is meaningless. The option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations. See the descriptions of A and WORK(1).
= 'N': the matrix U is not computed. However, see the description of A.

🎀 `JOBV` 🎀

JOBV is CHARACTER*1. Specifies whether to compute the right singular vectors, that is, the matrix V:

= 'V' or 'J': the matrix V is computed and returned in the array V;
= 'A': the Jacobi rotations are applied to the MV-by-N array V. In other words, the right singular vector matrix V is not computed explicitly; instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V.
= 'N': the matrix V is not computed and the array V is not referenced.

🐣 `M` 🐣

M is INTEGER. The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.

🐣 `N` 🐣

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

💜 `A` 💜

A is COMPLEX array, dimension (LDA,N). On entry, the M-by-N matrix A.

On exit,

If JOBU = 'U' .OR. JOBU = 'C':
- If INFO = 0:
  - RANKA orthonormal columns of U are returned in the leading RANKA columns of the array A. Here RANKA <= N is the number of computed singular values of A that are above the underflow threshold SLAMCH('S'). The singular vectors corresponding to underflowed or zero singular values are not computed. The value of RANKA is returned in the array RWORK as RANKA=NINT(RWORK(2)). Also see the descriptions of SVA and RWORK. The computed columns of U are mutually numerically orthogonal up to approximately TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), see the description of JOBU.
  - If INFO > 0, the procedure CGESVJ did not converge in the given number of iterations (sweeps). In that case, the computed columns of U may not be orthogonal up to TOL. The output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small.
If JOBU = 'N':
- If INFO = 0:
  - Note that the left singular vectors are 'for free' in the one-sided Jacobi SVD algorithm. However, if only the singular values are needed, the level of numerical orthogonality of U is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS. Thus, on exit, A contains the columns of U scaled with the corresponding singular values.
- If INFO > 0:
  - the procedure CGESVJ did not converge in the given number of iterations (sweeps).

🧚 `LDA` 🧚

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

🍓 `SVA` 🍓

SVA is REAL array, dimension (N).

On exit,

If INFO = 0:
- depending on the value SCALE = RWORK(1), we have:
  - If SCALE = ONE: SVA(1:N) contains the computed singular values of A. During the computation SVA contains the Euclidean column norms of the iterated matrices in the array A.
  - If SCALE .NE. ONE: the singular values of A are SCALE*SVA(1:N), and this factored representation is due to the fact that some of the singular values of A might underflow or overflow.
If INFO > 0: the procedure CGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.

🌈 `MV` 🌈

MV is INTEGER. If JOBV = 'A', then the product of Jacobi rotations in CGESVJ is applied to the first MV rows of V. See the description of JOBV.

🦄 `V` 🦄

V is COMPLEX array, dimension (LDV,N).

If JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors;
If JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V.
If JOBV = 'N', then V is not referenced.

🌟 `LDV` 🌟

LDV is INTEGER. The leading dimension of the array V, LDV >= 1.

If JOBV = 'V', then LDV >= max(1,N).
If JOBV = 'A', then LDV >= max(1,MV).

💠 `CWORK` 💠

CWORK is COMPLEX array, dimension (max(1,LWORK)). Used as workspace.

🔥 `LWORK` 🔥

LWORK is INTEGER. Length of CWORK. LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M+N, otherwise.

If on entry LWORK = -1, then a workspace query is assumed and no computation is done; CWORK(1) is set to the minial (and optimal) length of CWORK.

🎀 `RWORK` 🎀

RWORK — REAL array, dimension (max(6,LRWORK)). Okay so this one is doing double duty depending on whether you're entering or exiting, bestie.

On entry, if JOBU = 'C':

RWORK(1) = CTOL, where CTOL defines the threshold for convergence. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=SLAMCH('E').

🚨 It is required that CTOL >= ONE, i.e. it is not allowed to force the routine to obtain orthogonality below EPSILON.

On exit, slay, here's the receipts:

RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A. (See description of SVA().)
RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero singular values.
RWORK(3) = NINT(RWORK(3)) is the number of the computed singular values that are larger than the underflow threshold.
RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence.
RWORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. This is useful information in cases when CGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis.
RWORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep. It can be useful for a post festum analysis.

LRWORK — INTEGER. Length of RWORK. LRWORK >= 1, if MIN(M,N) = 0, and LRWORK >= MAX(6,N), otherwise.

If on entry LRWORK = -1, then a workspace query is assumed and no computation is done; RWORK(1) is set to the minial (and optimal) length of RWORK. Yes, minial — the typo stays, it's part of the lore 💅

INFO — INTEGER.

= 0: successful exit.
< 0: if INFO = -i, then the i-th argument had an illegal value.
> 0: CGESVJ did not converge in the maximal allowed number (NSWEEP=30) of sweeps. The output may still be useful. See the description of RWORK.

👑 Author

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

🧚 Further Details

Okay so the orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations. In the case of underflow of the tangent of the Jacobi angle, a modified Jacobi transformation of Drmac [3] is used. Pivot strategy uses column interchanges of de Rijk [1]. The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [2]. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the spectral condition number. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [4,5], and it is the kernel routine in the SIGMA library [6]. Some tuning parameters (marked with [TP]) are available for the implementer.

The computational range for the nonzero singular values is the machine number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even denormalized singular values can be computed with the corresponding gradual loss of accurate digits. Iconic, honestly.

🍓 Contributor

============

Zlatko Drmac (Zagreb, Croatia)

💜 References

[1] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer. SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
[2] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
[3] Z. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic. SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
[4] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. LAPACK Working note 169.
[5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. LAPACK Working note 170.
[6] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations. Department of Mathematics, University of Zagreb, 2008, 2015.

🎀 Bugs, examples and comments

===========================

Please report all bugs and send interesting test examples and comments to drmac@math.hr. Thank you. 💖

Definition at line 353 of file cgesvj.f.

   subroutine dgesvj (character*1 joba, character*1 jobu, character*1 jobv,
       integer m, integer n, double precision, dimension( lda, * ) a, integer
       lda, double precision, dimension( n ) sva, integer mv, double precision,
       dimension( ldv, * ) v, integer ldv, double precision, dimension( lwork )
       work, integer lwork, integer info)
       DGESVJ

✨ Purpose

DGESVJ 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^t,  [++] = [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. DGESVJ can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines, see below under Further Details.

🫶 Parameters

JOBA — CHARACTER*1. Specifies the structure of A.
- = 'L': The input matrix A is lower triangular;
- = 'U': The input matrix A is upper triangular;
- = 'G': The input matrix A is general M-by-N matrix, M >= N.
JOBU — CHARACTER*1. Specifies whether to compute the left singular vectors (columns of U):
- = 'U': The left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A. See more details in the description of A. The default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E').
- = 'C': Analogous to JOBU='U', except that user can control the level of numerical orthogonality of the computed left singular vectors. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK. No CTOL smaller than ONE is allowed. CTOL greater than 1 / EPS is meaningless. The option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations. See the descriptions of A and WORK(1).
- = 'N': The matrix U is not computed. However, see the description of A.
JOBV — CHARACTER*1. Specifies whether to compute the right singular vectors, that is, the matrix V:
- = 'V': the matrix V is computed and returned in the array V
- = 'A': the Jacobi rotations are applied to the MV-by-N array V. In other words, the right singular vector matrix V is not computed explicitly, instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V.
- = 'N': the matrix V is not computed and the array V is not referenced
M — INTEGER. The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
N — INTEGER. The number of columns of the input matrix A. M >= N >= 0.
A — the matrix itself; the rest of its description picks up in the next part, because this manpage is a whole six-part saga and we are not skipping a single crumb.

A is a DOUBLE PRECISION array, dimension (LDA,N) — okay so on entry this is literally just your M-by-N matrix A, the main character of the whole operation 💅.

On exit, though? It depends, bestie, so pay attention:

If JOBU = 'U' .OR. JOBU = 'C':
- If INFO = 0: the RANKA orthonormal columns of U get handed back to you in the leading RANKA columns of the array A. Here RANKA <= N is the number of computed singular values of A that are above the underflow threshold DLAMCH('S'). The singular vectors that correspond to underflowed or zero singular values? Not computed, periodt. The value of RANKA comes back to you in the array WORK as RANKA=NINT(WORK(2)) — also go peek at the descriptions of SVA and WORK. The computed columns of U are mutually numerically orthogonal up to approximately TOL=DSQRT(M)*EPS (that's the default) or TOL=CTOL*EPS (when JOBU = 'C'), see the description of JOBU ✨.
- If INFO > 0: the procedure DGESVJ did not converge in the given number of iterations (sweeps) 😭. In that case, the computed columns of U may not be orthogonal up to TOL. But like, the output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small. So not a total loss 🫶.
If JOBU = 'N':
- If INFO = 0: note that the left singular vectors are kinda 'for free' in the one-sided Jacobi SVD algorithm. However, if you only need the singular values, then the level of numerical orthogonality of U is not an issue at all, and the iterations stop when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS. So on exit, A contains the columns of U scaled with the corresponding singular values 🌸.
- If INFO > 0: the procedure DGESVJ did not converge in the given number of iterations (sweeps). Yeah. That one 💀.

LDA is INTEGER — the leading dimension of the array A. Rule of the day: LDA >= max(1,M).

SVA is a DOUBLE PRECISION array, dimension (N).

On exit:

If INFO = 0: depending on the value SCALE = WORK(1), we have:
- If SCALE = ONE: SVA(1:N) contains the computed singular values of A. And just so you know, during the computation SVA contains the Euclidean column norms of the iterated matrices in the array A 👀.
- If SCALE .NE. ONE: the singular values of A are SCALE*SVA(1:N), and this factored representation exists because some of the singular values of A might underflow or overflow. It's giving "we outsourced the scaling" 🎀.
If INFO > 0: the procedure DGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate 🙈.

MV is INTEGER — if JOBV = 'A', then the product of Jacobi rotations in DGESVJ is applied to the first MV rows of V. See the description of JOBV.

V is a DOUBLE PRECISION array, dimension (LDV,N).

If JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors 🦋.
If JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V.
If JOBV = 'N', then V is simply not referenced. Ghosted. Not even a vibe check 💀.

LDV is INTEGER — the leading dimension of the array V, so LDV >= 1.

If JOBV = 'V', then LDV >= max(1,N).
If JOBV = 'A', then LDV >= max(1,MV).

WORK is a DOUBLE PRECISION array, dimension (MAX(1,LWORK)) — this is the receipts array, bestie, it holds so much tea 🍓.

On entry:

If JOBU = 'C': WORK(1) = CTOL, where CTOL defines the threshold for convergence. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, with EPS=DLAMCH('E'). It is required that CTOL >= ONE, i.e. it is not allowed to force the routine to obtain orthogonality below EPS — nope, we don't do that here, the routine has boundaries 🚫.

On exit:

WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A. (See the description of SVA().)
WORK(2) = NINT(WORK(2)) is the number of the computed nonzero singular values.
WORK(3) = NINT(WORK(3)) is the number of the computed singular values that are larger than the underflow threshold.
WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence.
WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. This is useful info in cases when DGESVJ did not converge, because it can be used to estimate whether the output is still useful and for post festum analysis 👑.
WORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep. Also useful for a post festum analysis.

LWORK is INTEGER — the length of the array WORK.

LWORK >= 1, if MIN(M,N) = 0, and LWORK >= MAX(6,M+N), otherwise.

If on entry LWORK = -1, then a workspace query is assumed and no computation is done; WORK(1) is set to the minial (and optimal) length of WORK. (Yes, "minial" — that is fr how it's spelled in the source, bestie, I don't make the rules 😭)

INFO is INTEGER — the little status report at the end:

= 0: successful exit 🎉.
< 0: if INFO = -i, then the i-th argument had an illegal value 🙈.
> 0: DGESVJ did not converge in the maximal allowed number (30) of sweeps. The output may still be useful. See the description of WORK.

💖 Author 💖

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

✨ Further Details ✨

The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations. The rotations are implemented as fast scaled rotations of Anda and Park [1]. In the case of underflow of the Jacobi angle, a modified Jacobi transformation of Drmac [4] is used. The pivot strategy uses column interchanges of de Rijk [2]. The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [3] 🔥.

The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the spectral condition number.

The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. Some tuning parameters (marked with [TP]) are available for the implementer 🧚.

🌈 Contributors 🌈

============

Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) 🫶

🎀 References 🎀

[1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
    SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
[2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
    singular value decomposition on a vector computer.
    SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
[3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
[4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
    value computation in floating point arithmetic.
    SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
[5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
    SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
    LAPACK Working note 169.
[6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
    SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
    LAPACK Working note 170.
[7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
    QSVD, (H,K)-SVD computations.
    Department of Mathematics, University of Zagreb, 2008.

🦋 Bugs, examples and comments 🦋

===========================

Please report all bugs and send interesting test examples and comments to drmac@math.hr. Thank you 💌.

Definition at line 340 of file dgesvj.f.

subroutine sgesvj (character*1 joba, character*1 jobu, character*1 jobv,
    integer m, integer n, real, dimension( lda, * ) a, integer lda, real,
    dimension( n ) sva, integer mv, real, dimension( ldv, * ) v, integer
    ldv, real, dimension( lwork ) work, integer lwork, integer info)

⭐ `SGESVJ` ⭐

Purpose:

✨ Parameters ✨

this is the full parameter squad — every single one, no skipping, periodt 💅✨

🎀 JOBA 🎀

JOBA is CHARACTER*1 and it specifies the structure of A:

'L': the input matrix A is lower triangular;
'U': the input matrix A is upper triangular;
'G': the input matrix A is a general M-by-N matrix, M >= N.

💖 JOBU 💖

JOBU is CHARACTER*1 and it specifies whether to compute the left singular vectors — aka the columns of U:

'U': the left singular vectors corresponding to the nonzero singular values are computed and returned in the leading columns of A. see more details in the description of A. the default numerical orthogonality threshold is set to approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E').
'C': analogous to JOBU='U', except that the user can control the level of numerical orthogonality of the computed left singular vectors. TOL can be set to TOL = CTOL*EPS, where CTOL is given on input in the array WORK. no CTOL smaller than ONE is allowed. CTOL greater than 1 / EPS is meaningless. the option 'C' can be used if M*EPS is satisfactory orthogonality of the computed left singular vectors, so CTOL=M could save few sweeps of Jacobi rotations. see the descriptions of A and WORK(1).
'N': the matrix U is not computed. however, see the description of A.

🦋 JOBV 🦋

JOBV is CHARACTER*1 and it specifies whether to compute the right singular vectors, that is, the matrix V:

'V': the matrix V is computed and returned in the array V;
'A': the Jacobi rotations are applied to the MV-by-N array V. in other words, the right singular vector matrix V is not computed explicitly; instead it is applied to an MV-by-N matrix initially stored in the first MV rows of V.
'N': the matrix V is not computed and the array V is not referenced.

🌸 M 🌸

M is INTEGER. the number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.

🌈 N 🌈

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

👑 A 👑

A is a REAL array, dimension (LDA,N).

on entry, the M-by-N matrix A.

on exit:

if JOBU = 'U' .OR. JOBU = 'C':

if INFO = 0: RANKA orthonormal columns of U are returned in the leading RANKA columns of the array A. here RANKA <= N is the number of computed singular values of A that are above the underflow threshold SLAMCH('S'). the singular vectors corresponding to underflowed or zero singular values are not computed. the value of RANKA is returned in the array WORK as RANKA=NINT(WORK(2)). also see the descriptions of SVA and WORK. the computed columns of U are mutually numerically orthogonal up to approximately TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), see the description of JOBU.
if INFO > 0: the procedure SGESVJ did not converge in the given number of iterations (sweeps). in that case, the computed columns of U may not be orthogonal up to TOL. the output U (stored in A), SIGMA (given by the computed singular values in SVA(1:N)) and V is still a decomposition of the input matrix A in the sense that the residual ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.

if JOBU = 'N':

if INFO = 0: note that the left singular vectors are 'for free' in the one-sided Jacobi SVD algorithm. however, if only the singular values are needed, the level of numerical orthogonality of U is not an issue and iterations are stopped when the columns of the iterated matrix are numerically orthogonal up to approximately M*EPS. thus, on exit, A contains the columns of U scaled with the corresponding singular values.
if INFO > 0: the procedure SGESVJ did not converge in the given number of iterations (sweeps).

🍓 LDA 🍓

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

🧚 SVA 🧚

SVA is a REAL array, dimension (N).

on exit:

if INFO = 0: depending on the value SCALE = WORK(1), we have:

if SCALE = ONE: SVA(1:N) contains the computed singular values of A. during the computation SVA contains the Euclidean column norms of the iterated matrices in the array A.
if SCALE .NE. ONE: the singular values of A are SCALE*SVA(1:N), and this factored representation is due to the fact that some of the singular values of A might underflow or overflow.

if INFO > 0: the procedure SGESVJ did not converge in the given number of iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.

💜 MV 💜

MV is INTEGER. if JOBV = 'A', then the product of Jacobi rotations in SGESVJ is applied to the first MV rows of V. see the description of JOBV.

⭐ V ⭐

V is a REAL array, dimension (LDV,N).

if JOBV = 'V', then V contains on exit the N-by-N matrix of the right singular vectors;
if JOBV = 'A', then V contains the product of the computed right singular vector matrix and the initial matrix in the array V;
if JOBV = 'N', then V is not referenced.

🎀 LDV 🎀

LDV is INTEGER. the leading dimension of the array V, LDV >= 1. if JOBV = 'V', then LDV >= max(1,N). if JOBV = 'A', then LDV >= max(1,MV).

🔥 WORK 🔥

WORK is a REAL array, dimension (MAX(1,LWORK)).

on entry:

if JOBU = 'C': WORK(1) = CTOL, where CTOL defines the threshold for convergence. the process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=SLAMCH('E'). it is required that CTOL >= ONE, i.e. it is not allowed to force the routine to obtain orthogonality below EPSILON.

on exit:

WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A. (see description of SVA().)
WORK(2) = NINT(WORK(2)) is the number of the computed nonzero singular values.
WORK(3) = NINT(WORK(3)) is the number of the computed singular values that are larger than the underflow threshold.
WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence.
WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. this is useful information in cases when SGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis.
WORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep. it can be useful for a post festum analysis.

🫶 LWORK 🫶

LWORK is INTEGER. the length of WORK. LWORK >= 1, if MIN(M,N) = 0, and LWORK >= MAX(6,M+N), otherwise.

if on entry LWORK = -1, then a workspace query is assumed and no computation is done; WORK(1) is set to the minimal (and optimal) length of WORK.

😭 INFO 😭

INFO is INTEGER:

= 0: successful exit.
< 0: if INFO = -i, then the i-th argument had an illegal value.
> 0: SGESVJ did not converge in the maximal allowed number (30) of sweeps. the output may still be useful. see the description of WORK.

🌸 Author 🌸

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Further Details 🔬

ok so the tea on the actual math: the orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations. Those rotations are implemented as fast scaled rotations of Anda and Park [1]. And if the Jacobi angle underflows? Then a modified Jacobi transformation of Drmac [4] swoops in to save the day. The pivot strategy uses column interchanges of de Rijk [2]. 💫

The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [3]. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the spectral condition number. Bestie, say it with me: it's giving numerically stable. 👑

The computational range for the nonzero singular values is the machine number interval (UNDERFLOW, OVERFLOW). In extreme cases, even denormalized singular values can be computed with the corresponding gradual loss of accurate digits. 😭 (like yes she'll do it but the digits will be a little crunchy)

Contributors 🫶

Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

References 📚

[1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
[2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer. SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
[3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
[4] Z. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic. SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
[5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. LAPACK Working note 169.
[6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. LAPACK Working note 170.
[7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations. Department of Mathematics, University of Zagreb, 2008.

Bugs, Examples and Comments 🐛

Please report all bugs and send interesting test examples and comments to drmac@math.hr. Thank you. 💖

Definition at line 326 of file sgesvj.f.

subroutine zgesvj 💜

okok now for the complex cousin, and she is a whole moment. the signature, verbatim, no cap:

subroutine zgesvj (character*1 joba, character*1 jobu, character*1 jobv,
    integer m, integer n, complex*16, dimension( lda, * ) a, integer lda,
    double precision, dimension( n ) sva, integer mv, complex*16, dimension(
    ldv, * ) v, integer ldv, complex*16, dimension( lwork ) cwork, integer
    lwork, double precision, dimension( lrwork ) rwork, integer lrwork,
    integer info)

ZGESVJ

Purpose ✨

!>
!> ZGESVJ 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.
!>

so like, SIGMA is diagonal, U is M-by-N orthonormal, V is N-by-N unitary, and the diagonal girlies of SIGMA are literally the singular values. the columns of U and V are the left and right singular vectors respectively. iykyk 💅

Parameters 🌸

JOBA

!>          JOBA is CHARACTER*1
!>          Specifies the structure of A.
!>          = 'L': The input matrix A is lower triangular;
!>          = 'U': The input matrix A is upper triangular;
!>          = 'G': The input matrix A is general M-by-N matrix, M >= N.
!>

aka: what shape is A giving today?

JOBU

!>          JOBU is CHARACTER*1
!>          Specifies whether to compute the left singular vectors
!>          (columns of U):
!>          = 'U' or 'F': The left singular vectors corresponding to the nonzero
!>                 singular values are computed and returned in the leading
!>                 columns of A. See more details in the description of A.
!>                 The default numerical orthogonality threshold is set to
!>                 approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=DLAMCH('E').
!>          = 'C': Analogous to JOBU='U', except that user can control the
!>                 level of numerical orthogonality of the computed left
!>                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
!>                 CTOL is given on input in the array WORK.
!>                 No CTOL smaller than ONE is allowed. CTOL greater
!>                 than 1 / EPS is meaningless. The option 'C'
!>                 can be used if M*EPS is satisfactory orthogonality
!>                 of the computed left singular vectors, so CTOL=M could
!>                 save few sweeps of Jacobi rotations.
!>                 See the descriptions of A and WORK(1).
!>          = 'N': The matrix U is not computed. However, see the
!>                 description of A.
!>

careful with 'C': no CTOL smaller than ONE is allowed, and CTOL greater than 1 / EPS is meaningless. them's the rules 🫣

JOBV

!>          JOBV is CHARACTER*1
!>          Specifies whether to compute the right singular vectors, that
!>          is, the matrix V:
!>          = 'V' or 'J': the matrix V is computed and returned in the array V
!>          = 'A':  the Jacobi rotations are applied to the MV-by-N
!>                  array V. In other words, the right singular vector
!>                  matrix V is not computed explicitly; instead it is
!>                  applied to an MV-by-N matrix initially stored in the
!>                  first MV rows of V.
!>          = 'N':  the matrix V is not computed and the array V is not
!>                  referenced
!>

M

!>          M is INTEGER
!>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
!>

N

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

A

!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit,
!>          If JOBU = 'U' .OR. JOBU = 'C':
!>                 If INFO = 0 :
!>                 RANKA orthonormal columns of U are returned in the
!>                 leading RANKA columns of the array A. Here RANKA <= N
!>                 is the number of computed singular values of A that are
!>                 above the underflow threshold DLAMCH('S'). The singular
!>                 vectors corresponding to underflowed or zero singular
!>                 values are not computed. The value of RANKA is returned
!>                 in the array RWORK as RANKA=NINT(RWORK(2)). Also see the
!>                 descriptions of SVA and RWORK. The computed columns of U
!>                 are mutually numerically orthogonal up to approximately
!>                 TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
!>                 see the description of JOBU.
!>                 If INFO > 0,
!>                 the procedure ZGESVJ did not converge in the given number
!>                 of iterations (sweeps). In that case, the computed
!>                 columns of U may not be orthogonal up to TOL. The output
!>                 U (stored in A), SIGMA (given by the computed singular
!>                 values in SVA(1:N)) and V is still a decomposition of the
!>                 input matrix A in the sense that the residual
!>                 || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small.
!>          If JOBU = 'N':
!>                 If INFO = 0 :
!>                 Note that the left singular vectors are 'for free' in the
!>                 one-sided Jacobi SVD algorithm. However, if only the
!>                 singular values are needed, the level of numerical
!>                 orthogonality of U is not an issue and iterations are
!>                 stopped when the columns of the iterated matrix are
!>                 numerically orthogonal up to approximately M*EPS. Thus,
!>                 on exit, A contains the columns of U scaled with the
!>                 corresponding singular values.
!>                 If INFO > 0:
!>                 the procedure ZGESVJ did not converge in the given number
!>                 of iterations (sweeps).
!>

highkey love that for JOBU = 'N' — the left singular vectors come for free, bestie. no cap, free real estate. 🎁

LDA

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

SVA

!>          SVA is DOUBLE PRECISION array, dimension (N)
!>          On exit,
!>          If INFO = 0 :
!>          depending on the value SCALE = RWORK(1), we have:
!>                 If SCALE = ONE:
!>                 SVA(1:N) contains the computed singular values of A.
!>                 During the computation SVA contains the Euclidean column
!>                 norms of the iterated matrices in the array A.
!>                 If SCALE .NE. ONE:
!>                 The singular values of A are SCALE*SVA(1:N), and this
!>                 factored representation is due to the fact that some of the
!>                 singular values of A might underflow or overflow.
!>
!>          If INFO > 0:
!>          the procedure ZGESVJ did not converge in the given number of
!>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
!>

reminder that if SCALE .NE. ONE the actual singular values are SCALE*SVA(1:N) — do not forget the SCALE or you will be telling lies on the timeline 🙈

MV

!>          MV is INTEGER
!>          If JOBV = 'A', then the product of Jacobi rotations in ZGESVJ
!>          is applied to the first MV rows of V. See the description of JOBV.
!>

V

!>          V is COMPLEX*16 array, dimension (LDV,N)
!>          If JOBV = 'V', then V contains on exit the N-by-N matrix of
!>                         the right singular vectors;
!>          If JOBV = 'A', then V contains the product of the computed right
!>                         singular vector matrix and the initial matrix in
!>                         the array V.
!>          If JOBV = 'N', then V is not referenced.
!>

LDV

LDV — is INTEGER. The leading dimension of the array V, LDV >= 1. If JOBV = 'V', then LDV >= MAX(1,N). If JOBV = 'A', then LDV >= MAX(1,MV).
CWORK — is COMPLEX*16 array, dimension (MAX(1,LWORK)). Used as workspace. Slay, it holds the complex workspace tea.
LWORK — is INTEGER. Length of CWORK. LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M+N, otherwise.

If on entry LWORK = -1, then a workspace query is assumed and no computation is done; CWORK(1) is set to the minial (and optimal) length of CWORK. I know, “minial,” but that’s literally what the page says, bestie — no editing the lore.

RWORK — is DOUBLE PRECISION array, dimension (max(6,LRWORK)).

On entry, if JOBU = 'C': RWORK(1) = CTOL, where CTOL defines the threshold for convergence. The process stops if all columns of A are mutually orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). It is required that CTOL >= ONE, i.e. it is not allowed to force the routine to obtain orthogonality below EPSILON. Respect the threshold, okay?

On exit:

RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A. (See description of SVA().)
RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero singular values.
RWORK(3) = NINT(RWORK(3)) is the number of the computed singular values that are larger than the underflow threshold.
RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi rotations needed for numerical convergence.
RWORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. This is useful information in cases when ZGESVJ did not converge, as it can be used to estimate whether the output is still useful and for post festum analysis.
RWORK(6) = the largest absolute value over all sines of the Jacobi rotation angles in the last sweep. It can be useful for a post festum analysis.

LRWORK — is INTEGER. Length of RWORK. LRWORK >= 1, if MIN(M,N) = 0, and LRWORK >= MAX(6,N), otherwise.

If on entry LRWORK = -1, then a workspace query is assumed and no computation is done; RWORK(1) is set to the minial (and optimal) length of RWORK. Yes, again with “minial” — it’s giving source-code lore, not autocorrect.

INFO — is INTEGER
- = 0: successful exit. We love a clean win.
- < 0: if INFO = -i, then the i-th argument had an illegal value.
- > 0: ZGESVJ did not converge in the maximal allowed number (NSWEEP=30) of sweeps. The output may still be useful. See the description of RWORK. Not a total flop — but definitely a “read the notes” moment.

👑 Author 👑

Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

✨ Further Details ✨

The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane rotations. In the case of underflow of the tangent of the Jacobi angle, a modified Jacobi transformation of Drmac [3] is used. Pivot strategy uses column interchanges of de Rijk [1]. The relative accuracy of the computed singular values and the accuracy of the computed singular vectors (in angle metric) is as guaranteed by the theory of Demmel and Veselic [2]. The condition number that determines the accuracy in the full rank case is essentially min_{D=diag} kappa(A*D), where kappa(.) is the spectral condition number. The best performance of this Jacobi SVD procedure is achieved if used in an accelerated version of Drmac and Veselic [4,5], and it is the kernel routine in the SIGMA library [6]. Some tuning parameters (marked with [TP]) are available for the implementer. The computational range for the nonzero singular values is the machine number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even denormalized singular values can be computed with the corresponding gradual loss of accurate digits. Numerical accuracy lore, but make it fashion. ✨

💖 Contributor 💖

Zlatko Drmac (Zagreb, Croatia)

📚 References 📚

P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the singular value decomposition on a vector computer. SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
Z. Drmac: Implementation of Jacobi rotations for accurate singular value computation in floating point arithmetic. SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. LAPACK Working note 169.
Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. LAPACK Working note 170.
Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, QSVD, (H,K)-SVD computations. Department of Mathematics, University of Zagreb, 2008, 2015.

🐛 Bugs, examples and comments 🐛

===========================

Please report all bugs and send interesting test examples and comments to drmac@math.hr. Thank you. Bug reports are literally a love language here.

Definition at line 353 of file zgesvj.f.

🫶 Author 🫶

Generated automatically by Doxygen for LAPACK from the source code.

💖 cute sign-off time: the source/version energy is LAPACK, Version 3.12.0, page gesvj(3), and the footer date field stays exactly as this snippet gives it — not shown here, so I’m not inventing one, bestie. Periodt. 🫶

📦 where this page comes from

arch rolling: lapack-doc 3.12.1-2 · custom
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash-1017 on 2026-11-04: a language model's rewrite of the original, under its terms

gesvj(3)

🌸 NAME 🌸

✨ SYNOPSIS ✨

🦋 Functions 🦋

💫 DETAILED DESCRIPTION 💫

👑 Function Documentation 👑

🌟 CGESVJ 🌟

🌸 JOBA 🌸

🦋 JOBU 🦋

🎀 JOBV 🎀

🐣 M 🐣

🐣 N 🐣

💜 A 💜

🧚 LDA 🧚

🍓 SVA 🍓

🌈 MV 🌈

🦄 V 🦄

🌟 LDV 🌟

💠 CWORK 💠

🔥 LWORK 🔥

🎀 RWORK 🎀

👑 Author

🧚 Further Details

🍓 Contributor

💜 References

🎀 Bugs, examples and comments

✨ Purpose

🫶 Parameters

💖 Author 💖

✨ Further Details ✨

🌈 Contributors 🌈

🎀 References 🎀

🦋 Bugs, examples and comments 🦋

⭐ SGESVJ ⭐

Purpose:

✨ Parameters ✨

🎀 JOBA 🎀

💖 JOBU 💖

🦋 JOBV 🦋

🌸 M 🌸

🌈 N 🌈

👑 A 👑

🍓 LDA 🍓

🧚 SVA 🧚

💜 MV 💜

⭐ V ⭐

🎀 LDV 🎀

🔥 WORK 🔥

🫶 LWORK 🫶

😭 INFO 😭

🌸 Author 🌸

Further Details 🔬

Contributors 🫶

References 📚

Bugs, Examples and Comments 🐛

subroutine zgesvj 💜

Purpose ✨

Parameters 🌸

JOBA

JOBU

JOBV

M

N

A

LDA

SVA

MV

V

LDV

👑 Author 👑

✨ Further Details ✨

💖 Contributor 💖

📚 References 📚

🐛 Bugs, examples and comments 🐛

🫶 Author 🫶

🌟 `CGESVJ` 🌟

🌸 `JOBA` 🌸

🦋 `JOBU` 🦋

🎀 `JOBV` 🎀

🐣 `M` 🐣

🐣 `N` 🐣

💜 `A` 💜

🧚 `LDA` 🧚

🍓 `SVA` 🍓

🌈 `MV` 🌈

🦄 `V` 🦄

🌟 `LDV` 🌟

💠 `CWORK` 💠

🔥 `LWORK` 🔥

🎀 `RWORK` 🎀

⭐ `SGESVJ` ⭐