- Generated by
1.15.0
Functions | |
| subroutine | zgejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, cwork, lwork, rwork, lrwork, iwork, info) |
| ZGEJSV | |
| subroutine | zgesdd (jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info) |
| ZGESDD | |
| subroutine | zgesvd (jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info) |
| ZGESVD computes the singular value decomposition (SVD) 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 | |
| subroutine | zgesvdx (jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info) |
| ZGESVDX computes the singular value decomposition (SVD) for GE matrices | |
| subroutine | zggsvd3 (jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork, rwork, iwork, info) |
| ZGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices | |
This is the group of complex16 singular value driver functions for GE matrices
| subroutine zgejsv | ( | character*1 | joba, |
| character*1 | jobu, | ||
| character*1 | jobv, | ||
| character*1 | jobr, | ||
| character*1 | jobt, | ||
| character*1 | jobp, | ||
| integer | m, | ||
| integer | n, | ||
| complex*16, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| double precision, dimension( n ) | sva, | ||
| complex*16, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| complex*16, dimension( ldv, * ) | v, | ||
| integer | ldv, | ||
| complex*16, dimension( lwork ) | cwork, | ||
| integer | lwork, | ||
| double precision, dimension( lrwork ) | rwork, | ||
| integer | lrwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
ZGEJSV
Download ZGEJSV + dependencies [TGZ] [ZIP] [TXT]
!> !> ZGEJSV 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 (M-by-N) matrix which is zero except for its N !> diagonal elements, [U] is an M-by-N (or M-by-M) unitary 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. The matrices [U] and [V] !> are computed and stored in the arrays U and V, respectively. The diagonal !> of [SIGMA] is computed and stored in the array SVA. !>
| [in] | JOBA | !> JOBA is CHARACTER*1
!> Specifies the level of accuracy:
!> = 'C': This option works well (high relative accuracy) if A = B * D,
!> with well-conditioned B and arbitrary diagonal matrix D.
!> The accuracy cannot be spoiled by COLUMN scaling. The
!> accuracy of the computed output depends on the condition of
!> B, and the procedure aims at the best theoretical accuracy.
!> The relative error max_{i=1:N}|d sigma_i| / sigma_i is
!> bounded by f(M,N)*epsilon* cond(B), independent of D.
!> The input matrix is preprocessed with the QRF with column
!> pivoting. This initial preprocessing and preconditioning by
!> a rank revealing QR factorization is common for all values of
!> JOBA. Additional actions are specified as follows:
!> = 'E': Computation as with 'C' with an additional estimate of the
!> condition number of B. It provides a realistic error bound.
!> = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
!> D1, D2, and well-conditioned matrix C, this option gives
!> higher accuracy than the 'C' option. If the structure of the
!> input matrix is not known, and relative accuracy is
!> desirable, then this option is advisable. The input matrix A
!> is preprocessed with QR factorization with FULL (row and
!> column) pivoting.
!> = 'G': Computation as with 'F' with an additional estimate of the
!> condition number of B, where A=B*D. If A has heavily weighted
!> rows, then using this condition number gives too pessimistic
!> error bound.
!> = 'A': Small singular values are not well determined by the data
!> and are considered as noisy; the matrix is treated as
!> numerically rank deficient. The error in the computed
!> singular values is bounded by f(m,n)*epsilon*||A||.
!> The computed SVD A = U * S * V^* restores A up to
!> f(m,n)*epsilon*||A||.
!> This gives the procedure the licence to discard (set to zero)
!> all singular values below N*epsilon*||A||.
!> = 'R': Similar as in 'A'. Rank revealing property of the initial
!> QR factorization is used do reveal (using triangular factor)
!> a gap sigma_{r+1} < epsilon * sigma_r in which case the
!> numerical RANK is declared to be r. The SVD is computed with
!> absolute error bounds, but more accurately than with 'A'.
!> |
| [in] | JOBU | !> JOBU is CHARACTER*1 !> Specifies whether to compute the columns of U: !> = 'U': N columns of U are returned in the array U. !> = 'F': full set of M left sing. vectors is returned in the array U. !> = 'W': U may be used as workspace of length M*N. See the description !> of U. !> = 'N': U is not computed. !> |
| [in] | JOBV | !> JOBV is CHARACTER*1 !> Specifies whether to compute the matrix V: !> = 'V': N columns of V are returned in the array V; Jacobi rotations !> are not explicitly accumulated. !> = 'J': N columns of V are returned in the array V, but they are !> computed as the product of Jacobi rotations, if JOBT = 'N'. !> = 'W': V may be used as workspace of length N*N. See the description !> of V. !> = 'N': V is not computed. !> |
| [in] | JOBR | !> JOBR is CHARACTER*1
!> Specifies the RANGE for the singular values. Issues the licence to
!> set to zero small positive singular values if they are outside
!> specified range. If A .NE. 0 is scaled so that the largest singular
!> value of c*A is around SQRT(BIG), BIG=DLAMCH('O'), then JOBR issues
!> the licence to kill columns of A whose norm in c*A is less than
!> SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
!> where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E').
!> = 'N': Do not kill small columns of c*A. This option assumes that
!> BLAS and QR factorizations and triangular solvers are
!> implemented to work in that range. If the condition of A
!> is greater than BIG, use ZGESVJ.
!> = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)]
!> (roughly, as described above). This option is recommended.
!> ===========================
!> For computing the singular values in the FULL range [SFMIN,BIG]
!> use ZGESVJ.
!> |
| [in] | JOBT | !> JOBT is CHARACTER*1 !> If the matrix is square then the procedure may determine to use !> transposed A if A^* seems to be better with respect to convergence. !> If the matrix is not square, JOBT is ignored. !> The decision is based on two values of entropy over the adjoint !> orbit of A^* * A. See the descriptions of WORK(6) and WORK(7). !> = 'T': transpose if entropy test indicates possibly faster !> convergence of Jacobi process if A^* is taken as input. If A is !> replaced with A^*, then the row pivoting is included automatically. !> = 'N': do not speculate. !> The option 'T' can be used to compute only the singular values, or !> the full SVD (U, SIGMA and V). For only one set of singular vectors !> (U or V), the caller should provide both U and V, as one of the !> matrices is used as workspace if the matrix A is transposed. !> The implementer can easily remove this constraint and make the !> code more complicated. See the descriptions of U and V. !> In general, this option is considered experimental, and 'N'; should !> be preferred. This is subject to changes in the future. !> |
| [in] | JOBP | !> JOBP is CHARACTER*1 !> Issues the licence to introduce structured perturbations to drown !> denormalized numbers. This licence should be active if the !> denormals are poorly implemented, causing slow computation, !> especially in cases of fast convergence (!). For details see [1,2]. !> For the sake of simplicity, this perturbations are included only !> when the full SVD or only the singular values are requested. The !> implementer/user can easily add the perturbation for the cases of !> computing one set of singular vectors. !> = 'P': introduce perturbation !> = 'N': do not perturb !> |
| [in] | M | !> M is INTEGER !> The number of rows of the input matrix A. M >= 0. !> |
| [in] | N | !> N is INTEGER !> The number of columns of the input matrix A. M >= N >= 0. !> |
| [in,out] | A | !> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> |
| [in] | LDA | !> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [out] | SVA | !> SVA is DOUBLE PRECISION array, dimension (N) !> On exit, !> - For WORK(1)/WORK(2) = ONE: The singular values of A. During the !> computation SVA contains Euclidean column norms of the !> iterated matrices in the array A. !> - For WORK(1) .NE. WORK(2): The singular values of A are !> (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if !> sigma_max(A) overflows or if small singular values have been !> saved from underflow by scaling the input matrix A. !> - If JOBR='R' then some of the singular values may be returned !> as exact zeros obtained by because they are !> below the numerical rank threshold or are denormalized numbers. !> |
| [out] | U | !> U is COMPLEX*16 array, dimension ( LDU, N ) !> If JOBU = 'U', then U contains on exit the M-by-N matrix of !> the left singular vectors. !> If JOBU = 'F', then U contains on exit the M-by-M matrix of !> the left singular vectors, including an ONB !> of the orthogonal complement of the Range(A). !> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N), !> then U is used as workspace if the procedure !> replaces A with A^*. In that case, [V] is computed !> in U as left singular vectors of A^* and then !> copied back to the V array. This 'W' option is just !> a reminder to the caller that in this case U is !> reserved as workspace of length N*N. !> If JOBU = 'N' U is not referenced, unless JOBT='T'. !> |
| [in] | LDU | !> LDU is INTEGER !> The leading dimension of the array U, LDU >= 1. !> IF JOBU = 'U' or 'F' or 'W', then LDU >= M. !> |
| [out] | V | !> V is COMPLEX*16 array, dimension ( LDV, N ) !> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of !> the right singular vectors; !> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N), !> then V is used as workspace if the pprocedure !> replaces A with A^*. In that case, [U] is computed !> in V as right singular vectors of A^* and then !> copied back to the U array. This 'W' option is just !> a reminder to the caller that in this case V is !> reserved as workspace of length N*N. !> If JOBV = 'N' V is not referenced, unless JOBT='T'. !> |
| [in] | LDV | !> LDV is INTEGER !> The leading dimension of the array V, LDV >= 1. !> If JOBV = 'V' or 'J' or 'W', then LDV >= N. !> |
| [out] | CWORK | !> CWORK is COMPLEX*16 array, dimension (MAX(2,LWORK)) !> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or !> LRWORK=-1), then on exit CWORK(1) contains the required length of !> CWORK for the job parameters used in the call. !> |
| [in] | LWORK | !> LWORK is INTEGER !> Length of CWORK to confirm proper allocation of workspace. !> LWORK depends on the job: !> !> 1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and !> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'): !> LWORK >= 2*N+1. This is the minimal requirement. !> ->> For optimal performance (blocked code) the optimal value !> is LWORK >= N + (N+1)*NB. Here NB is the optimal !> block size for ZGEQP3 and ZGEQRF. !> In general, optimal LWORK is computed as !> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ)). !> 1.2. .. an estimate of the scaled condition number of A is !> required (JOBA='E', or 'G'). In this case, LWORK the minimal !> requirement is LWORK >= N*N + 2*N. !> ->> For optimal performance (blocked code) the optimal value !> is LWORK >= max(N+(N+1)*NB, N*N+2*N)=N**2+2*N. !> In general, the optimal length LWORK is computed as !> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ), !> N*N+LWORK(ZPOCON)). !> 2. If SIGMA and the right singular vectors are needed (JOBV = 'V'), !> (JOBU = 'N') !> 2.1 .. no scaled condition estimate requested (JOBE = 'N'): !> -> the minimal requirement is LWORK >= 3*N. !> -> For optimal performance, !> LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, !> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ, !> ZUNMLQ. In general, the optimal length LWORK is computed as !> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZGESVJ), !> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)). !> 2.2 .. an estimate of the scaled condition number of A is !> required (JOBA='E', or 'G'). !> -> the minimal requirement is LWORK >= 3*N. !> -> For optimal performance, !> LWORK >= max(N+(N+1)*NB, 2*N,2*N+N*NB)=2*N+N*NB, !> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ, !> ZUNMLQ. In general, the optimal length LWORK is computed as !> LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ), !> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)). !> 3. If SIGMA and the left singular vectors are needed !> 3.1 .. no scaled condition estimate requested (JOBE = 'N'): !> -> the minimal requirement is LWORK >= 3*N. !> -> For optimal performance: !> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, !> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR. !> In general, the optimal length LWORK is computed as !> LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). !> 3.2 .. an estimate of the scaled condition number of A is !> required (JOBA='E', or 'G'). !> -> the minimal requirement is LWORK >= 3*N. !> -> For optimal performance: !> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, !> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR. !> In general, the optimal length LWORK is computed as !> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON), !> 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). !> 4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and !> 4.1. if JOBV = 'V' !> the minimal requirement is LWORK >= 5*N+2*N*N. !> 4.2. if JOBV = 'J' the minimal requirement is !> LWORK >= 4*N+N*N. !> In both cases, the allocated CWORK can accommodate blocked runs !> of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ. !> !> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or !> LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the !> minimal length of CWORK for the job parameters used in the call. !> |
| [out] | RWORK | !> RWORK is DOUBLE PRECISION array, dimension (MAX(7,LWORK)) !> On exit, !> RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1) !> such that SCALE*SVA(1:N) are the computed singular values !> of A. (See the description of SVA().) !> RWORK(2) = See the description of RWORK(1). !> RWORK(3) = SCONDA is an estimate for the condition number of !> column equilibrated A. (If JOBA = 'E' or 'G') !> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). !> It is computed using ZPOCON. It holds !> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA !> where R is the triangular factor from the QRF of A. !> However, if R is truncated and the numerical rank is !> determined to be strictly smaller than N, SCONDA is !> returned as -1, thus indicating that the smallest !> singular values might be lost. !> !> If full SVD is needed, the following two condition numbers are !> useful for the analysis of the algorithm. They are provided for !> a developer/implementer who is familiar with the details of !> the method. !> !> RWORK(4) = an estimate of the scaled condition number of the !> triangular factor in the first QR factorization. !> RWORK(5) = an estimate of the scaled condition number of the !> triangular factor in the second QR factorization. !> The following two parameters are computed if JOBT = 'T'. !> They are provided for a developer/implementer who is familiar !> with the details of the method. !> RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy !> of diag(A^* * A) / Trace(A^* * A) taken as point in the !> probability simplex. !> RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).) !> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or !> LRWORK=-1), then on exit RWORK(1) contains the required length of !> RWORK for the job parameters used in the call. !> |
| [in] | LRWORK | !> LRWORK is INTEGER !> Length of RWORK to confirm proper allocation of workspace. !> LRWORK depends on the job: !> !> 1. If only the singular values are requested i.e. if !> LSAME(JOBU,'N') .AND. LSAME(JOBV,'N') !> then: !> 1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), !> then: LRWORK = max( 7, 2 * M ). !> 1.2. Otherwise, LRWORK = max( 7, N ). !> 2. If singular values with the right singular vectors are requested !> i.e. if !> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND. !> .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) !> then: !> 2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), !> then LRWORK = max( 7, 2 * M ). !> 2.2. Otherwise, LRWORK = max( 7, N ). !> 3. If singular values with the left singular vectors are requested, i.e. if !> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. !> .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) !> then: !> 3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), !> then LRWORK = max( 7, 2 * M ). !> 3.2. Otherwise, LRWORK = max( 7, N ). !> 4. If singular values with both the left and the right singular vectors !> are requested, i.e. if !> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. !> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) !> then: !> 4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), !> then LRWORK = max( 7, 2 * M ). !> 4.2. Otherwise, LRWORK = max( 7, N ). !> !> If, on entry, LRWORK = -1 or LWORK=-1, a workspace query is assumed and !> the length of RWORK is returned in RWORK(1). !> |
| [out] | IWORK | !> IWORK is INTEGER array, of dimension at least 4, that further depends !> on the job: !> !> 1. If only the singular values are requested then: !> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) !> then the length of IWORK is N+M; otherwise the length of IWORK is N. !> 2. If the singular values and the right singular vectors are requested then: !> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) !> then the length of IWORK is N+M; otherwise the length of IWORK is N. !> 3. If the singular values and the left singular vectors are requested then: !> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) !> then the length of IWORK is N+M; otherwise the length of IWORK is N. !> 4. If the singular values with both the left and the right singular vectors !> are requested, then: !> 4.1. If LSAME(JOBV,'J') the length of IWORK is determined as follows: !> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) !> then the length of IWORK is N+M; otherwise the length of IWORK is N. !> 4.2. If LSAME(JOBV,'V') the length of IWORK is determined as follows: !> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) !> then the length of IWORK is 2*N+M; otherwise the length of IWORK is 2*N. !> !> On exit, !> IWORK(1) = the numerical rank determined after the initial !> QR factorization with pivoting. See the descriptions !> of JOBA and JOBR. !> IWORK(2) = the number of the computed nonzero singular values !> IWORK(3) = if nonzero, a warning message: !> If IWORK(3) = 1 then some of the column norms of A !> were denormalized floats. The requested high accuracy !> is not warranted by the data. !> IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to !> do the job as specified by the JOB parameters. !> If the call to ZGEJSV is a workspace query (indicated by LWORK = -1 or !> LRWORK = -1), then on exit IWORK(1) contains the required length of !> IWORK for the job parameters used in the call. !> |
| [out] | INFO | !> INFO is INTEGER !> < 0: if INFO = -i, then the i-th argument had an illegal value. !> = 0: successful exit; !> > 0: ZGEJSV did not converge in the maximal allowed number !> of sweeps. The computed values may be inaccurate. !> |
!>
!> ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3,
!> ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an
!> additional row pivoting can be used as a preprocessor, which in some
!> cases results in much higher accuracy. An example is matrix A with the
!> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
!> diagonal matrices and C is well-conditioned matrix. In that case, complete
!> pivoting in the first QR factorizations provides accuracy dependent on the
!> condition number of C, and independent of D1, D2. Such higher accuracy is
!> not completely understood theoretically, but it works well in practice.
!> Further, if A can be written as A = B*D, with well-conditioned B and some
!> diagonal D, then the high accuracy is guaranteed, both theoretically and
!> in software, independent of D. For more details see [1], [2].
!> The computational range for the singular values can be the full range
!> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
!> & LAPACK routines called by ZGEJSV are implemented to work in that range.
!> If that is not the case, then the restriction for safe computation with
!> the singular values in the range of normalized IEEE numbers is that the
!> spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
!> overflow. This code (ZGEJSV) is best used in this restricted range,
!> meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
!> returned as zeros. See JOBR for details on this.
!> Further, this implementation is somewhat slower than the one described
!> in [1,2] due to replacement of some non-LAPACK components, and because
!> the choice of some tuning parameters in the iterative part (ZGESVJ) is
!> left to the implementer on a particular machine.
!> The rank revealing QR factorization (in this code: ZGEQP3) should be
!> implemented as in [3]. We have a new version of ZGEQP3 under development
!> that is more robust than the current one in LAPACK, with a cleaner cut in
!> rank deficient cases. It will be available in the SIGMA library [4].
!> If M is much larger than N, it is obvious that the initial QRF with
!> column pivoting can be preprocessed by the QRF without pivoting. That
!> well known trick is not used in ZGEJSV because in some cases heavy row
!> weighting can be treated with complete pivoting. The overhead in cases
!> M much larger than N is then only due to pivoting, but the benefits in
!> terms of accuracy have prevailed. The implementer/user can incorporate
!> this extra QRF step easily. The implementer can also improve data movement
!> (matrix transpose, matrix copy, matrix transposed copy) - this
!> implementation of ZGEJSV uses only the simplest, naive data movement.
!> !> !> [1] 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. !> [2] 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. !> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR !> factorization software - a case study. !> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. !> LAPACK Working note 176. !> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, !> QSVD, (H,K)-SVD computations. !> Department of Mathematics, University of Zagreb, 2008, 2016. !>
Definition at line 566 of file zgejsv.f.
| subroutine zgesdd | ( | character | jobz, |
| 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( ldvt, * ) | vt, | ||
| integer | ldvt, | ||
| complex*16, dimension( * ) | work, | ||
| integer | lwork, | ||
| double precision, dimension( * ) | rwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
ZGESDD
Download ZGESDD + dependencies [TGZ] [ZIP] [TXT]
!> !> ZGESDD computes the singular value decomposition (SVD) of a complex !> M-by-N matrix A, optionally computing the left and/or right singular !> vectors, by using divide-and-conquer method. The SVD is written !> !> A = U * SIGMA * conjugate-transpose(V) !> !> where SIGMA is an M-by-N matrix which is zero except for its !> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and !> V is an N-by-N unitary matrix. The diagonal elements of SIGMA !> are the singular values of A; they are real and non-negative, and !> are returned in descending order. The first min(m,n) columns of !> U and V are the left and right singular vectors of A. !> !> Note that the routine returns VT = V**H, not V. !> !> The divide and conquer algorithm makes very mild assumptions about !> floating point arithmetic. It will work on machines with a guard !> digit in add/subtract, or on those binary machines without guard !> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or !> Cray-2. It could conceivably fail on hexadecimal or decimal machines !> without guard digits, but we know of none. !>
| [in] | JOBZ | !> JOBZ is CHARACTER*1 !> Specifies options for computing all or part of the matrix U: !> = 'A': all M columns of U and all N rows of V**H are !> returned in the arrays U and VT; !> = 'S': the first min(M,N) columns of U and the first !> min(M,N) rows of V**H are returned in the arrays U !> and VT; !> = 'O': If M >= N, the first N columns of U are overwritten !> in the array A and all rows of V**H are returned in !> the array VT; !> otherwise, all columns of U are returned in the !> array U and the first M rows of V**H are overwritten !> in the array A; !> = 'N': no columns of U or rows of V**H are computed. !> |
| [in] | M | !> M is INTEGER !> The number of rows of the input matrix A. M >= 0. !> |
| [in] | N | !> N is INTEGER !> The number of columns of the input matrix A. N >= 0. !> |
| [in,out] | A | !> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> On exit, !> if JOBZ = 'O', A is overwritten with the first N columns !> of U (the left singular vectors, stored !> columnwise) if M >= N; !> A is overwritten with the first M rows !> of V**H (the right singular vectors, stored !> rowwise) otherwise. !> if JOBZ .ne. 'O', the contents of A are destroyed. !> |
| [in] | LDA | !> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [out] | S | !> S is DOUBLE PRECISION array, dimension (min(M,N)) !> The singular values of A, sorted so that S(i) >= S(i+1). !> |
| [out] | U | !> U is COMPLEX*16 array, dimension (LDU,UCOL) !> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; !> UCOL = min(M,N) if JOBZ = 'S'. !> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M !> unitary matrix U; !> if JOBZ = 'S', U contains the first min(M,N) columns of U !> (the left singular vectors, stored columnwise); !> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. !> |
| [in] | LDU | !> LDU is INTEGER !> The leading dimension of the array U. LDU >= 1; !> if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. !> |
| [out] | VT | !> VT is COMPLEX*16 array, dimension (LDVT,N) !> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the !> N-by-N unitary matrix V**H; !> if JOBZ = 'S', VT contains the first min(M,N) rows of !> V**H (the right singular vectors, stored rowwise); !> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. !> |
| [in] | LDVT | !> LDVT is INTEGER !> The leading dimension of the array VT. LDVT >= 1; !> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; !> if JOBZ = 'S', LDVT >= min(M,N). !> |
| [out] | WORK | !> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. !> |
| [in] | LWORK | !> LWORK is INTEGER !> The dimension of the array WORK. LWORK >= 1. !> If LWORK = -1, a workspace query is assumed. The optimal !> size for the WORK array is calculated and stored in WORK(1), !> and no other work except argument checking is performed. !> !> Let mx = max(M,N) and mn = min(M,N). !> If JOBZ = 'N', LWORK >= 2*mn + mx. !> If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx. !> If JOBZ = 'S', LWORK >= mn*mn + 3*mn. !> If JOBZ = 'A', LWORK >= mn*mn + 2*mn + mx. !> These are not tight minimums in all cases; see comments inside code. !> For good performance, LWORK should generally be larger; !> a query is recommended. !> |
| [out] | RWORK | !> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) !> Let mx = max(M,N) and mn = min(M,N). !> If JOBZ = 'N', LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn); !> else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn; !> else LRWORK >= max( 5*mn*mn + 5*mn, !> 2*mx*mn + 2*mn*mn + mn ). !> |
| [out] | IWORK | !> IWORK is INTEGER array, dimension (8*min(M,N)) !> |
| [out] | INFO | !> INFO is INTEGER !> < 0: if INFO = -i, the i-th argument had an illegal value. !> = -4: if A had a NAN entry. !> > 0: The updating process of DBDSDC did not converge. !> = 0: successful exit. !> |
Definition at line 225 of file zgesdd.f.
| subroutine zgesvd | ( | character | jobu, |
| character | jobvt, | ||
| 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( ldvt, * ) | vt, | ||
| integer | ldvt, | ||
| complex*16, dimension( * ) | work, | ||
| integer | lwork, | ||
| double precision, dimension( * ) | rwork, | ||
| integer | info ) |
ZGESVD computes the singular value decomposition (SVD) for GE matrices
Download ZGESVD + dependencies [TGZ] [ZIP] [TXT]
!> !> ZGESVD computes the singular value decomposition (SVD) of a complex !> M-by-N matrix A, optionally computing the left and/or right singular !> vectors. The SVD is written !> !> A = U * SIGMA * conjugate-transpose(V) !> !> where SIGMA is an M-by-N matrix which is zero except for its !> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and !> V is an N-by-N unitary matrix. The diagonal elements of SIGMA !> are the singular values of A; they are real and non-negative, and !> are returned in descending order. The first min(m,n) columns of !> U and V are the left and right singular vectors of A. !> !> Note that the routine returns V**H, not V. !>
| [in] | JOBU | !> JOBU is CHARACTER*1 !> Specifies options for computing all or part of the matrix U: !> = 'A': all M columns of U are returned in array U: !> = 'S': the first min(m,n) columns of U (the left singular !> vectors) are returned in the array U; !> = 'O': the first min(m,n) columns of U (the left singular !> vectors) are overwritten on the array A; !> = 'N': no columns of U (no left singular vectors) are !> computed. !> |
| [in] | JOBVT | !> JOBVT is CHARACTER*1 !> Specifies options for computing all or part of the matrix !> V**H: !> = 'A': all N rows of V**H are returned in the array VT; !> = 'S': the first min(m,n) rows of V**H (the right singular !> vectors) are returned in the array VT; !> = 'O': the first min(m,n) rows of V**H (the right singular !> vectors) are overwritten on the array A; !> = 'N': no rows of V**H (no right singular vectors) are !> computed. !> !> JOBVT and JOBU cannot both be 'O'. !> |
| [in] | M | !> M is INTEGER !> The number of rows of the input matrix A. M >= 0. !> |
| [in] | N | !> N is INTEGER !> The number of columns of the input matrix A. N >= 0. !> |
| [in,out] | A | !> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> On exit, !> if JOBU = 'O', A is overwritten with the first min(m,n) !> columns of U (the left singular vectors, !> stored columnwise); !> if JOBVT = 'O', A is overwritten with the first min(m,n) !> rows of V**H (the right singular vectors, !> stored rowwise); !> if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A !> are destroyed. !> |
| [in] | LDA | !> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [out] | S | !> S is DOUBLE PRECISION array, dimension (min(M,N)) !> The singular values of A, sorted so that S(i) >= S(i+1). !> |
| [out] | U | !> U is COMPLEX*16 array, dimension (LDU,UCOL) !> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. !> If JOBU = 'A', U contains the M-by-M unitary matrix U; !> if JOBU = 'S', U contains the first min(m,n) columns of U !> (the left singular vectors, stored columnwise); !> if JOBU = 'N' or 'O', U is not referenced. !> |
| [in] | LDU | !> LDU is INTEGER !> The leading dimension of the array U. LDU >= 1; if !> JOBU = 'S' or 'A', LDU >= M. !> |
| [out] | VT | !> VT is COMPLEX*16 array, dimension (LDVT,N) !> If JOBVT = 'A', VT contains the N-by-N unitary matrix !> V**H; !> if JOBVT = 'S', VT contains the first min(m,n) rows of !> V**H (the right singular vectors, stored rowwise); !> if JOBVT = 'N' or 'O', VT is not referenced. !> |
| [in] | LDVT | !> LDVT is INTEGER !> The leading dimension of the array VT. LDVT >= 1; if !> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). !> |
| [out] | WORK | !> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. !> |
| [in] | LWORK | !> LWORK is INTEGER !> The dimension of the array WORK. !> LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)). !> For good performance, LWORK should generally be larger. !> !> If LWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal size of the WORK array, returns !> this value as the first entry of the WORK array, and no error !> message related to LWORK is issued by XERBLA. !> |
| [out] | RWORK | !> RWORK is DOUBLE PRECISION array, dimension (5*min(M,N)) !> On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the !> unconverged superdiagonal elements of an upper bidiagonal !> matrix B whose diagonal is in S (not necessarily sorted). !> B satisfies A = U * B * VT, so it has the same singular !> values as A, and singular vectors related by U and VT. !> |
| [out] | 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 !> did not converge to zero. See the description of RWORK !> above for details. !> |
Definition at line 212 of file zgesvd.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
Download ZGESVDQ + dependencies [TGZ] [ZIP] [TXT]
!> !> 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. !>
| [in] | 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 = 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)
!> |
| [in] | 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. !> |
| [in] | 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. !> |
| [in] | 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. !> |
| [in] | 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. !> |
| [in] | M | !> M is INTEGER !> The number of rows of the input matrix A. M >= 0. !> |
| [in] | N | !> N is INTEGER !> The number of columns of the input matrix A. M >= N >= 0. !> |
| [in,out] | 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. !> |
| [in] | LDA | !> LDA is INTEGER. !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [out] | S | !> S is DOUBLE PRECISION array of dimension N. !> The singular values of A, ordered so that S(i) >= S(i+1). !> |
| [out] | 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. !> |
| [in] | 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. !> |
| [out] | 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. !> |
| [in] | 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. !> |
| [out] | 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. !> |
| [out] | 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. !> |
| [in] | 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. !> |
| [out] | 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. !> |
| [in,out] | LCWORK | !> LCWORK is INTEGER
!> The dimension of the array CWORK. It is determined as follows:
!> 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, 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.
!> |
| [out] | RWORK | !> RWORK is 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. !> |
| [in] | LRWORK | !> LRWORK is INTEGER. !> 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. !> |
| [out] | 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. !> |
!> !> 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, enlosed in [[double brackets]]. !>
!> Please report all bugs and send interesting examples and/or comments to !> drmac@math.hr. Thank you. !>
!> [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 !>
!> Developed and coded by Zlatko Drmac, Department of Mathematics !> University of Zagreb, Croatia, drmac@math.hr !>
Definition at line 410 of file zgesvdq.f.
| subroutine zgesvdx | ( | character | jobu, |
| character | jobvt, | ||
| character | range, | ||
| integer | m, | ||
| integer | n, | ||
| complex*16, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| double precision | vl, | ||
| double precision | vu, | ||
| integer | il, | ||
| integer | iu, | ||
| integer | ns, | ||
| double precision, dimension( * ) | s, | ||
| complex*16, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| complex*16, dimension( ldvt, * ) | vt, | ||
| integer | ldvt, | ||
| complex*16, dimension( * ) | work, | ||
| integer | lwork, | ||
| double precision, dimension( * ) | rwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
ZGESVDX computes the singular value decomposition (SVD) for GE matrices
Download ZGESVDX + dependencies [TGZ] [ZIP] [TXT]
!> !> ZGESVDX computes the singular value decomposition (SVD) of a complex !> M-by-N matrix A, optionally computing the left and/or right singular !> vectors. The SVD is written !> !> A = U * SIGMA * transpose(V) !> !> where SIGMA is an M-by-N matrix which is zero except for its !> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and !> V is an N-by-N unitary matrix. The diagonal elements of SIGMA !> are the singular values of A; they are real and non-negative, and !> are returned in descending order. The first min(m,n) columns of !> U and V are the left and right singular vectors of A. !> !> ZGESVDX uses an eigenvalue problem for obtaining the SVD, which !> allows for the computation of a subset of singular values and !> vectors. See DBDSVDX for details. !> !> Note that the routine returns V**T, not V. !>
| [in] | JOBU | !> JOBU is CHARACTER*1 !> Specifies options for computing all or part of the matrix U: !> = 'V': the first min(m,n) columns of U (the left singular !> vectors) or as specified by RANGE are returned in !> the array U; !> = 'N': no columns of U (no left singular vectors) are !> computed. !> |
| [in] | JOBVT | !> JOBVT is CHARACTER*1 !> Specifies options for computing all or part of the matrix !> V**T: !> = 'V': the first min(m,n) rows of V**T (the right singular !> vectors) or as specified by RANGE are returned in !> the array VT; !> = 'N': no rows of V**T (no right singular vectors) are !> computed. !> |
| [in] | RANGE | !> RANGE is CHARACTER*1 !> = 'A': all singular values will be found. !> = 'V': all singular values in the half-open interval (VL,VU] !> will be found. !> = 'I': the IL-th through IU-th singular values will be found. !> |
| [in] | M | !> M is INTEGER !> The number of rows of the input matrix A. M >= 0. !> |
| [in] | N | !> N is INTEGER !> The number of columns of the input matrix A. N >= 0. !> |
| [in,out] | A | !> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> On exit, the contents of A are destroyed. !> |
| [in] | LDA | !> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [in] | VL | !> VL is DOUBLE PRECISION !> If RANGE='V', the lower bound of the interval to !> be searched for singular values. VU > VL. !> Not referenced if RANGE = 'A' or 'I'. !> |
| [in] | VU | !> VU is DOUBLE PRECISION !> If RANGE='V', the upper bound of the interval to !> be searched for singular values. VU > VL. !> Not referenced if RANGE = 'A' or 'I'. !> |
| [in] | IL | !> IL is INTEGER !> If RANGE='I', the index of the !> smallest singular value to be returned. !> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. !> Not referenced if RANGE = 'A' or 'V'. !> |
| [in] | IU | !> IU is INTEGER !> If RANGE='I', the index of the !> largest singular value to be returned. !> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. !> Not referenced if RANGE = 'A' or 'V'. !> |
| [out] | NS | !> NS is INTEGER !> The total number of singular values found, !> 0 <= NS <= min(M,N). !> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1. !> |
| [out] | S | !> S is DOUBLE PRECISION array, dimension (min(M,N)) !> The singular values of A, sorted so that S(i) >= S(i+1). !> |
| [out] | U | !> U is COMPLEX*16 array, dimension (LDU,UCOL) !> If JOBU = 'V', U contains columns of U (the left singular !> vectors, stored columnwise) as specified by RANGE; if !> JOBU = 'N', U is not referenced. !> Note: The user must ensure that UCOL >= NS; if RANGE = 'V', !> the exact value of NS is not known in advance and an upper !> bound must be used. !> |
| [in] | LDU | !> LDU is INTEGER !> The leading dimension of the array U. LDU >= 1; if !> JOBU = 'V', LDU >= M. !> |
| [out] | VT | !> VT is COMPLEX*16 array, dimension (LDVT,N) !> If JOBVT = 'V', VT contains the rows of V**T (the right singular !> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N', !> VT is not referenced. !> Note: The user must ensure that LDVT >= NS; if RANGE = 'V', !> the exact value of NS is not known in advance and an upper !> bound must be used. !> |
| [in] | LDVT | !> LDVT is INTEGER !> The leading dimension of the array VT. LDVT >= 1; if !> JOBVT = 'V', LDVT >= NS (see above). !> |
| [out] | WORK | !> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK; !> |
| [in] | LWORK | !> LWORK is INTEGER !> The dimension of the array WORK. !> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see !> comments inside the code): !> - PATH 1 (M much larger than N) !> - PATH 1t (N much larger than M) !> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths. !> For good performance, LWORK should generally be larger. !> !> If LWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal size of the WORK array, returns !> this value as the first entry of the WORK array, and no error !> message related to LWORK is issued by XERBLA. !> |
| [out] | RWORK | !> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) !> LRWORK >= MIN(M,N)*(MIN(M,N)*2+15*MIN(M,N)). !> |
| [out] | IWORK | !> IWORK is INTEGER array, dimension (12*MIN(M,N)) !> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0, !> then IWORK contains the indices of the eigenvectors that failed !> to converge in DBDSVDX/DSTEVX. !> |
| [out] | INFO | !> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !> > 0: if INFO = i, then i eigenvectors failed to converge !> in DBDSVDX/DSTEVX. !> if INFO = N*2 + 1, an internal error occurred in !> DBDSVDX !> |
Definition at line 267 of file zgesvdx.f.
| subroutine zggsvd3 | ( | character | jobu, |
| character | jobv, | ||
| character | jobq, | ||
| integer | m, | ||
| integer | n, | ||
| integer | p, | ||
| integer | k, | ||
| integer | l, | ||
| complex*16, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| complex*16, dimension( ldb, * ) | b, | ||
| integer | ldb, | ||
| double precision, dimension( * ) | alpha, | ||
| double precision, dimension( * ) | beta, | ||
| complex*16, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| complex*16, dimension( ldv, * ) | v, | ||
| integer | ldv, | ||
| complex*16, dimension( ldq, * ) | q, | ||
| integer | ldq, | ||
| complex*16, dimension( * ) | work, | ||
| integer | lwork, | ||
| double precision, dimension( * ) | rwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
ZGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Download ZGGSVD3 + dependencies [TGZ] [ZIP] [TXT]
!> !> ZGGSVD3 computes the generalized singular value decomposition (GSVD) !> of an M-by-N complex matrix A and P-by-N complex matrix B: !> !> U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R ) !> !> where U, V and Q are unitary matrices. !> Let K+L = the effective numerical rank of the !> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper !> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) !> matrices and of the following structures, respectively: !> !> If M-K-L >= 0, !> !> K L !> D1 = K ( I 0 ) !> L ( 0 C ) !> M-K-L ( 0 0 ) !> !> K L !> D2 = L ( 0 S ) !> P-L ( 0 0 ) !> !> N-K-L K L !> ( 0 R ) = K ( 0 R11 R12 ) !> L ( 0 0 R22 ) !> where !> !> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), !> S = diag( BETA(K+1), ... , BETA(K+L) ), !> C**2 + S**2 = I. !> !> R is stored in A(1:K+L,N-K-L+1:N) on exit. !> !> If M-K-L < 0, !> !> K M-K K+L-M !> D1 = K ( I 0 0 ) !> M-K ( 0 C 0 ) !> !> K M-K K+L-M !> D2 = M-K ( 0 S 0 ) !> K+L-M ( 0 0 I ) !> P-L ( 0 0 0 ) !> !> N-K-L K M-K K+L-M !> ( 0 R ) = K ( 0 R11 R12 R13 ) !> M-K ( 0 0 R22 R23 ) !> K+L-M ( 0 0 0 R33 ) !> !> where !> !> C = diag( ALPHA(K+1), ... , ALPHA(M) ), !> S = diag( BETA(K+1), ... , BETA(M) ), !> C**2 + S**2 = I. !> !> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored !> ( 0 R22 R23 ) !> in B(M-K+1:L,N+M-K-L+1:N) on exit. !> !> The routine computes C, S, R, and optionally the unitary !> transformation matrices U, V and Q. !> !> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of !> A and B implicitly gives the SVD of A*inv(B): !> A*inv(B) = U*(D1*inv(D2))*V**H. !> If ( A**H,B**H)**H has orthonormal columns, then the GSVD of A and B is also !> equal to the CS decomposition of A and B. Furthermore, the GSVD can !> be used to derive the solution of the eigenvalue problem: !> A**H*A x = lambda* B**H*B x. !> In some literature, the GSVD of A and B is presented in the form !> U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 ) !> where U and V are orthogonal and X is nonsingular, and D1 and D2 are !> ``diagonal''. The former GSVD form can be converted to the latter !> form by taking the nonsingular matrix X as !> !> X = Q*( I 0 ) !> ( 0 inv(R) ) !>
| [in] | JOBU | !> JOBU is CHARACTER*1 !> = 'U': Unitary matrix U is computed; !> = 'N': U is not computed. !> |
| [in] | JOBV | !> JOBV is CHARACTER*1 !> = 'V': Unitary matrix V is computed; !> = 'N': V is not computed. !> |
| [in] | JOBQ | !> JOBQ is CHARACTER*1 !> = 'Q': Unitary matrix Q is computed; !> = 'N': Q is not computed. !> |
| [in] | M | !> M is INTEGER !> The number of rows of the matrix A. M >= 0. !> |
| [in] | N | !> N is INTEGER !> The number of columns of the matrices A and B. N >= 0. !> |
| [in] | P | !> P is INTEGER !> The number of rows of the matrix B. P >= 0. !> |
| [out] | K | !> K is INTEGER !> |
| [out] | L | !> L is INTEGER !> !> On exit, K and L specify the dimension of the subblocks !> described in Purpose. !> K + L = effective numerical rank of (A**H,B**H)**H. !> |
| [in,out] | A | !> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> On exit, A contains the triangular matrix R, or part of R. !> See Purpose for details. !> |
| [in] | LDA | !> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [in,out] | B | !> B is COMPLEX*16 array, dimension (LDB,N) !> On entry, the P-by-N matrix B. !> On exit, B contains part of the triangular matrix R if !> M-K-L < 0. See Purpose for details. !> |
| [in] | LDB | !> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,P). !> |
| [out] | ALPHA | !> ALPHA is DOUBLE PRECISION array, dimension (N) !> |
| [out] | BETA | !> BETA is DOUBLE PRECISION array, dimension (N) !> !> On exit, ALPHA and BETA contain the generalized singular !> value pairs of A and B; !> ALPHA(1:K) = 1, !> BETA(1:K) = 0, !> and if M-K-L >= 0, !> ALPHA(K+1:K+L) = C, !> BETA(K+1:K+L) = S, !> or if M-K-L < 0, !> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 !> BETA(K+1:M) =S, BETA(M+1:K+L) =1 !> and !> ALPHA(K+L+1:N) = 0 !> BETA(K+L+1:N) = 0 !> |
| [out] | U | !> U is COMPLEX*16 array, dimension (LDU,M) !> If JOBU = 'U', U contains the M-by-M unitary matrix U. !> If JOBU = 'N', U is not referenced. !> |
| [in] | LDU | !> LDU is INTEGER !> The leading dimension of the array U. LDU >= max(1,M) if !> JOBU = 'U'; LDU >= 1 otherwise. !> |
| [out] | V | !> V is COMPLEX*16 array, dimension (LDV,P) !> If JOBV = 'V', V contains the P-by-P unitary matrix V. !> If JOBV = 'N', V is not referenced. !> |
| [in] | LDV | !> LDV is INTEGER !> The leading dimension of the array V. LDV >= max(1,P) if !> JOBV = 'V'; LDV >= 1 otherwise. !> |
| [out] | Q | !> Q is COMPLEX*16 array, dimension (LDQ,N) !> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. !> If JOBQ = 'N', Q is not referenced. !> |
| [in] | LDQ | !> LDQ is INTEGER !> The leading dimension of the array Q. LDQ >= max(1,N) if !> JOBQ = 'Q'; LDQ >= 1 otherwise. !> |
| [out] | WORK | !> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. !> |
| [in] | LWORK | !> LWORK is INTEGER !> The dimension of the array WORK. !> !> If LWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal size of the WORK array, returns !> this value as the first entry of the WORK array, and no error !> message related to LWORK is issued by XERBLA. !> |
| [out] | RWORK | !> RWORK is DOUBLE PRECISION array, dimension (2*N) !> |
| [out] | IWORK | !> IWORK is INTEGER array, dimension (N) !> On exit, IWORK stores the sorting information. More !> precisely, the following loop will sort ALPHA !> for I = K+1, min(M,K+L) !> swap ALPHA(I) and ALPHA(IWORK(I)) !> endfor !> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). !> |
| [out] | INFO | !> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: if INFO = 1, the Jacobi-type procedure failed to !> converge. For further details, see subroutine ZTGSJA. !> |
!> TOLA DOUBLE PRECISION !> TOLB DOUBLE PRECISION !> TOLA and TOLB are the thresholds to determine the effective !> rank of (A**H,B**H)**H. Generally, they are set to !> TOLA = MAX(M,N)*norm(A)*MACHEPS, !> TOLB = MAX(P,N)*norm(B)*MACHEPS. !> The size of TOLA and TOLB may affect the size of backward !> errors of the decomposition. !>
Definition at line 350 of file zggsvd3.f.
1.15.0