Functions | |
| subroutine | sgejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, a, lda, sva, u, ldu, v, ldv, work, lwork, iwork, info) |
| SGEJSV | |
| subroutine | sgesdd (jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info) |
| SGESDD | |
| subroutine | sgesvd (jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info) |
| SGESVD computes the singular value decomposition (SVD) for GE matrices | |
| subroutine | sgesvdq (joba, jobp, jobr, jobu, jobv, m, n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work, lwork, rwork, lrwork, info) |
| SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices | |
| subroutine | sgesvdx (jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, iwork, info) |
| SGESVDX computes the singular value decomposition (SVD) for GE matrices | |
| subroutine | sggsvd3 (jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, lwork, iwork, info) |
| SGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices | |
This is the group of real singular value driver functions for GE matrices
| subroutine sgejsv | ( | character*1 | joba, |
| character*1 | jobu, | ||
| character*1 | jobv, | ||
| character*1 | jobr, | ||
| character*1 | jobt, | ||
| character*1 | jobp, | ||
| integer | m, | ||
| integer | n, | ||
| real, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| real, dimension( n ) | sva, | ||
| real, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| real, dimension( ldv, * ) | v, | ||
| integer | ldv, | ||
| real, dimension( lwork ) | work, | ||
| integer | lwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
SGEJSV
Download SGEJSV + dependencies [TGZ] [ZIP] [TXT]
!> !> SGEJSV computes the singular value decomposition (SVD) of a real M-by-N !> matrix [A], where M >= N. The SVD of [A] is written as !> !> [A] = [U] * [SIGMA] * [V]^t, !> !> 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) 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. 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. !> SGEJSV can sometimes compute tiny singular values and their singular vectors much !> more accurately than other SVD routines, see below under Further Details. !>
| [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=D*B. If A has heavily weighted
!> rows, then using this condition number gives too pessimistic
!> error bound.
!> = 'A': Small singular values are the noise and 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^t 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. This option is !> allowed only if JOBU .NE. 'N', i.e. in computing the full SVD. !> = '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=SLAMCH('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=SLAMCH('S'), EPSLN=SLAMCH('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 SGESVJ.
!> = '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 SGESVJ.
!> |
| [in] | JOBT | !> JOBT is CHARACTER*1 !> If the matrix is square then the procedure may determine to use !> transposed A if A^t seems to be better with respect to convergence. !> If the matrix is not square, JOBT is ignored. This is subject to !> changes in the future. !> The decision is based on two values of entropy over the adjoint !> orbit of A^t * A. See the descriptions of WORK(6) and WORK(7). !> = 'T': transpose if entropy test indicates possibly faster !> convergence of Jacobi process if A^t is taken as input. If A is !> replaced with A^t, then the row pivoting is included automatically. !> = 'N': do not speculate. !> This option 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] | 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 REAL 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 REAL 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 REAL 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^t. In that case, [V] is computed !> in U as left singular vectors of A^t 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 REAL 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^t. In that case, [U] is computed !> in V as right singular vectors of A^t 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] | WORK | !> WORK is REAL array, dimension (LWORK) !> On exit, !> WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such !> that SCALE*SVA(1:N) are the computed singular values !> of A. (See the description of SVA().) !> WORK(2) = See the description of WORK(1). !> WORK(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^t * R)^(-1)||_1). !> It is computed using SPOCON. 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. !> !> WORK(4) = an estimate of the scaled condition number of the !> triangular factor in the first QR factorization. !> WORK(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. !> !> WORK(6) = the entropy of A^t*A :: this is the Shannon entropy !> of diag(A^t*A) / Trace(A^t*A) taken as point in the !> probability simplex. !> WORK(7) = the entropy of A*A^t. !> |
| [in] | LWORK | !> LWORK is INTEGER !> Length of WORK to confirm proper allocation of work space. !> LWORK depends on the job: !> !> If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and !> -> .. no scaled condition estimate required (JOBE = 'N'): !> LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement. !> ->> For optimal performance (blocked code) the optimal value !> is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal !> block size for DGEQP3 and DGEQRF. !> In general, optimal LWORK is computed as !> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7). !> -> .. an estimate of the scaled condition number of A is !> required (JOBA='E', 'G'). In this case, LWORK is the maximum !> of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7). !> ->> For optimal performance (blocked code) the optimal value !> is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7). !> In general, the optimal length LWORK is computed as !> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), !> N+N*N+LWORK(DPOCON),7). !> !> If SIGMA and the right singular vectors are needed (JOBV = 'V'), !> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). !> -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7), !> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ, !> DORMLQ. In general, the optimal length LWORK is computed as !> LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), !> N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)). !> !> If SIGMA and the left singular vectors are needed !> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). !> -> For optimal performance: !> if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7), !> if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7), !> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR. !> In general, the optimal length LWORK is computed as !> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON), !> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). !> Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or !> M*NB (for JOBU = 'F'). !> !> If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and !> -> if JOBV = 'V' !> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). !> -> if JOBV = 'J' the minimal requirement is !> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6). !> -> For optimal performance, LWORK should be additionally !> larger than N+M*NB, where NB is the optimal block size !> for DORMQR. !> |
| [out] | IWORK | !> IWORK is INTEGER array, dimension (M+3*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. !> |
| [out] | INFO | !> INFO is INTEGER !> < 0: if INFO = -i, then the i-th argument had an illegal value. !> = 0: successful exit; !> > 0: SGEJSV did not converge in the maximal allowed number !> of sweeps. The computed values may be inaccurate. !> |
!>
!> SGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3,
!> SGEQRF, and SGELQF 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 SGEJSV 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 (SGEJSV) is best used in this restricted range,
!> meaning that singular values of magnitude below ||A||_2 / SLAMCH('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 (SGESVJ) is
!> left to the implementer on a particular machine.
!> The rank revealing QR factorization (in this code: SGEQP3) should be
!> implemented as in [3]. We have a new version of SGEQP3 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 SGEJSV 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 SGEJSV 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. !>
Definition at line 473 of file sgejsv.f.
| subroutine sgesdd | ( | character | jobz, |
| integer | m, | ||
| integer | n, | ||
| real, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| real, dimension( * ) | s, | ||
| real, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| real, dimension( ldvt, * ) | vt, | ||
| integer | ldvt, | ||
| real, dimension( * ) | work, | ||
| integer | lwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
SGESDD
Download SGESDD + dependencies [TGZ] [ZIP] [TXT]
!> !> SGESDD computes the singular value decomposition (SVD) of a real !> M-by-N matrix A, optionally computing the left and right singular !> vectors. If singular vectors are desired, it uses a !> divide-and-conquer algorithm. !> !> 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 orthogonal matrix, and !> V is an N-by-N orthogonal 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**T, 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**T 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**T are returned in the arrays U !> and VT; !> = 'O': If M >= N, the first N columns of U are overwritten !> on the array A and all rows of V**T are returned in !> the array VT; !> otherwise, all columns of U are returned in the !> array U and the first M rows of V**T are overwritten !> in the array A; !> = 'N': no columns of U or rows of V**T 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 REAL 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**T (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 REAL array, dimension (min(M,N)) !> The singular values of A, sorted so that S(i) >= S(i+1). !> |
| [out] | U | !> U is REAL 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 !> orthogonal 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 REAL array, dimension (LDVT,N) !> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the !> N-by-N orthogonal matrix V**T; !> if JOBZ = 'S', VT contains the first min(M,N) rows of !> V**T (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 REAL 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 >= 3*mn + max( mx, 7*mn ). !> If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ). !> If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn. !> If JOBZ = 'A', LWORK >= 4*mn*mn + 6*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] | 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: SBDSDC did not converge, updating process failed. !> = 0: successful exit. !> |
Definition at line 217 of file sgesdd.f.
| subroutine sgesvd | ( | character | jobu, |
| character | jobvt, | ||
| integer | m, | ||
| integer | n, | ||
| real, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| real, dimension( * ) | s, | ||
| real, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| real, dimension( ldvt, * ) | vt, | ||
| integer | ldvt, | ||
| real, dimension( * ) | work, | ||
| integer | lwork, | ||
| integer | info ) |
SGESVD computes the singular value decomposition (SVD) for GE matrices
Download SGESVD + dependencies [TGZ] [ZIP] [TXT]
!> !> SGESVD computes the singular value decomposition (SVD) of a real !> 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 orthogonal matrix, and !> V is an N-by-N orthogonal 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**T, 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**T: !> = 'A': all N rows of V**T are returned in the array VT; !> = 'S': the first min(m,n) rows of V**T (the right singular !> vectors) are returned in the array VT; !> = 'O': the first min(m,n) rows of V**T (the right singular !> vectors) are overwritten on the array A; !> = 'N': no rows of V**T (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 REAL 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**T (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 REAL array, dimension (min(M,N)) !> The singular values of A, sorted so that S(i) >= S(i+1). !> |
| [out] | U | !> U is REAL 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 orthogonal 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 REAL array, dimension (LDVT,N) !> If JOBVT = 'A', VT contains the N-by-N orthogonal matrix !> V**T; !> if JOBVT = 'S', VT contains the first min(m,n) rows of !> V**T (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 REAL array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK; !> if INFO > 0, WORK(2:MIN(M,N)) 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. !> |
| [in] | LWORK | !> LWORK is INTEGER !> The dimension of the array WORK. !> LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code): !> - PATH 1 (M much larger than N, JOBU='N') !> - PATH 1t (N much larger than M, JOBVT='N') !> LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(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] | INFO | !> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: if SBDSQR did not converge, INFO specifies how many !> superdiagonals of an intermediate bidiagonal form B !> did not converge to zero. See the description of WORK !> above for details. !> |
Definition at line 209 of file sgesvd.f.
| subroutine sgesvdq | ( | character | joba, |
| character | jobp, | ||
| character | jobr, | ||
| character | jobu, | ||
| character | jobv, | ||
| integer | m, | ||
| integer | n, | ||
| real, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| real, dimension( * ) | s, | ||
| real, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| real, dimension( ldv, * ) | v, | ||
| integer | ldv, | ||
| integer | numrank, | ||
| integer, dimension( * ) | iwork, | ||
| integer | liwork, | ||
| real, dimension( * ) | work, | ||
| integer | lwork, | ||
| real, dimension( * ) | rwork, | ||
| integer | lrwork, | ||
| integer | info ) |
SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
Download SGESVDQ + dependencies [TGZ] [ZIP] [TXT]
!> !> SGESVDQ computes the singular value decomposition (SVD) of a real !> M-by-N matrix A, where M >= N. The SVD of A is written as !> [++] [xx] [x0] [xx] !> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx] !> [++] [xx] !> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal !> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements !> of SIGMA are the singular values of A. The columns of U and V are the !> left and the right singular vectors of A, respectively. !>
| [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 = SLAMCH('Epsilon'). This authorises CGESVDQ to
!> truncate the computed triangular factor in a rank revealing
!> QR factorization whenever the truncated part is below the
!> threshold of the order of EPS * ||A||_F. This is aggressive
!> truncation level.
!> = 'M' Similarly as with 'A', but the truncation is more gentle: it
!> is allowed only when there is a drop on the diagonal of the
!> triangular factor in the QR factorization. This is medium
!> truncation level.
!> = 'H' High accuracy requested. No numerical rank determination based
!> on the rank revealing QR factorization is attempted.
!> = 'E' Same as 'H', and in addition the condition number of column
!> scaled A is estimated and returned in RWORK(1).
!> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
!> |
| [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, SGESVD is applied to !> the transposed R**T of the computed triangular factor R. This involves !> some extra data movement (matrix transpositions). Useful for !> experiments, research and development. !> = 'N' The triangular factor R is given as input to SGESVD. This may be !> preferred as it involves less data movement. !> |
| [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**T , 0)**T. If row pivoting is used, !> then the necessary information on the row pivoting is stored in !> IWORK(N+1:N+M-1). !> = 'N' The left singular vectors are not computed. !> |
| [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 REAL array of dimensions LDA x N !> On entry, the input matrix A. !> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains !> the Householder vectors as stored by SGEQP3. If JOBU = 'F', these Householder !> vectors together with WORK(1:N) can be used to restore the Q factors from !> the initial pivoted QR factorization of A. See the description of U. !> |
| [in] | LDA | !> LDA is INTEGER. !> The leading dimension of the array A. LDA >= max(1,M). !> |
| [out] | S | !> S is REAL array of dimension N. !> The singular values of A, ordered so that S(i) >= S(i+1). !> |
| [out] | U | !> U is REAL array, dimension !> LDU x M if JOBU = 'A'; see the description of LDU. In this case, !> on exit, U contains the M left singular vectors. !> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this !> case, U contains the leading N or the leading NUMRANK left singular vectors. !> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U !> contains N x N orthogonal matrix that can be used to form the left !> singular vectors. !> If JOBU = 'N', U is not referenced. !> |
| [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 REAL array, dimension !> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' . !> If JOBV = 'A', or 'V', V contains the N-by-N orthogonal matrix V**T; !> If JOBV = 'R', V contains the first NUMRANK rows of V**T (the right !> singular vectors, stored rowwise, of the NUMRANK largest singular values). !> If JOBV = 'N' and JOBA = 'E', V is used as a workspace. !> If JOBV = 'N', and JOBA.NE.'E', V is not referenced. !> |
| [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 SGESVD. 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, LWORK, 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' and JOBA .NE. 'E'; !> LIWORK >= N if JOBP = 'N' and JOBA .NE. 'E'; !> LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E'; !> LIWORK >= N + N if JOBP = 'N' and JOBA = 'E'. !> !> If LIWORK = -1, then a workspace query is assumed; the routine !> only calculates and returns the optimal and minimal sizes !> for the WORK, IWORK, and RWORK arrays, and no error !> message related to LWORK is issued by XERBLA. !> |
| [out] | WORK | !> WORK is REAL array, dimension (max(2, LWORK)), used as a workspace. !> On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters !> needed to recover the Q factor from the QR factorization computed by !> SGEQP3. !> !> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, !> WORK(1) returns the optimal LWORK, and !> WORK(2) returns the minimal LWORK. !> |
| [in,out] | LWORK | !> LWORK is INTEGER
!> The dimension of the array WORK. It is determined as follows:
!> Let LWQP3 = 3*N+1, LWCON = 3*N, and let
!> LWORQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U'
!> { MAX( M, 1 ), if JOBU = 'A'
!> LWSVD = MAX( 5*N, 1 )
!> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
!> LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
!> Then the minimal value of LWORK is:
!> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed;
!> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
!> and a scaled condition estimate requested;
!>
!> = N + MAX( LWQP3, LWSVD, LWORQ ) if the singular values and the left
!> singular vectors are requested;
!> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular values and the left
!> singular vectors are requested, and also
!> a scaled condition estimate requested;
!>
!> = N + MAX( LWQP3, LWSVD ) if the singular values and the right
!> singular vectors are requested;
!> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
!> singular vectors are requested, and also
!> a scaled condition etimate requested;
!>
!> = N + MAX( LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV = 'R';
!> independent of JOBR;
!> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
!> JOBV = 'R' and, also a scaled condition
!> estimate requested; independent of JOBR;
!> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
!> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
!> full SVD is requested with JOBV = 'A' or 'V', and
!> JOBR ='N'
!> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
!> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
!> if the full SVD is requested with JOBV = 'A' or 'V', and
!> JOBR ='N', and also a scaled condition number estimate
!> requested.
!> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
!> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
!> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
!> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
!> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
!> if the full SVD is requested with JOBV = 'A' or 'V', and
!> JOBR ='T', and also a scaled condition number estimate
!> requested.
!> Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).
!>
!> If LWORK = -1, then a workspace query is assumed; the routine
!> only calculates and returns the optimal and minimal sizes
!> for the WORK, IWORK, and RWORK arrays, and no error
!> message related to LWORK is issued by XERBLA.
!> |
| [out] | RWORK | !> RWORK is REAL array, dimension (max(1, LRWORK)). !> On exit, !> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition !> number of column scaled A. If A = C * D where D is diagonal and C !> has unit columns in the Euclidean norm, then, assuming full column rank, !> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1). !> Otherwise, RWORK(1) = -1. !> 2. RWORK(2) contains the number of singular values computed as !> exact zeros in SGESVD applied to the upper triangular or trapezoidal !> R (from the initial QR factorization). In case of early exit (no call to !> SGESVD, such as in the case of zero matrix) RWORK(2) = -1. !> !> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, !> RWORK(1) returns the minimal LRWORK. !> |
| [in] | LRWORK | !> LRWORK is INTEGER. !> The dimension of the array RWORK. !> If JOBP ='P', then LRWORK >= MAX(2, M). !> Otherwise, LRWORK >= 2 !> !> If LRWORK = -1, then a workspace query is assumed; the routine !> only calculates and returns the optimal and minimal sizes !> for the WORK, IWORK, and RWORK arrays, and no error !> message related to LWORK is issued by XERBLA. !> |
| [out] | INFO | !> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: if SBDSQR did not converge, INFO specifies how many superdiagonals !> of an intermediate bidiagonal form B (computed in SGESVD) did not !> converge to zero. !> |
!> !> 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 412 of file sgesvdq.f.
| subroutine sgesvdx | ( | character | jobu, |
| character | jobvt, | ||
| character | range, | ||
| integer | m, | ||
| integer | n, | ||
| real, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| real | vl, | ||
| real | vu, | ||
| integer | il, | ||
| integer | iu, | ||
| integer | ns, | ||
| real, dimension( * ) | s, | ||
| real, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| real, dimension( ldvt, * ) | vt, | ||
| integer | ldvt, | ||
| real, dimension( * ) | work, | ||
| integer | lwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
SGESVDX computes the singular value decomposition (SVD) for GE matrices
Download SGESVDX + dependencies [TGZ] [ZIP] [TXT]
!> !> SGESVDX computes the singular value decomposition (SVD) of a real !> 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 orthogonal matrix, and !> V is an N-by-N orthogonal 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. !> !> SGESVDX uses an eigenvalue problem for obtaining the SVD, which !> allows for the computation of a subset of singular values and !> vectors. See SBDSVDX 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 REAL 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 REAL !> 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 REAL !> 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 REAL array, dimension (min(M,N)) !> The singular values of A, sorted so that S(i) >= S(i+1). !> |
| [out] | U | !> U is REAL 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 REAL 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 REAL 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] | 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 SBDSVDX/SSTEVX. !> |
| [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 SBDSVDX/SSTEVX. !> if INFO = N*2 + 1, an internal error occurred in !> SBDSVDX !> |
Definition at line 260 of file sgesvdx.f.
| subroutine sggsvd3 | ( | character | jobu, |
| character | jobv, | ||
| character | jobq, | ||
| integer | m, | ||
| integer | n, | ||
| integer | p, | ||
| integer | k, | ||
| integer | l, | ||
| real, dimension( lda, * ) | a, | ||
| integer | lda, | ||
| real, dimension( ldb, * ) | b, | ||
| integer | ldb, | ||
| real, dimension( * ) | alpha, | ||
| real, dimension( * ) | beta, | ||
| real, dimension( ldu, * ) | u, | ||
| integer | ldu, | ||
| real, dimension( ldv, * ) | v, | ||
| integer | ldv, | ||
| real, dimension( ldq, * ) | q, | ||
| integer | ldq, | ||
| real, dimension( * ) | work, | ||
| integer | lwork, | ||
| integer, dimension( * ) | iwork, | ||
| integer | info ) |
SGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
Download SGGSVD3 + dependencies [TGZ] [ZIP] [TXT]
!> !> SGGSVD3 computes the generalized singular value decomposition (GSVD) !> of an M-by-N real matrix A and P-by-N real matrix B: !> !> U**T*A*Q = D1*( 0 R ), V**T*B*Q = D2*( 0 R ) !> !> where U, V and Q are orthogonal matrices. !> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T, !> 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 orthogonal !> 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**T. !> If ( A**T,B**T)**T 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**T*A x = lambda* B**T*B x. !> In some literature, the GSVD of A and B is presented in the form !> U**T*A*X = ( 0 D1 ), V**T*B*X = ( 0 D2 ) !> where U and V are orthogonal and X is nonsingular, 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': Orthogonal matrix U is computed; !> = 'N': U is not computed. !> |
| [in] | JOBV | !> JOBV is CHARACTER*1 !> = 'V': Orthogonal matrix V is computed; !> = 'N': V is not computed. !> |
| [in] | JOBQ | !> JOBQ is CHARACTER*1 !> = 'Q': Orthogonal 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**T,B**T)**T. !> |
| [in,out] | A | !> A is REAL 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 REAL array, dimension (LDB,N) !> On entry, the P-by-N matrix B. !> On exit, B contains 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 REAL array, dimension (N) !> |
| [out] | BETA | !> BETA is REAL 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 REAL array, dimension (LDU,M) !> If JOBU = 'U', U contains the M-by-M orthogonal 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 REAL array, dimension (LDV,P) !> If JOBV = 'V', V contains the P-by-P orthogonal 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 REAL array, dimension (LDQ,N) !> If JOBQ = 'Q', Q contains the N-by-N orthogonal 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 REAL 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] | 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 STGSJA. !> |
!> TOLA REAL !> TOLB REAL !> TOLA and TOLB are the thresholds to determine the effective !> rank of (A**T,B**T)**T. 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 346 of file sggsvd3.f.