gesvj#

sgesvj

Functions

void sgesvj(const char *joba, const char *jobu, const char *jobv, const INT m, const INT n, f32 *restrict A, const INT lda, f32 *restrict SVA, const INT mv, f32 *restrict V, const INT ldv, f32 *restrict work, const INT lwork, INT *info)#

dgesvj

Functions

void dgesvj(const char *joba, const char *jobu, const char *jobv, const INT m, const INT n, f64 *restrict A, const INT lda, f64 *restrict SVA, const INT mv, f64 *restrict V, const INT ldv, f64 *restrict work, const INT lwork, INT *info)#

cgesvj

Functions

void cgesvj(
    const char*          joba,
    const char*          jobu,
    const char*          jobv,
    const INT            m,
    const INT            n,
          c64*  restrict A,
    const INT            lda,
          f32*  restrict SVA,
    const INT            mv,
          c64*  restrict V,
    const INT            ldv,
          c64*  restrict cwork,
    const INT            lwork,
          f32*  restrict rwork,
    const INT            lrwork,
          INT*           info
);

void cgesvj(const char *joba, const char *jobu, const char *jobv, const INT m, const INT n, c64 *restrict A, const INT lda, f32 *restrict SVA, const INT mv, c64 *restrict V, const INT ldv, c64 *restrict cwork, const INT lwork, f32 *restrict rwork, const INT lrwork, INT *info)#

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

         A = U * SIGMA * V^*,

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. CGESVJ can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines, see below under Further Details.

Parameters

in

joba

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.

in

jobu

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 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 RWORK. 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 RWORK(1). = ‘N’: The matrix U is not computed. However, see the description of A.

in

jobv

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

in

m

The number of rows of the input matrix A. 1/SLAMCH(‘E’) > M >= 0.

in

n

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

inout

A

Complex*16 array, dimension (lda, n). On entry, the M-by-N matrix A. On exit: If JOBU = ‘U’ or JOBU = ‘F’ 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).

in

lda

The leading dimension of the array A. lda >= max(1, m).

out

SVA

Single 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 /= 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.

in

mv

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.

inout

V

Complex*16 array, dimension (ldv, n). If JOBV = ‘V’ or ‘J’, 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.

in

ldv

The leading dimension of the array V, ldv >= 1. If JOBV = ‘V’ or ‘J’, then ldv >= max(1, n). If JOBV = ‘A’, then ldv >= max(1, mv).

inout

cwork

Complex*16 array, dimension (max(1, lwork)). Used as workspace.

in

lwork

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 minimal (and optimal) length of CWORK.

inout

rwork

Single precision array, dimension (max(6, lrwork)). On entry: If JOBU = ‘C’: rwork[0] = 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: rwork[0] = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A. (See description of SVA.) rwork[1] = NINT(rwork[1]) is the number of the computed nonzero singular values. rwork[2] = NINT(rwork[2]) is the number of the computed singular values that are larger than the underflow threshold. rwork[3] = NINT(rwork[3]) is the number of sweeps of Jacobi rotations needed for numerical convergence. rwork[4] = max_{i/=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[5] = 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.

in

lrwork

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 minimal (and optimal) length of RWORK.

out

info

= 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 (30) of sweeps. The output may still be useful. See the description of rwork.

zgesvj

Functions

void zgesvj(
    const char*          joba,
    const char*          jobu,
    const char*          jobv,
    const INT            m,
    const INT            n,
          c128* restrict A,
    const INT            lda,
          f64*  restrict SVA,
    const INT            mv,
          c128* restrict V,
    const INT            ldv,
          c128* restrict cwork,
    const INT            lwork,
          f64*  restrict rwork,
    const INT            lrwork,
          INT*           info
);

void zgesvj(const char *joba, const char *jobu, const char *jobv, const INT m, const INT n, c128 *restrict A, const INT lda, f64 *restrict SVA, const INT mv, c128 *restrict V, const INT ldv, c128 *restrict cwork, const INT lwork, f64 *restrict rwork, const INT lrwork, INT *info)#

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

         A = U * SIGMA * V^*,

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. ZGESVJ can sometimes compute tiny singular values and their singular vectors much more accurately than other SVD routines, see below under Further Details.

Parameters

in

joba

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.

in

jobu

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 RWORK. 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 RWORK(1). = ‘N’: The matrix U is not computed. However, see the description of A.

in

jobv

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

in

m

The number of rows of the input matrix A. 1/DLAMCH(‘E’) > M >= 0.

in

n

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

inout

A

Complex*16 array, dimension (lda, n). On entry, the M-by-N matrix A. On exit: If JOBU = ‘U’ or JOBU = ‘F’ 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).

in

lda

The leading dimension of the array A. lda >= max(1, m).

out

SVA

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

in

mv

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.

inout

V

Complex*16 array, dimension (ldv, n). If JOBV = ‘V’ or ‘J’, 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.

in

ldv

The leading dimension of the array V, ldv >= 1. If JOBV = ‘V’ or ‘J’, then ldv >= max(1, n). If JOBV = ‘A’, then ldv >= max(1, mv).

inout

cwork

Complex*16 array, dimension (max(1, lwork)). Used as workspace.

in

lwork

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 minimal (and optimal) length of CWORK.

inout

rwork

Double precision array, dimension (max(6, lrwork)). On entry: If JOBU = ‘C’: rwork[0] = 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. On exit: rwork[0] = SCALE is the scaling factor such that SCALE*SVA(1:N) are the computed singular values of A. (See description of SVA.) rwork[1] = NINT(rwork[1]) is the number of the computed nonzero singular values. rwork[2] = NINT(rwork[2]) is the number of the computed singular values that are larger than the underflow threshold. rwork[3] = NINT(rwork[3]) is the number of sweeps of Jacobi rotations needed for numerical convergence. rwork[4] = max_{i/=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[5] = 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.

in

lrwork

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 minimal (and optimal) length of RWORK.

out

info

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

gesvj#

This Page