|
version 1.3, 2016/08/27 15:34:45
|
version 1.4, 2017/06/17 10:54:09
|
|
Line 1
|
Line 1
|
| *> \brief \b ZGEJSV |
*> \brief \b ZGEJSV
|
| * |
*
|
| * =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION ===========
|
| * |
*
|
| * Online html documentation available at |
* Online html documentation available at
|
| * http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/
|
| * |
*
|
| *> \htmlonly |
*> \htmlonly
|
| *> Download ZGEJSV + dependencies |
*> Download ZGEJSV + dependencies
|
| *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgejsv.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgejsv.f">
|
| *> [TGZ]</a> |
*> [TGZ]</a>
|
| *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgejsv.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgejsv.f">
|
| *> [ZIP]</a> |
*> [ZIP]</a>
|
| *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgejsv.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgejsv.f">
|
| *> [TXT]</a> |
*> [TXT]</a>
|
| *> \endhtmlonly |
*> \endhtmlonly
|
| * |
*
|
| * Definition: |
* Definition:
|
| * =========== |
* ===========
|
| * |
*
|
| * SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, |
* SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
|
| * M, N, A, LDA, SVA, U, LDU, V, LDV, |
* M, N, A, LDA, SVA, U, LDU, V, LDV,
|
| * CWORK, LWORK, RWORK, LRWORK, IWORK, INFO ) |
* CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
|
| * |
*
|
| * .. Scalar Arguments .. |
* .. Scalar Arguments ..
|
| * IMPLICIT NONE |
* IMPLICIT NONE
|
| * INTEGER INFO, LDA, LDU, LDV, LWORK, M, N |
* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
|
| * .. |
* ..
|
| * .. Array Arguments .. |
* .. Array Arguments ..
|
| * COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( LWORK ) |
* COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( LWORK )
|
| * DOUBLE PRECISION SVA( N ), RWORK( LRWORK ) |
* DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
|
| * INTEGER IWORK( * ) |
* INTEGER IWORK( * )
|
| * CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV |
* CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
|
| * .. |
* ..
|
| * |
*
|
| * |
*
|
| *> \par Purpose: |
*> \par Purpose:
|
| * ============= |
* =============
|
| *> |
*>
|
| *> \verbatim |
*> \verbatim
|
| *> |
*>
|
| *> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N |
*> 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 |
*> matrix [A], where M >= N. The SVD of [A] is written as
|
| *> |
*>
|
| *> [A] = [U] * [SIGMA] * [V]^*, |
*> [A] = [U] * [SIGMA] * [V]^*,
|
| *> |
*>
|
| *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N |
*> 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 |
*> 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 |
*> [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 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] |
*> 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 |
*> are computed and stored in the arrays U and V, respectively. The diagonal
|
| *> of [SIGMA] is computed and stored in the array SVA. |
*> of [SIGMA] is computed and stored in the array SVA.
|
| *> \endverbatim |
*> \endverbatim
|
| *> |
*>
|
| *> Arguments: |
*> Arguments:
|
| *> ========== |
*> ==========
|
| *> |
*>
|
| *> \param[in] JOBA |
*> \param[in] JOBA
|
| *> \verbatim |
*> \verbatim
|
| *> JOBA is CHARACTER*1 |
*> JOBA is CHARACTER*1
|
| *> Specifies the level of accuracy: |
*> Specifies the level of accuracy:
|
| *> = 'C': This option works well (high relative accuracy) if A = B * D, |
*> = 'C': This option works well (high relative accuracy) if A = B * D,
|
| *> with well-conditioned B and arbitrary diagonal matrix D. |
*> with well-conditioned B and arbitrary diagonal matrix D.
|
| *> The accuracy cannot be spoiled by COLUMN scaling. The |
*> The accuracy cannot be spoiled by COLUMN scaling. The
|
| *> accuracy of the computed output depends on the condition of |
*> accuracy of the computed output depends on the condition of
|
| *> B, and the procedure aims at the best theoretical accuracy. |
*> B, and the procedure aims at the best theoretical accuracy.
|
| *> The relative error max_{i=1:N}|d sigma_i| / sigma_i is |
*> The relative error max_{i=1:N}|d sigma_i| / sigma_i is
|
| *> bounded by f(M,N)*epsilon* cond(B), independent of D. |
*> bounded by f(M,N)*epsilon* cond(B), independent of D.
|
| *> The input matrix is preprocessed with the QRF with column |
*> The input matrix is preprocessed with the QRF with column
|
| *> pivoting. This initial preprocessing and preconditioning by |
*> pivoting. This initial preprocessing and preconditioning by
|
| *> a rank revealing QR factorization is common for all values of |
*> a rank revealing QR factorization is common for all values of
|
| *> JOBA. Additional actions are specified as follows: |
*> JOBA. Additional actions are specified as follows:
|
| *> = 'E': Computation as with 'C' with an additional estimate of the |
*> = 'E': Computation as with 'C' with an additional estimate of the
|
| *> condition number of B. It provides a realistic error bound. |
*> condition number of B. It provides a realistic error bound.
|
| *> = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings |
*> = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
|
| *> D1, D2, and well-conditioned matrix C, this option gives |
*> D1, D2, and well-conditioned matrix C, this option gives
|
| *> higher accuracy than the 'C' option. If the structure of the |
*> higher accuracy than the 'C' option. If the structure of the
|
| *> input matrix is not known, and relative accuracy is |
*> input matrix is not known, and relative accuracy is
|
| *> desirable, then this option is advisable. The input matrix A |
*> desirable, then this option is advisable. The input matrix A
|
| *> is preprocessed with QR factorization with FULL (row and |
*> is preprocessed with QR factorization with FULL (row and
|
| *> column) pivoting. |
*> column) pivoting.
|
| *> = 'G' Computation as with 'F' with an additional estimate of the |
*> = 'G' Computation as with 'F' with an additional estimate of the
|
| *> condition number of B, where A=D*B. If A has heavily weighted |
*> condition number of B, where A=B*D. If A has heavily weighted
|
| *> rows, then using this condition number gives too pessimistic |
*> rows, then using this condition number gives too pessimistic
|
| *> error bound. |
*> error bound.
|
| *> = 'A': Small singular values are the noise and the matrix is treated |
*> = 'A': Small singular values are not well determined by the data
|
| *> as numerically rank defficient. The error in the computed |
*> and are considered as noisy; the matrix is treated as
|
| *> singular values is bounded by f(m,n)*epsilon*||A||. |
*> numerically rank defficient. The error in the computed
|
| *> The computed SVD A = U * S * V^* restores A up to |
*> singular values is bounded by f(m,n)*epsilon*||A||.
|
| *> f(m,n)*epsilon*||A||. |
*> The computed SVD A = U * S * V^* restores A up to
|
| *> This gives the procedure the licence to discard (set to zero) |
*> f(m,n)*epsilon*||A||.
|
| *> all singular values below N*epsilon*||A||. |
*> This gives the procedure the licence to discard (set to zero)
|
| *> = 'R': Similar as in 'A'. Rank revealing property of the initial |
*> all singular values below N*epsilon*||A||.
|
| *> QR factorization is used do reveal (using triangular factor) |
*> = 'R': Similar as in 'A'. Rank revealing property of the initial
|
| *> a gap sigma_{r+1} < epsilon * sigma_r in which case the |
*> QR factorization is used do reveal (using triangular factor)
|
| *> numerical RANK is declared to be r. The SVD is computed with |
*> a gap sigma_{r+1} < epsilon * sigma_r in which case the
|
| *> absolute error bounds, but more accurately than with 'A'. |
*> numerical RANK is declared to be r. The SVD is computed with
|
| *> \endverbatim |
*> absolute error bounds, but more accurately than with 'A'.
|
| *> |
*> \endverbatim
|
| *> \param[in] JOBU |
*>
|
| *> \verbatim |
*> \param[in] JOBU
|
| *> JOBU is CHARACTER*1 |
*> \verbatim
|
| *> Specifies whether to compute the columns of U: |
*> JOBU is CHARACTER*1
|
| *> = 'U': N columns of U are returned in the array U. |
*> Specifies whether to compute the columns of U:
|
| *> = 'F': full set of M left sing. vectors is returned in the array U. |
*> = 'U': N columns of U are returned in the array U.
|
| *> = 'W': U may be used as workspace of length M*N. See the description |
*> = 'F': full set of M left sing. vectors is returned in the array U.
|
| *> of U. |
*> = 'W': U may be used as workspace of length M*N. See the description
|
| *> = 'N': U is not computed. |
*> of U.
|
| *> \endverbatim |
*> = 'N': U is not computed.
|
| *> |
*> \endverbatim
|
| *> \param[in] JOBV |
*>
|
| *> \verbatim |
*> \param[in] JOBV
|
| *> JOBV is CHARACTER*1 |
*> \verbatim
|
| *> Specifies whether to compute the matrix V: |
*> JOBV is CHARACTER*1
|
| *> = 'V': N columns of V are returned in the array V; Jacobi rotations |
*> Specifies whether to compute the matrix V:
|
| *> are not explicitly accumulated. |
*> = 'V': N columns of V are returned in the array V; Jacobi rotations
|
| *> = 'J': N columns of V are returned in the array V, but they are |
*> are not explicitly accumulated.
|
| *> computed as the product of Jacobi rotations. This option is |
*> = 'J': N columns of V are returned in the array V, but they are
|
| *> allowed only if JOBU .NE. 'N', i.e. in computing the full SVD. |
*> computed as the product of Jacobi rotations, if JOBT .EQ. 'N'.
|
| *> = 'W': V may be used as workspace of length N*N. See the description |
*> = 'W': V may be used as workspace of length N*N. See the description
|
| *> of V. |
*> of V.
|
| *> = 'N': V is not computed. |
*> = 'N': V is not computed.
|
| *> \endverbatim |
*> \endverbatim
|
| *> |
*>
|
| *> \param[in] JOBR |
*> \param[in] JOBR
|
| *> \verbatim |
*> \verbatim
|
| *> JOBR is CHARACTER*1 |
*> JOBR is CHARACTER*1
|
| *> Specifies the RANGE for the singular values. Issues the licence to |
*> Specifies the RANGE for the singular values. Issues the licence to
|
| *> set to zero small positive singular values if they are outside |
*> set to zero small positive singular values if they are outside
|
| *> specified range. If A .NE. 0 is scaled so that the largest singular |
*> 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 |
*> 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 |
*> the licence to kill columns of A whose norm in c*A is less than
|
| *> SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN, |
*> SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
|
| *> where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E'). |
*> where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E').
|
| *> = 'N': Do not kill small columns of c*A. This option assumes that |
*> = 'N': Do not kill small columns of c*A. This option assumes that
|
| *> BLAS and QR factorizations and triangular solvers are |
*> BLAS and QR factorizations and triangular solvers are
|
| *> implemented to work in that range. If the condition of A |
*> implemented to work in that range. If the condition of A
|
| *> is greater than BIG, use ZGESVJ. |
*> is greater than BIG, use ZGESVJ.
|
| *> = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)] |
*> = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)]
|
| *> (roughly, as described above). This option is recommended. |
*> (roughly, as described above). This option is recommended.
|
| *> =========================== |
*> ===========================
|
| *> For computing the singular values in the FULL range [SFMIN,BIG] |
*> For computing the singular values in the FULL range [SFMIN,BIG]
|
| *> use ZGESVJ. |
*> use ZGESVJ.
|
| *> \endverbatim |
*> \endverbatim
|
| *> |
*>
|
| *> \param[in] JOBT |
*> \param[in] JOBT
|
| *> \verbatim |
*> \verbatim
|
| *> JOBT is CHARACTER*1 |
*> JOBT is CHARACTER*1
|
| *> If the matrix is square then the procedure may determine to use |
*> If the matrix is square then the procedure may determine to use
|
| *> transposed A if A^* seems to be better with respect to convergence. |
*> transposed A if A^* seems to be better with respect to convergence.
|
| *> If the matrix is not square, JOBT is ignored. This is subject to |
*> If the matrix is not square, JOBT is ignored.
|
| *> changes in the future. |
*> The decision is based on two values of entropy over the adjoint
|
| *> 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).
|
| *> orbit of A^* * A. See the descriptions of WORK(6) and WORK(7). |
*> = 'T': transpose if entropy test indicates possibly faster
|
| *> = 'T': transpose if entropy test indicates possibly faster |
*> convergence of Jacobi process if A^* is taken as input. If A is
|
| *> convergence of Jacobi process if A^* is taken as input. If A is |
*> replaced with A^*, then the row pivoting is included automatically.
|
| *> replaced with A^*, then the row pivoting is included automatically. |
*> = 'N': do not speculate.
|
| *> = 'N': do not speculate. |
*> The option 'T' can be used to compute only the singular values, or
|
| *> This option can be used to compute only the singular values, or the |
*> the full SVD (U, SIGMA and V). For only one set of singular vectors
|
| *> 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
|
| *> (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.
|
| *> matrices is used as workspace if the matrix A is transposed. |
*> The implementer can easily remove this constraint and make the
|
| *> The implementer can easily remove this constraint and make the |
*> code more complicated. See the descriptions of U and V.
|
| *> code more complicated. See the descriptions of U and V. |
*> In general, this option is considered experimental, and 'N'; should
|
| *> \endverbatim |
*> be preferred. This is subject to changes in the future.
|
| *> |
*> \endverbatim
|
| *> \param[in] JOBP |
*>
|
| *> \verbatim |
*> \param[in] JOBP
|
| *> JOBP is CHARACTER*1 |
*> \verbatim
|
| *> Issues the licence to introduce structured perturbations to drown |
*> JOBP is CHARACTER*1
|
| *> denormalized numbers. This licence should be active if the |
*> Issues the licence to introduce structured perturbations to drown
|
| *> denormals are poorly implemented, causing slow computation, |
*> denormalized numbers. This licence should be active if the
|
| *> especially in cases of fast convergence (!). For details see [1,2]. |
*> denormals are poorly implemented, causing slow computation,
|
| *> For the sake of simplicity, this perturbations are included only |
*> especially in cases of fast convergence (!). For details see [1,2].
|
| *> when the full SVD or only the singular values are requested. The |
*> For the sake of simplicity, this perturbations are included only
|
| *> implementer/user can easily add the perturbation for the cases of |
*> when the full SVD or only the singular values are requested. The
|
| *> computing one set of singular vectors. |
*> implementer/user can easily add the perturbation for the cases of
|
| *> = 'P': introduce perturbation |
*> computing one set of singular vectors.
|
| *> = 'N': do not perturb |
*> = 'P': introduce perturbation
|
| *> \endverbatim |
*> = 'N': do not perturb
|
| *> |
*> \endverbatim
|
| *> \param[in] M |
*>
|
| *> \verbatim |
*> \param[in] M
|
| *> M is INTEGER |
*> \verbatim
|
| *> The number of rows of the input matrix A. M >= 0. |
*> M is INTEGER
|
| *> \endverbatim |
*> The number of rows of the input matrix A. M >= 0.
|
| *> |
*> \endverbatim
|
| *> \param[in] N |
*>
|
| *> \verbatim |
*> \param[in] N
|
| *> N is INTEGER |
*> \verbatim
|
| *> The number of columns of the input matrix A. M >= N >= 0. |
*> N is INTEGER
|
| *> \endverbatim |
*> The number of columns of the input matrix A. M >= N >= 0.
|
| *> |
*> \endverbatim
|
| *> \param[in,out] A |
*>
|
| *> \verbatim |
*> \param[in,out] A
|
| *> A is COMPLEX*16 array, dimension (LDA,N) |
*> \verbatim
|
| *> On entry, the M-by-N matrix A. |
*> A is COMPLEX*16 array, dimension (LDA,N)
|
| *> \endverbatim |
*> On entry, the M-by-N matrix A.
|
| *> |
*> \endverbatim
|
| *> \param[in] LDA |
*>
|
| *> \verbatim |
*> \param[in] LDA
|
| *> LDA is INTEGER |
*> \verbatim
|
| *> The leading dimension of the array A. LDA >= max(1,M). |
*> LDA is INTEGER
|
| *> \endverbatim |
*> The leading dimension of the array A. LDA >= max(1,M).
|
| *> |
*> \endverbatim
|
| *> \param[out] SVA |
*>
|
| *> \verbatim |
*> \param[out] SVA
|
| *> SVA is DOUBLE PRECISION array, dimension (N) |
*> \verbatim
|
| *> On exit, |
*> SVA is DOUBLE PRECISION array, dimension (N)
|
| *> - For WORK(1)/WORK(2) = ONE: The singular values of A. During the |
*> On exit,
|
| *> computation SVA contains Euclidean column norms of the |
*> - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
|
| *> iterated matrices in the array A. |
*> computation SVA contains Euclidean column norms of the
|
| *> - For WORK(1) .NE. WORK(2): The singular values of A are |
*> iterated matrices in the array A.
|
| *> (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if |
*> - For WORK(1) .NE. WORK(2): The singular values of A are
|
| *> sigma_max(A) overflows or if small singular values have been |
*> (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
|
| *> saved from underflow by scaling the input matrix A. |
*> sigma_max(A) overflows or if small singular values have been
|
| *> - If JOBR='R' then some of the singular values may be returned |
*> saved from underflow by scaling the input matrix A.
|
| *> as exact zeros obtained by "set to zero" because they are |
*> - If JOBR='R' then some of the singular values may be returned
|
| *> below the numerical rank threshold or are denormalized numbers. |
*> as exact zeros obtained by "set to zero" because they are
|
| *> \endverbatim |
*> below the numerical rank threshold or are denormalized numbers.
|
| *> |
*> \endverbatim
|
| *> \param[out] U |
*>
|
| *> \verbatim |
*> \param[out] U
|
| *> U is COMPLEX*16 array, dimension ( LDU, N ) |
*> \verbatim
|
| *> If JOBU = 'U', then U contains on exit the M-by-N matrix of |
*> U is COMPLEX*16 array, dimension ( LDU, N )
|
| *> the left singular vectors. |
*> If JOBU = 'U', then U contains on exit the M-by-N matrix of
|
| *> If JOBU = 'F', then U contains on exit the M-by-M matrix of |
*> the left singular vectors.
|
| *> the left singular vectors, including an ONB |
*> If JOBU = 'F', then U contains on exit the M-by-M matrix of
|
| *> of the orthogonal complement of the Range(A). |
*> the left singular vectors, including an ONB
|
| *> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N), |
*> of the orthogonal complement of the Range(A).
|
| *> then U is used as workspace if the procedure |
*> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
|
| *> replaces A with A^*. In that case, [V] is computed |
*> then U is used as workspace if the procedure
|
| *> in U as left singular vectors of A^* and then |
*> replaces A with A^*. In that case, [V] is computed
|
| *> copied back to the V array. This 'W' option is just |
*> in U as left singular vectors of A^* and then
|
| *> a reminder to the caller that in this case U is |
*> copied back to the V array. This 'W' option is just
|
| *> reserved as workspace of length N*N. |
*> a reminder to the caller that in this case U is
|
| *> If JOBU = 'N' U is not referenced, unless JOBT='T'. |
*> reserved as workspace of length N*N.
|
| *> \endverbatim |
*> If JOBU = 'N' U is not referenced, unless JOBT='T'.
|
| *> |
*> \endverbatim
|
| *> \param[in] LDU |
*>
|
| *> \verbatim |
*> \param[in] LDU
|
| *> LDU is INTEGER |
*> \verbatim
|
| *> The leading dimension of the array U, LDU >= 1. |
*> LDU is INTEGER
|
| *> IF JOBU = 'U' or 'F' or 'W', then LDU >= M. |
*> The leading dimension of the array U, LDU >= 1.
|
| *> \endverbatim |
*> IF JOBU = 'U' or 'F' or 'W', then LDU >= M.
|
| *> |
*> \endverbatim
|
| *> \param[out] V |
*>
|
| *> \verbatim |
*> \param[out] V
|
| *> V is COMPLEX*16 array, dimension ( LDV, N ) |
*> \verbatim
|
| *> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of |
*> V is COMPLEX*16 array, dimension ( LDV, N )
|
| *> the right singular vectors; |
*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
|
| *> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N), |
*> the right singular vectors;
|
| *> then V is used as workspace if the pprocedure |
*> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
|
| *> replaces A with A^*. In that case, [U] is computed |
*> then V is used as workspace if the pprocedure
|
| *> in V as right singular vectors of A^* and then |
*> replaces A with A^*. In that case, [U] is computed
|
| *> copied back to the U array. This 'W' option is just |
*> in V as right singular vectors of A^* and then
|
| *> a reminder to the caller that in this case V is |
*> copied back to the U array. This 'W' option is just
|
| *> reserved as workspace of length N*N. |
*> a reminder to the caller that in this case V is
|
| *> If JOBV = 'N' V is not referenced, unless JOBT='T'. |
*> reserved as workspace of length N*N.
|
| *> \endverbatim |
*> If JOBV = 'N' V is not referenced, unless JOBT='T'.
|
| *> |
*> \endverbatim
|
| *> \param[in] LDV |
*>
|
| *> \verbatim |
*> \param[in] LDV
|
| *> LDV is INTEGER |
*> \verbatim
|
| *> The leading dimension of the array V, LDV >= 1. |
*> LDV is INTEGER
|
| *> If JOBV = 'V' or 'J' or 'W', then LDV >= N. |
*> The leading dimension of the array V, LDV >= 1.
|
| *> \endverbatim |
*> If JOBV = 'V' or 'J' or 'W', then LDV >= N.
|
| *> |
*> \endverbatim
|
| *> \param[out] CWORK |
*>
|
| *> \verbatim |
*> \param[out] CWORK
|
| *> CWORK is COMPLEX*16 array, dimension at least LWORK. |
*> \verbatim
|
| *> \endverbatim |
*> CWORK is COMPLEX*16 array, dimension (MAX(2,LWORK))
|
| *> |
*> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
|
| *> \param[in] LWORK |
*> LRWORK=-1), then on exit CWORK(1) contains the required length of
|
| *> \verbatim |
*> CWORK for the job parameters used in the call.
|
| *> LWORK is INTEGER |
*> \endverbatim
|
| *> Length of CWORK to confirm proper allocation of workspace. |
*>
|
| *> LWORK depends on the job: |
*> \param[in] LWORK
|
| *> |
*> \verbatim
|
| *> 1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and |
*> LWORK is INTEGER
|
| *> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'): |
*> Length of CWORK to confirm proper allocation of workspace.
|
| *> LWORK >= 2*N+1. This is the minimal requirement. |
*> LWORK depends on the job:
|
| *> ->> For optimal performance (blocked code) the optimal value |
*>
|
| *> is LWORK >= N + (N+1)*NB. Here NB is the optimal |
*> 1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
|
| *> block size for ZGEQP3 and ZGEQRF. |
*> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):
|
| *> In general, optimal LWORK is computed as |
*> LWORK >= 2*N+1. This is the minimal requirement.
|
| *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF)). |
*> ->> For optimal performance (blocked code) the optimal value
|
| *> 1.2. .. an estimate of the scaled condition number of A is |
*> is LWORK >= N + (N+1)*NB. Here NB is the optimal
|
| *> required (JOBA='E', or 'G'). In this case, LWORK the minimal |
*> block size for ZGEQP3 and ZGEQRF.
|
| *> requirement is LWORK >= N*N + 3*N. |
*> In general, optimal LWORK is computed as
|
| *> ->> For optimal performance (blocked code) the optimal value |
*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ)).
|
| *> is LWORK >= max(N+(N+1)*NB, N*N+3*N). |
*> 1.2. .. an estimate of the scaled condition number of A is
|
| *> In general, the optimal length LWORK is computed as |
*> required (JOBA='E', or 'G'). In this case, LWORK the minimal
|
| *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), |
*> requirement is LWORK >= N*N + 2*N.
|
| *> N+N*N+LWORK(ZPOCON)). |
*> ->> For optimal performance (blocked code) the optimal value
|
| *> |
*> is LWORK >= max(N+(N+1)*NB, N*N+2*N)=N**2+2*N.
|
| *> 2. If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), |
*> In general, the optimal length LWORK is computed as
|
| *> (JOBU.EQ.'N') |
*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ),
|
| *> -> the minimal requirement is LWORK >= 3*N. |
*> N*N+LWORK(ZPOCON)).
|
| *> -> For optimal performance, LWORK >= max(N+(N+1)*NB, 3*N,2*N+N*NB), |
*> 2. If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
|
| *> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQF, |
*> (JOBU.EQ.'N')
|
| *> ZUNMLQ. In general, the optimal length LWORK is computed as |
*> 2.1 .. no scaled condition estimate requested (JOBE.EQ.'N'):
|
| *> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZPOCON), N+LWORK(ZGESVJ), |
*> -> the minimal requirement is LWORK >= 3*N.
|
| *> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)). |
*> -> For optimal performance,
|
| *> |
*> LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
|
| *> 3. If SIGMA and the left singular vectors are needed |
*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
|
| *> -> the minimal requirement is LWORK >= 3*N. |
*> ZUNMLQ. In general, the optimal length LWORK is computed as
|
| *> -> For optimal performance: |
*> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZGESVJ),
|
| *> if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB), |
*> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
|
| *> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR. |
*> 2.2 .. an estimate of the scaled condition number of A is
|
| *> In general, the optimal length LWORK is computed as |
*> required (JOBA='E', or 'G').
|
| *> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON), |
*> -> the minimal requirement is LWORK >= 3*N.
|
| *> 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). |
*> -> For optimal performance,
|
| *> |
*> LWORK >= max(N+(N+1)*NB, 2*N,2*N+N*NB)=2*N+N*NB,
|
| *> 4. If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and |
*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,
|
| *> 4.1. if JOBV.EQ.'V' |
*> ZUNMLQ. In general, the optimal length LWORK is computed as
|
| *> the minimal requirement is LWORK >= 5*N+2*N*N. |
*> LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ),
|
| *> 4.2. if JOBV.EQ.'J' the minimal requirement is |
*> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).
|
| *> LWORK >= 4*N+N*N. |
*> 3. If SIGMA and the left singular vectors are needed
|
| *> In both cases, the allocated CWORK can accommodate blocked runs |
*> 3.1 .. no scaled condition estimate requested (JOBE.EQ.'N'):
|
| *> of ZGEQP3, ZGEQRF, ZGELQF, ZUNMQR, ZUNMLQ. |
*> -> the minimal requirement is LWORK >= 3*N.
|
| *> \endverbatim |
*> -> For optimal performance:
|
| *> |
*> if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
|
| *> \param[out] RWORK |
*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
|
| *> \verbatim |
*> In general, the optimal length LWORK is computed as
|
| *> RWORK is DOUBLE PRECISION array, dimension at least LRWORK. |
*> LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
|
| *> On exit, |
*> 3.2 .. an estimate of the scaled condition number of A is
|
| *> RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1) |
*> required (JOBA='E', or 'G').
|
| *> such that SCALE*SVA(1:N) are the computed singular values |
*> -> the minimal requirement is LWORK >= 3*N.
|
| *> of A. (See the description of SVA().) |
*> -> For optimal performance:
|
| *> RWORK(2) = See the description of RWORK(1). |
*> if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
|
| *> RWORK(3) = SCONDA is an estimate for the condition number of |
*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
|
| *> column equilibrated A. (If JOBA .EQ. 'E' or 'G') |
*> In general, the optimal length LWORK is computed as
|
| *> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). |
*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON),
|
| *> It is computed using SPOCON. It holds |
*> 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
|
| *> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA |
*> 4. If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and
|
| *> where R is the triangular factor from the QRF of A. |
*> 4.1. if JOBV.EQ.'V'
|
| *> However, if R is truncated and the numerical rank is |
*> the minimal requirement is LWORK >= 5*N+2*N*N.
|
| *> determined to be strictly smaller than N, SCONDA is |
*> 4.2. if JOBV.EQ.'J' the minimal requirement is
|
| *> returned as -1, thus indicating that the smallest |
*> LWORK >= 4*N+N*N.
|
| *> singular values might be lost. |
*> In both cases, the allocated CWORK can accomodate blocked runs
|
| *> |
*> of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ.
|
| *> If full SVD is needed, the following two condition numbers are |
*>
|
| *> useful for the analysis of the algorithm. They are provied for |
*> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
|
| *> a developer/implementer who is familiar with the details of |
*> LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the
|
| *> the method. |
*> minimal length of CWORK for the job parameters used in the call.
|
| *> |
*> \endverbatim
|
| *> RWORK(4) = an estimate of the scaled condition number of the |
*>
|
| *> triangular factor in the first QR factorization. |
*> \param[out] RWORK
|
| *> RWORK(5) = an estimate of the scaled condition number of the |
*> \verbatim
|
| *> triangular factor in the second QR factorization. |
*> RWORK is DOUBLE PRECISION array, dimension (MAX(7,LWORK))
|
| *> The following two parameters are computed if JOBT .EQ. 'T'. |
*> On exit,
|
| *> They are provided for a developer/implementer who is familiar |
*> RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1)
|
| *> with the details of the method. |
*> such that SCALE*SVA(1:N) are the computed singular values
|
| *> RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy |
*> of A. (See the description of SVA().)
|
| *> of diag(A^* * A) / Trace(A^* * A) taken as point in the |
*> RWORK(2) = See the description of RWORK(1).
|
| *> probability simplex. |
*> RWORK(3) = SCONDA is an estimate for the condition number of
|
| *> RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).) |
*> column equilibrated A. (If JOBA .EQ. 'E' or 'G')
|
| *> \endverbatim |
*> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
|
| *> |
*> It is computed using SPOCON. It holds
|
| *> \param[in] LRWORK |
*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
|
| *> \verbatim |
*> where R is the triangular factor from the QRF of A.
|
| *> LRWORK is INTEGER |
*> However, if R is truncated and the numerical rank is
|
| *> Length of RWORK to confirm proper allocation of workspace. |
*> determined to be strictly smaller than N, SCONDA is
|
| *> LRWORK depends on the job: |
*> returned as -1, thus indicating that the smallest
|
| *> |
*> singular values might be lost.
|
| *> 1. If only singular values are requested i.e. if |
*>
|
| *> LSAME(JOBU,'N') .AND. LSAME(JOBV,'N') |
*> If full SVD is needed, the following two condition numbers are
|
| *> then: |
*> useful for the analysis of the algorithm. They are provied for
|
| *> 1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), |
*> a developer/implementer who is familiar with the details of
|
| *> then LRWORK = max( 7, N + 2 * M ). |
*> the method.
|
| *> 1.2. Otherwise, LRWORK = max( 7, 2 * N ). |
*>
|
| *> 2. If singular values with the right singular vectors are requested |
*> RWORK(4) = an estimate of the scaled condition number of the
|
| *> i.e. if |
*> triangular factor in the first QR factorization.
|
| *> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND. |
*> RWORK(5) = an estimate of the scaled condition number of the
|
| *> .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) |
*> triangular factor in the second QR factorization.
|
| *> then: |
*> The following two parameters are computed if JOBT .EQ. 'T'.
|
| *> 2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), |
*> They are provided for a developer/implementer who is familiar
|
| *> then LRWORK = max( 7, N + 2 * M ). |
*> with the details of the method.
|
| *> 2.2. Otherwise, LRWORK = max( 7, 2 * N ). |
*> RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy
|
| *> 3. If singular values with the left singular vectors are requested, i.e. if |
*> of diag(A^* * A) / Trace(A^* * A) taken as point in the
|
| *> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. |
*> probability simplex.
|
| *> .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) |
*> RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).)
|
| *> then: |
*> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or
|
| *> 3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), |
*> LRWORK=-1), then on exit RWORK(1) contains the required length of
|
| *> then LRWORK = max( 7, N + 2 * M ). |
*> RWORK for the job parameters used in the call.
|
| *> 3.2. Otherwise, LRWORK = max( 7, 2 * N ). |
*> \endverbatim
|
| *> 4. If singular values with both the left and the right singular vectors |
*>
|
| *> are requested, i.e. if |
*> \param[in] LRWORK
|
| *> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. |
*> \verbatim
|
| *> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) |
*> LRWORK is INTEGER
|
| *> then: |
*> Length of RWORK to confirm proper allocation of workspace.
|
| *> 4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), |
*> LRWORK depends on the job:
|
| *> then LRWORK = max( 7, N + 2 * M ). |
*>
|
| *> 4.2. Otherwise, LRWORK = max( 7, 2 * N ). |
*> 1. If only the singular values are requested i.e. if
|
| *> \endverbatim |
*> LSAME(JOBU,'N') .AND. LSAME(JOBV,'N')
|
| *> |
*> then:
|
| *> \param[out] IWORK |
*> 1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
|
| *> \verbatim |
*> then: LRWORK = max( 7, 2 * M ).
|
| *> IWORK is INTEGER array, of dimension: |
*> 1.2. Otherwise, LRWORK = max( 7, N ).
|
| *> If LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), then |
*> 2. If singular values with the right singular vectors are requested
|
| *> the dimension of IWORK is max( 3, 2 * N + M ). |
*> i.e. if
|
| *> Otherwise, the dimension of IWORK is |
*> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND.
|
| *> -> max( 3, 2*N ) for full SVD |
*> .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F'))
|
| *> -> max( 3, N ) for singular values only or singular |
*> then:
|
| *> values with one set of singular vectors (left or right) |
*> 2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
|
| *> On exit, |
*> then LRWORK = max( 7, 2 * M ).
|
| *> IWORK(1) = the numerical rank determined after the initial |
*> 2.2. Otherwise, LRWORK = max( 7, N ).
|
| *> QR factorization with pivoting. See the descriptions |
*> 3. If singular values with the left singular vectors are requested, i.e. if
|
| *> of JOBA and JOBR. |
*> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND.
|
| *> IWORK(2) = the number of the computed nonzero singular values |
*> .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J'))
|
| *> IWORK(3) = if nonzero, a warning message: |
*> then:
|
| *> If IWORK(3).EQ.1 then some of the column norms of A |
*> 3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
|
| *> were denormalized floats. The requested high accuracy |
*> then LRWORK = max( 7, 2 * M ).
|
| *> is not warranted by the data. |
*> 3.2. Otherwise, LRWORK = max( 7, N ).
|
| *> \endverbatim |
*> 4. If singular values with both the left and the right singular vectors
|
| *> |
*> are requested, i.e. if
|
| *> \param[out] INFO |
*> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND.
|
| *> \verbatim |
*> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J'))
|
| *> INFO is INTEGER |
*> then:
|
| *> < 0 : if INFO = -i, then the i-th argument had an illegal value. |
*> 4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'),
|
| *> = 0 : successfull exit; |
*> then LRWORK = max( 7, 2 * M ).
|
| *> > 0 : ZGEJSV did not converge in the maximal allowed number |
*> 4.2. Otherwise, LRWORK = max( 7, N ).
|
| *> of sweeps. The computed values may be inaccurate. |
*>
|
| *> \endverbatim |
*> If, on entry, LRWORK = -1 ot LWORK=-1, a workspace query is assumed and
|
| * |
*> the length of RWORK is returned in RWORK(1).
|
| * Authors: |
*> \endverbatim
|
| * ======== |
*>
|
| * |
*> \param[out] IWORK
|
| *> \author Univ. of Tennessee |
*> \verbatim
|
| *> \author Univ. of California Berkeley |
*> IWORK is INTEGER array, of dimension at least 4, that further depends
|
| *> \author Univ. of Colorado Denver |
*> on the job:
|
| *> \author NAG Ltd. |
*>
|
| * |
*> 1. If only the singular values are requested then:
|
| *> \date June 2016 |
*> 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.
|
| *> \ingroup complex16GEsing |
*> 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') )
|
| *> \par Further Details: |
*> 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') )
|
| *> \verbatim |
*> 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
|
| *> ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3, |
*> are requested, then:
|
| *> ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an |
*> 4.1. If LSAME(JOBV,'J') the length of IWORK is determined as follows:
|
| *> additional row pivoting can be used as a preprocessor, which in some |
*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
|
| *> cases results in much higher accuracy. An example is matrix A with the |
*> then the length of IWORK is N+M; otherwise the length of IWORK is N.
|
| *> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned |
*> 4.2. If LSAME(JOBV,'V') the length of IWORK is determined as follows:
|
| *> diagonal matrices and C is well-conditioned matrix. In that case, complete |
*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') )
|
| *> pivoting in the first QR factorizations provides accuracy dependent on the |
*> then the length of IWORK is 2*N+M; otherwise the length of IWORK is 2*N.
|
| *> condition number of C, and independent of D1, D2. Such higher accuracy is |
*>
|
| *> not completely understood theoretically, but it works well in practice. |
*> On exit,
|
| *> Further, if A can be written as A = B*D, with well-conditioned B and some |
*> IWORK(1) = the numerical rank determined after the initial
|
| *> diagonal D, then the high accuracy is guaranteed, both theoretically and |
*> QR factorization with pivoting. See the descriptions
|
| *> in software, independent of D. For more details see [1], [2]. |
*> of JOBA and JOBR.
|
| *> The computational range for the singular values can be the full range |
*> IWORK(2) = the number of the computed nonzero singular values
|
| *> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS |
*> IWORK(3) = if nonzero, a warning message:
|
| *> & LAPACK routines called by ZGEJSV are implemented to work in that range. |
*> If IWORK(3).EQ.1 then some of the column norms of A
|
| *> If that is not the case, then the restriction for safe computation with |
*> were denormalized floats. The requested high accuracy
|
| *> the singular values in the range of normalized IEEE numbers is that the |
*> is not warranted by the data.
|
| *> spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not |
*> IWORK(4) = 1 or -1. If IWORK(4) .EQ. 1, then the procedure used A^* to
|
| *> overflow. This code (ZGEJSV) is best used in this restricted range, |
*> do the job as specified by the JOB parameters.
|
| *> meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are |
*> If the call to ZGEJSV is a workspace query (indicated by LWORK .EQ. -1 or
|
| *> returned as zeros. See JOBR for details on this. |
*> LRWORK .EQ. -1), then on exit IWORK(1) contains the required length of
|
| *> Further, this implementation is somewhat slower than the one described |
*> IWORK for the job parameters used in the call.
|
| *> in [1,2] due to replacement of some non-LAPACK components, and because |
*> \endverbatim
|
| *> the choice of some tuning parameters in the iterative part (ZGESVJ) is |
*>
|
| *> left to the implementer on a particular machine. |
*> \param[out] INFO
|
| *> The rank revealing QR factorization (in this code: ZGEQP3) should be |
*> \verbatim
|
| *> implemented as in [3]. We have a new version of ZGEQP3 under development |
*> INFO is INTEGER
|
| *> that is more robust than the current one in LAPACK, with a cleaner cut in |
*> < 0 : if INFO = -i, then the i-th argument had an illegal value.
|
| *> rank defficient cases. It will be available in the SIGMA library [4]. |
*> = 0 : successful exit;
|
| *> If M is much larger than N, it is obvious that the inital QRF with |
*> > 0 : ZGEJSV did not converge in the maximal allowed number
|
| *> column pivoting can be preprocessed by the QRF without pivoting. That |
*> of sweeps. The computed values may be inaccurate.
|
| *> well known trick is not used in ZGEJSV because in some cases heavy row |
*> \endverbatim
|
| *> 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 |
* Authors:
|
| *> 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 |
*> \author Univ. of Tennessee
|
| *> implementation of ZGEJSV uses only the simplest, naive data movement. |
*> \author Univ. of California Berkeley
|
| * |
*> \author Univ. of Colorado Denver
|
| *> \par Contributors: |
*> \author NAG Ltd.
|
| * ================== |
*
|
| *> |
*> \date June 2016
|
| *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) |
*
|
| * |
*> \ingroup complex16GEsing
|
| *> \par References: |
*
|
| * ================ |
*> \par Further Details:
|
| *> |
* =====================
|
| *> \verbatim |
*>
|
| *> |
*> \verbatim
|
| *> [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. |
*> ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3,
|
| *> LAPACK Working note 169. |
*> ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an
|
| *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. |
*> additional row pivoting can be used as a preprocessor, which in some
|
| *> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. |
*> cases results in much higher accuracy. An example is matrix A with the
|
| *> LAPACK Working note 170. |
*> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
|
| *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR |
*> diagonal matrices and C is well-conditioned matrix. In that case, complete
|
| *> factorization software - a case study. |
*> pivoting in the first QR factorizations provides accuracy dependent on the
|
| *> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. |
*> condition number of C, and independent of D1, D2. Such higher accuracy is
|
| *> LAPACK Working note 176. |
*> not completely understood theoretically, but it works well in practice.
|
| *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, |
*> Further, if A can be written as A = B*D, with well-conditioned B and some
|
| *> QSVD, (H,K)-SVD computations. |
*> diagonal D, then the high accuracy is guaranteed, both theoretically and
|
| *> Department of Mathematics, University of Zagreb, 2008. |
*> in software, independent of D. For more details see [1], [2].
|
| *> \endverbatim |
*> The computational range for the singular values can be the full range
|
| * |
*> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
|
| *> \par Bugs, examples and comments: |
*> & 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
|
| *> Please report all bugs and send interesting examples and/or comments to |
*> spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
|
| *> drmac@math.hr. Thank you. |
*> 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.
|
| SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, |
*> Further, this implementation is somewhat slower than the one described
|
| $ M, N, A, LDA, SVA, U, LDU, V, LDV, |
*> in [1,2] due to replacement of some non-LAPACK components, and because
|
| $ CWORK, LWORK, RWORK, LRWORK, IWORK, INFO ) |
*> the choice of some tuning parameters in the iterative part (ZGESVJ) is
|
| * |
*> left to the implementer on a particular machine.
|
| * -- LAPACK computational routine (version 3.6.1) -- |
*> The rank revealing QR factorization (in this code: ZGEQP3) should be
|
| * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
*> implemented as in [3]. We have a new version of ZGEQP3 under development
|
| * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
*> that is more robust than the current one in LAPACK, with a cleaner cut in
|
| * June 2016 |
*> 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
|
| * .. Scalar Arguments .. |
*> column pivoting can be preprocessed by the QRF without pivoting. That
|
| IMPLICIT NONE |
*> well known trick is not used in ZGEJSV because in some cases heavy row
|
| INTEGER INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N |
*> 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
|
| * .. Array Arguments .. |
*> terms of accuracy have prevailed. The implementer/user can incorporate
|
| COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), |
*> this extra QRF step easily. The implementer can also improve data movement
|
| $ CWORK( LWORK ) |
*> (matrix transpose, matrix copy, matrix transposed copy) - this
|
| DOUBLE PRECISION SVA( N ), RWORK( * ) |
*> implementation of ZGEJSV uses only the simplest, naive data movement.
|
| INTEGER IWORK( * ) |
*> \endverbatim
|
| CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV |
*
|
| * .. |
*> \par Contributor:
|
| * |
* ==================
|
| * =========================================================================== |
*>
|
| * |
*> Zlatko Drmac, Department of Mathematics, Faculty of Science,
|
| * .. Local Parameters .. |
*> University of Zagreb (Zagreb, Croatia); drmac@math.hr
|
| DOUBLE PRECISION ZERO, ONE |
*
|
| PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) |
*> \par References:
|
| COMPLEX*16 CZERO, CONE |
* ================
|
| PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) ) |
*>
|
| * .. |
*> \verbatim
|
| * .. Local Scalars .. |
*>
|
| COMPLEX*16 CTEMP |
*> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
|
| DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, |
*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
|
| $ COND_OK, CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, |
*> LAPACK Working note 169.
|
| $ MAXPRJ, SCALEM, SCONDA, SFMIN, SMALL, TEMP1, |
*> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
|
| $ USCAL1, USCAL2, XSC |
*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
|
| INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING |
*> LAPACK Working note 170.
|
| LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC, |
*> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
|
| $ L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, |
*> factorization software - a case study.
|
| $ NOSCAL, ROWPIV, RSVEC, TRANSP |
*> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
|
| * .. |
*> LAPACK Working note 176.
|
| * .. Intrinsic Functions .. |
*> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
|
| INTRINSIC ABS, DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DBLE, |
*> QSVD, (H,K)-SVD computations.
|
| $ MAX0, MIN0, NINT, DSQRT |
*> Department of Mathematics, University of Zagreb, 2008, 2016.
|
| * .. |
*> \endverbatim
|
| * .. External Functions .. |
*
|
| DOUBLE PRECISION DLAMCH, DZNRM2 |
*> \par Bugs, examples and comments:
|
| INTEGER IDAMAX, IZAMAX |
* =================================
|
| LOGICAL LSAME |
*>
|
| EXTERNAL IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2 |
*> Please report all bugs and send interesting examples and/or comments to
|
| * .. |
*> drmac@math.hr. Thank you.
|
| * .. External Subroutines .. |
*>
|
| EXTERNAL DLASSQ, ZCOPY, ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL, |
* =====================================================================
|
| $ DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ, |
SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
|
| $ ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP, ZTRSM, XERBLA |
$ M, N, A, LDA, SVA, U, LDU, V, LDV,
|
| * |
$ CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
|
| EXTERNAL ZGESVJ |
*
|
| * .. |
* -- LAPACK computational routine (version 3.7.0) --
|
| * |
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
| * Test the input arguments |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
| * |
* December 2016
|
| |
*
|
| LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' ) |
* .. Scalar Arguments ..
|
| JRACC = LSAME( JOBV, 'J' ) |
IMPLICIT NONE
|
| RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC |
INTEGER INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N
|
| ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' ) |
* ..
|
| L2RANK = LSAME( JOBA, 'R' ) |
* .. Array Arguments ..
|
| L2ABER = LSAME( JOBA, 'A' ) |
COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ),
|
| ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' ) |
$ CWORK( LWORK )
|
| L2TRAN = LSAME( JOBT, 'T' ) |
DOUBLE PRECISION SVA( N ), RWORK( LRWORK )
|
| L2KILL = LSAME( JOBR, 'R' ) |
INTEGER IWORK( * )
|
| DEFR = LSAME( JOBR, 'N' ) |
CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
|
| L2PERT = LSAME( JOBP, 'P' ) |
* ..
|
| * |
*
|
| IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR. |
* ===========================================================================
|
| $ ERREST .OR. LSAME( JOBA, 'C' ) )) THEN |
*
|
| INFO = - 1 |
* .. Local Parameters ..
|
| ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR. |
DOUBLE PRECISION ZERO, ONE
|
| $ LSAME( JOBU, 'W' )) ) THEN |
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
|
| INFO = - 2 |
COMPLEX*16 CZERO, CONE
|
| ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR. |
PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) )
|
| $ LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN |
* ..
|
| INFO = - 3 |
* .. Local Scalars ..
|
| ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN |
COMPLEX*16 CTEMP
|
| INFO = - 4 |
DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1,
|
| ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN |
$ COND_OK, CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN,
|
| INFO = - 5 |
$ MAXPRJ, SCALEM, SCONDA, SFMIN, SMALL, TEMP1,
|
| ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN |
$ USCAL1, USCAL2, XSC
|
| INFO = - 6 |
INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
|
| ELSE IF ( M .LT. 0 ) THEN |
LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LQUERY,
|
| INFO = - 7 |
$ LSVEC, L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, NOSCAL,
|
| ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN |
$ ROWPIV, RSVEC, TRANSP
|
| INFO = - 8 |
*
|
| ELSE IF ( LDA .LT. M ) THEN |
INTEGER OPTWRK, MINWRK, MINRWRK, MINIWRK
|
| INFO = - 10 |
INTEGER LWCON, LWLQF, LWQP3, LWQRF, LWUNMLQ, LWUNMQR, LWUNMQRM,
|
| ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN |
$ LWSVDJ, LWSVDJV, LRWQP3, LRWCON, LRWSVDJ, IWOFF
|
| INFO = - 13 |
INTEGER LWRK_ZGELQF, LWRK_ZGEQP3, LWRK_ZGEQP3N, LWRK_ZGEQRF,
|
| ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN |
$ LWRK_ZGESVJ, LWRK_ZGESVJV, LWRK_ZGESVJU, LWRK_ZUNMLQ,
|
| INFO = - 15 |
$ LWRK_ZUNMQR, LWRK_ZUNMQRM
|
| ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. |
* ..
|
| $ (LWORK .LT. 2*N+1)) .OR. |
* .. Local Arrays
|
| $ (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND. |
COMPLEX*16 CDUMMY(1)
|
| $ (LWORK .LT. N*N+3*N)) .OR. |
DOUBLE PRECISION RDUMMY(1)
|
| $ (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. 3*N)) |
*
|
| $ .OR. |
* .. Intrinsic Functions ..
|
| $ (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. 3*N)) |
INTRINSIC ABS, DCMPLX, CONJG, DLOG, MAX, MIN, DBLE, NINT, SQRT
|
| $ .OR. |
* ..
|
| $ (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. |
* .. External Functions ..
|
| $ (LWORK.LT.5*N+2*N*N)) |
DOUBLE PRECISION DLAMCH, DZNRM2
|
| $ .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND. |
INTEGER IDAMAX, IZAMAX
|
| $ LWORK.LT.4*N+N*N)) |
LOGICAL LSAME
|
| $ THEN |
EXTERNAL IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2
|
| INFO = - 17 |
* ..
|
| ELSE IF ( LRWORK.LT. MAX0(N+2*M,7)) THEN |
* .. External Subroutines ..
|
| INFO = -19 |
EXTERNAL DLASSQ, ZCOPY, ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLAPMR,
|
| ELSE |
$ ZLASCL, DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
|
| * #:) |
$ ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP, ZTRSM, ZLACGV,
|
| INFO = 0 |
$ XERBLA
|
| END IF |
*
|
| * |
EXTERNAL ZGESVJ
|
| IF ( INFO .NE. 0 ) THEN |
* ..
|
| * #:( |
*
|
| CALL XERBLA( 'ZGEJSV', - INFO ) |
* Test the input arguments
|
| RETURN |
*
|
| END IF |
LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
|
| * |
JRACC = LSAME( JOBV, 'J' )
|
| * Quick return for void matrix (Y3K safe) |
RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC
|
| * #:) |
ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
|
| IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN |
L2RANK = LSAME( JOBA, 'R' )
|
| IWORK(1:3) = 0 |
L2ABER = LSAME( JOBA, 'A' )
|
| RWORK(1:7) = 0 |
ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
|
| RETURN |
L2TRAN = LSAME( JOBT, 'T' ) .AND. ( M .EQ. N )
|
| ENDIF |
L2KILL = LSAME( JOBR, 'R' )
|
| * |
DEFR = LSAME( JOBR, 'N' )
|
| * Determine whether the matrix U should be M x N or M x M |
L2PERT = LSAME( JOBP, 'P' )
|
| * |
*
|
| IF ( LSVEC ) THEN |
LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
|
| N1 = N |
*
|
| IF ( LSAME( JOBU, 'F' ) ) N1 = M |
IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
|
| END IF |
$ ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
|
| * |
INFO = - 1
|
| * Set numerical parameters |
ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR.
|
| * |
$ ( LSAME( JOBU, 'W' ) .AND. RSVEC .AND. L2TRAN ) ) ) THEN
|
| *! NOTE: Make sure DLAMCH() does not fail on the target architecture. |
INFO = - 2
|
| * |
ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
|
| EPSLN = DLAMCH('Epsilon') |
$ ( LSAME( JOBV, 'W' ) .AND. LSVEC .AND. L2TRAN ) ) ) THEN
|
| SFMIN = DLAMCH('SafeMinimum') |
INFO = - 3
|
| SMALL = SFMIN / EPSLN |
ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN
|
| BIG = DLAMCH('O') |
INFO = - 4
|
| * BIG = ONE / SFMIN |
ELSE IF ( .NOT. ( LSAME(JOBT,'T') .OR. LSAME(JOBT,'N') ) ) THEN
|
| * |
INFO = - 5
|
| * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N |
ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
|
| * |
INFO = - 6
|
| *(!) If necessary, scale SVA() to protect the largest norm from |
ELSE IF ( M .LT. 0 ) THEN
|
| * overflow. It is possible that this scaling pushes the smallest |
INFO = - 7
|
| * column norm left from the underflow threshold (extreme case). |
ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
|
| * |
INFO = - 8
|
| SCALEM = ONE / DSQRT(DBLE(M)*DBLE(N)) |
ELSE IF ( LDA .LT. M ) THEN
|
| NOSCAL = .TRUE. |
INFO = - 10
|
| GOSCAL = .TRUE. |
ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
|
| DO 1874 p = 1, N |
INFO = - 13
|
| AAPP = ZERO |
ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
|
| AAQQ = ONE |
INFO = - 15
|
| CALL ZLASSQ( M, A(1,p), 1, AAPP, AAQQ ) |
ELSE
|
| IF ( AAPP .GT. BIG ) THEN |
* #:)
|
| INFO = - 9 |
INFO = 0
|
| CALL XERBLA( 'ZGEJSV', -INFO ) |
END IF
|
| RETURN |
*
|
| END IF |
IF ( INFO .EQ. 0 ) THEN
|
| AAQQ = DSQRT(AAQQ) |
* .. compute the minimal and the optimal workspace lengths
|
| IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN |
* [[The expressions for computing the minimal and the optimal
|
| SVA(p) = AAPP * AAQQ |
* values of LCWORK, LRWORK are written with a lot of redundancy and
|
| ELSE |
* can be simplified. However, this verbose form is useful for
|
| NOSCAL = .FALSE. |
* maintenance and modifications of the code.]]
|
| SVA(p) = AAPP * ( AAQQ * SCALEM ) |
*
|
| IF ( GOSCAL ) THEN |
* .. minimal workspace length for ZGEQP3 of an M x N matrix,
|
| GOSCAL = .FALSE. |
* ZGEQRF of an N x N matrix, ZGELQF of an N x N matrix,
|
| CALL DSCAL( p-1, SCALEM, SVA, 1 ) |
* ZUNMLQ for computing N x N matrix, ZUNMQR for computing N x N
|
| END IF |
* matrix, ZUNMQR for computing M x N matrix, respectively.
|
| END IF |
LWQP3 = N+1
|
| 1874 CONTINUE |
LWQRF = MAX( 1, N )
|
| * |
LWLQF = MAX( 1, N )
|
| IF ( NOSCAL ) SCALEM = ONE |
LWUNMLQ = MAX( 1, N )
|
| * |
LWUNMQR = MAX( 1, N )
|
| AAPP = ZERO |
LWUNMQRM = MAX( 1, M )
|
| AAQQ = BIG |
* .. minimal workspace length for ZPOCON of an N x N matrix
|
| DO 4781 p = 1, N |
LWCON = 2 * N
|
| AAPP = DMAX1( AAPP, SVA(p) ) |
* .. minimal workspace length for ZGESVJ of an N x N matrix,
|
| IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) ) |
* without and with explicit accumulation of Jacobi rotations
|
| 4781 CONTINUE |
LWSVDJ = MAX( 2 * N, 1 )
|
| * |
LWSVDJV = MAX( 2 * N, 1 )
|
| * Quick return for zero M x N matrix |
* .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ
|
| * #:) |
LRWQP3 = N
|
| IF ( AAPP .EQ. ZERO ) THEN |
LRWCON = N
|
| IF ( LSVEC ) CALL ZLASET( 'G', M, N1, CZERO, CONE, U, LDU ) |
LRWSVDJ = N
|
| IF ( RSVEC ) CALL ZLASET( 'G', N, N, CZERO, CONE, V, LDV ) |
IF ( LQUERY ) THEN
|
| RWORK(1) = ONE |
CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,
|
| RWORK(2) = ONE |
$ RDUMMY, IERR )
|
| IF ( ERREST ) RWORK(3) = ONE |
LWRK_ZGEQP3 = CDUMMY(1)
|
| IF ( LSVEC .AND. RSVEC ) THEN |
CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
|
| RWORK(4) = ONE |
LWRK_ZGEQRF = CDUMMY(1)
|
| RWORK(5) = ONE |
CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
|
| END IF |
LWRK_ZGELQF = CDUMMY(1)
|
| IF ( L2TRAN ) THEN |
END IF
|
| RWORK(6) = ZERO |
MINWRK = 2
|
| RWORK(7) = ZERO |
OPTWRK = 2
|
| END IF |
MINIWRK = N
|
| IWORK(1) = 0 |
IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN
|
| IWORK(2) = 0 |
* .. minimal and optimal sizes of the complex workspace if
|
| IWORK(3) = 0 |
* only the singular values are requested
|
| RETURN |
IF ( ERREST ) THEN
|
| END IF |
MINWRK = MAX( N+LWQP3, N**2+LWCON, N+LWQRF, LWSVDJ )
|
| * |
ELSE
|
| * Issue warning if denormalized column norms detected. Override the |
MINWRK = MAX( N+LWQP3, N+LWQRF, LWSVDJ )
|
| * high relative accuracy request. Issue licence to kill columns |
END IF
|
| * (set them to zero) whose norm is less than sigma_max / BIG (roughly). |
IF ( LQUERY ) THEN
|
| * #:( |
CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V,
|
| WARNING = 0 |
$ LDV, CDUMMY, -1, RDUMMY, -1, IERR )
|
| IF ( AAQQ .LE. SFMIN ) THEN |
LWRK_ZGESVJ = CDUMMY(1)
|
| L2RANK = .TRUE. |
IF ( ERREST ) THEN
|
| L2KILL = .TRUE. |
OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON,
|
| WARNING = 1 |
$ N+LWRK_ZGEQRF, LWRK_ZGESVJ )
|
| END IF |
ELSE
|
| * |
OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWRK_ZGEQRF,
|
| * Quick return for one-column matrix |
$ LWRK_ZGESVJ )
|
| * #:) |
END IF
|
| IF ( N .EQ. 1 ) THEN |
END IF
|
| * |
IF ( L2TRAN .OR. ROWPIV ) THEN
|
| IF ( LSVEC ) THEN |
IF ( ERREST ) THEN
|
| CALL ZLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR ) |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWCON, LRWSVDJ )
|
| CALL ZLACPY( 'A', M, 1, A, LDA, U, LDU ) |
ELSE
|
| * computing all M left singular vectors of the M x 1 matrix |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ )
|
| IF ( N1 .NE. N ) THEN |
END IF
|
| CALL ZGEQRF( M, N, U,LDU, CWORK, CWORK(N+1),LWORK-N,IERR ) |
ELSE
|
| CALL ZUNGQR( M,N1,1, U,LDU,CWORK,CWORK(N+1),LWORK-N,IERR ) |
IF ( ERREST ) THEN
|
| CALL ZCOPY( M, A(1,1), 1, U(1,1), 1 ) |
MINRWRK = MAX( 7, LRWQP3, LRWCON, LRWSVDJ )
|
| END IF |
ELSE
|
| END IF |
MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
|
| IF ( RSVEC ) THEN |
END IF
|
| V(1,1) = CONE |
END IF
|
| END IF |
IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
|
| IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN |
ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
|
| SVA(1) = SVA(1) / SCALEM |
* .. minimal and optimal sizes of the complex workspace if the
|
| SCALEM = ONE |
* singular values and the right singular vectors are requested
|
| END IF |
IF ( ERREST ) THEN
|
| RWORK(1) = ONE / SCALEM |
MINWRK = MAX( N+LWQP3, LWCON, LWSVDJ, N+LWLQF,
|
| RWORK(2) = ONE |
$ 2*N+LWQRF, N+LWSVDJ, N+LWUNMLQ )
|
| IF ( SVA(1) .NE. ZERO ) THEN |
ELSE
|
| IWORK(1) = 1 |
MINWRK = MAX( N+LWQP3, LWSVDJ, N+LWLQF, 2*N+LWQRF,
|
| IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN |
$ N+LWSVDJ, N+LWUNMLQ )
|
| IWORK(2) = 1 |
END IF
|
| ELSE |
IF ( LQUERY ) THEN
|
| IWORK(2) = 0 |
CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
|
| END IF |
$ LDA, CDUMMY, -1, RDUMMY, -1, IERR )
|
| ELSE |
LWRK_ZGESVJ = CDUMMY(1)
|
| IWORK(1) = 0 |
CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
|
| IWORK(2) = 0 |
$ V, LDV, CDUMMY, -1, IERR )
|
| END IF |
LWRK_ZUNMLQ = CDUMMY(1)
|
| IWORK(3) = 0 |
IF ( ERREST ) THEN
|
| IF ( ERREST ) RWORK(3) = ONE |
OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ,
|
| IF ( LSVEC .AND. RSVEC ) THEN |
$ N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF,
|
| RWORK(4) = ONE |
$ N+LWRK_ZGESVJ, N+LWRK_ZUNMLQ )
|
| RWORK(5) = ONE |
ELSE
|
| END IF |
OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVJ,N+LWRK_ZGELQF,
|
| IF ( L2TRAN ) THEN |
$ 2*N+LWRK_ZGEQRF, N+LWRK_ZGESVJ,
|
| RWORK(6) = ZERO |
$ N+LWRK_ZUNMLQ )
|
| RWORK(7) = ZERO |
END IF
|
| END IF |
END IF
|
| RETURN |
IF ( L2TRAN .OR. ROWPIV ) THEN
|
| * |
IF ( ERREST ) THEN
|
| END IF |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON )
|
| * |
ELSE
|
| TRANSP = .FALSE. |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ )
|
| L2TRAN = L2TRAN .AND. ( M .EQ. N ) |
END IF
|
| * |
ELSE
|
| AATMAX = -ONE |
IF ( ERREST ) THEN
|
| AATMIN = BIG |
MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
|
| IF ( ROWPIV .OR. L2TRAN ) THEN |
ELSE
|
| * |
MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
|
| * Compute the row norms, needed to determine row pivoting sequence |
END IF
|
| * (in the case of heavily row weighted A, row pivoting is strongly |
END IF
|
| * advised) and to collect information needed to compare the |
IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
|
| * structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.). |
ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
|
| * |
* .. minimal and optimal sizes of the complex workspace if the
|
| IF ( L2TRAN ) THEN |
* singular values and the left singular vectors are requested
|
| DO 1950 p = 1, M |
IF ( ERREST ) THEN
|
| XSC = ZERO |
MINWRK = N + MAX( LWQP3,LWCON,N+LWQRF,LWSVDJ,LWUNMQRM )
|
| TEMP1 = ONE |
ELSE
|
| CALL ZLASSQ( N, A(p,1), LDA, XSC, TEMP1 ) |
MINWRK = N + MAX( LWQP3, N+LWQRF, LWSVDJ, LWUNMQRM )
|
| * ZLASSQ gets both the ell_2 and the ell_infinity norm |
END IF
|
| * in one pass through the vector |
IF ( LQUERY ) THEN
|
| RWORK(M+N+p) = XSC * SCALEM |
CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
|
| RWORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1)) |
$ LDA, CDUMMY, -1, RDUMMY, -1, IERR )
|
| AATMAX = DMAX1( AATMAX, RWORK(N+p) ) |
LWRK_ZGESVJ = CDUMMY(1)
|
| IF (RWORK(N+p) .NE. ZERO) |
CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
|
| $ AATMIN = DMIN1(AATMIN,RWORK(N+p)) |
$ LDU, CDUMMY, -1, IERR )
|
| 1950 CONTINUE |
LWRK_ZUNMQRM = CDUMMY(1)
|
| ELSE |
IF ( ERREST ) THEN
|
| DO 1904 p = 1, M |
OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF,
|
| RWORK(M+N+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) ) |
$ LWRK_ZGESVJ, LWRK_ZUNMQRM )
|
| AATMAX = DMAX1( AATMAX, RWORK(M+N+p) ) |
ELSE
|
| AATMIN = DMIN1( AATMIN, RWORK(M+N+p) ) |
OPTWRK = N + MAX( LWRK_ZGEQP3, N+LWRK_ZGEQRF,
|
| 1904 CONTINUE |
$ LWRK_ZGESVJ, LWRK_ZUNMQRM )
|
| END IF |
END IF
|
| * |
END IF
|
| END IF |
IF ( L2TRAN .OR. ROWPIV ) THEN
|
| * |
IF ( ERREST ) THEN
|
| * For square matrix A try to determine whether A^* would be better |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON )
|
| * input for the preconditioned Jacobi SVD, with faster convergence. |
ELSE
|
| * The decision is based on an O(N) function of the vector of column |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ )
|
| * and row norms of A, based on the Shannon entropy. This should give |
END IF
|
| * the right choice in most cases when the difference actually matters. |
ELSE
|
| * It may fail and pick the slower converging side. |
IF ( ERREST ) THEN
|
| * |
MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
|
| ENTRA = ZERO |
ELSE
|
| ENTRAT = ZERO |
MINRWRK = MAX( 7, LRWQP3, LRWSVDJ )
|
| IF ( L2TRAN ) THEN |
END IF
|
| * |
END IF
|
| XSC = ZERO |
IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
|
| TEMP1 = ONE |
ELSE
|
| CALL DLASSQ( N, SVA, 1, XSC, TEMP1 ) |
* .. minimal and optimal sizes of the complex workspace if the
|
| TEMP1 = ONE / TEMP1 |
* full SVD is requested
|
| * |
IF ( .NOT. JRACC ) THEN
|
| ENTRA = ZERO |
IF ( ERREST ) THEN
|
| DO 1113 p = 1, N |
MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+N**2+LWCON,
|
| BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1 |
$ 2*N+LWQRF, 2*N+LWQP3,
|
| IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1) |
$ 2*N+N**2+N+LWLQF, 2*N+N**2+N+N**2+LWCON,
|
| 1113 CONTINUE |
$ 2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
|
| ENTRA = - ENTRA / DLOG(DBLE(N)) |
$ 2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
|
| * |
$ N+N**2+LWSVDJ, N+LWUNMQRM )
|
| * Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex. |
ELSE
|
| * It is derived from the diagonal of A^* * A. Do the same with the |
MINWRK = MAX( N+LWQP3, 2*N+N**2+LWCON,
|
| * diagonal of A * A^*, compute the entropy of the corresponding |
$ 2*N+LWQRF, 2*N+LWQP3,
|
| * probability distribution. Note that A * A^* and A^* * A have the |
$ 2*N+N**2+N+LWLQF, 2*N+N**2+N+N**2+LWCON,
|
| * same trace. |
$ 2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV,
|
| * |
$ 2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ,
|
| ENTRAT = ZERO |
$ N+N**2+LWSVDJ, N+LWUNMQRM )
|
| DO 1114 p = N+1, N+M |
END IF
|
| BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1 |
MINIWRK = MINIWRK + N
|
| IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1) |
IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
|
| 1114 CONTINUE |
ELSE
|
| ENTRAT = - ENTRAT / DLOG(DBLE(M)) |
IF ( ERREST ) THEN
|
| * |
MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+LWQRF,
|
| * Analyze the entropies and decide A or A^*. Smaller entropy |
$ 2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
|
| * usually means better input for the algorithm. |
$ N+LWUNMQRM )
|
| * |
ELSE
|
| TRANSP = ( ENTRAT .LT. ENTRA ) |
MINWRK = MAX( N+LWQP3, 2*N+LWQRF,
|
| TRANSP = .TRUE. |
$ 2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR,
|
| * |
$ N+LWUNMQRM )
|
| * If A^* is better than A, take the adjoint of A. |
END IF
|
| * |
IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M
|
| IF ( TRANSP ) THEN |
END IF
|
| * In an optimal implementation, this trivial transpose |
IF ( LQUERY ) THEN
|
| * should be replaced with faster transpose. |
CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
|
| DO 1115 p = 1, N - 1 |
$ LDU, CDUMMY, -1, IERR )
|
| A(p,p) = DCONJG(A(p,p)) |
LWRK_ZUNMQRM = CDUMMY(1)
|
| DO 1116 q = p + 1, N |
CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U,
|
| CTEMP = DCONJG(A(q,p)) |
$ LDU, CDUMMY, -1, IERR )
|
| A(q,p) = DCONJG(A(p,q)) |
LWRK_ZUNMQR = CDUMMY(1)
|
| A(p,q) = CTEMP |
IF ( .NOT. JRACC ) THEN
|
| 1116 CONTINUE |
CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1,
|
| 1115 CONTINUE |
$ RDUMMY, IERR )
|
| A(N,N) = DCONJG(A(N,N)) |
LWRK_ZGEQP3N = CDUMMY(1)
|
| DO 1117 p = 1, N |
CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA,
|
| RWORK(M+N+p) = SVA(p) |
$ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
|
| SVA(p) = RWORK(N+p) |
LWRK_ZGESVJ = CDUMMY(1)
|
| * previously computed row 2-norms are now column 2-norms |
CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA,
|
| * of the transposed matrix |
$ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
|
| 1117 CONTINUE |
LWRK_ZGESVJU = CDUMMY(1)
|
| TEMP1 = AAPP |
CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
|
| AAPP = AATMAX |
$ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
|
| AATMAX = TEMP1 |
LWRK_ZGESVJV = CDUMMY(1)
|
| TEMP1 = AAQQ |
CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
|
| AAQQ = AATMIN |
$ V, LDV, CDUMMY, -1, IERR )
|
| AATMIN = TEMP1 |
LWRK_ZUNMLQ = CDUMMY(1)
|
| KILL = LSVEC |
IF ( ERREST ) THEN
|
| LSVEC = RSVEC |
OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
|
| RSVEC = KILL |
$ 2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
|
| IF ( LSVEC ) N1 = N |
$ 2*N+LWRK_ZGEQP3N,
|
| * |
$ 2*N+N**2+N+LWRK_ZGELQF,
|
| ROWPIV = .TRUE. |
$ 2*N+N**2+N+N**2+LWCON,
|
| END IF |
$ 2*N+N**2+N+LWRK_ZGESVJ,
|
| * |
$ 2*N+N**2+N+LWRK_ZGESVJV,
|
| END IF |
$ 2*N+N**2+N+LWRK_ZUNMQR,
|
| * END IF L2TRAN |
$ 2*N+N**2+N+LWRK_ZUNMLQ,
|
| * |
$ N+N**2+LWRK_ZGESVJU,
|
| * Scale the matrix so that its maximal singular value remains less |
$ N+LWRK_ZUNMQRM )
|
| * than SQRT(BIG) -- the matrix is scaled so that its maximal column |
ELSE
|
| * has Euclidean norm equal to SQRT(BIG/N). The only reason to keep |
OPTWRK = MAX( N+LWRK_ZGEQP3,
|
| * SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and |
$ 2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,
|
| * BLAS routines that, in some implementations, are not capable of |
$ 2*N+LWRK_ZGEQP3N,
|
| * working in the full interval [SFMIN,BIG] and that they may provoke |
$ 2*N+N**2+N+LWRK_ZGELQF,
|
| * overflows in the intermediate results. If the singular values spread |
$ 2*N+N**2+N+N**2+LWCON,
|
| * from SFMIN to BIG, then ZGESVJ will compute them. So, in that case, |
$ 2*N+N**2+N+LWRK_ZGESVJ,
|
| * one should use ZGESVJ instead of ZGEJSV. |
$ 2*N+N**2+N+LWRK_ZGESVJV,
|
| * |
$ 2*N+N**2+N+LWRK_ZUNMQR,
|
| BIG1 = DSQRT( BIG ) |
$ 2*N+N**2+N+LWRK_ZUNMLQ,
|
| TEMP1 = DSQRT( BIG / DBLE(N) ) |
$ N+N**2+LWRK_ZGESVJU,
|
| * |
$ N+LWRK_ZUNMQRM )
|
| CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR ) |
END IF
|
| IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN |
ELSE
|
| AAQQ = ( AAQQ / AAPP ) * TEMP1 |
CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
|
| ELSE |
$ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
|
| AAQQ = ( AAQQ * TEMP1 ) / AAPP |
LWRK_ZGESVJV = CDUMMY(1)
|
| END IF |
CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY,
|
| TEMP1 = TEMP1 * SCALEM |
$ V, LDV, CDUMMY, -1, IERR )
|
| CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR ) |
LWRK_ZUNMQR = CDUMMY(1)
|
| * |
CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
|
| * To undo scaling at the end of this procedure, multiply the |
$ LDU, CDUMMY, -1, IERR )
|
| * computed singular values with USCAL2 / USCAL1. |
LWRK_ZUNMQRM = CDUMMY(1)
|
| * |
IF ( ERREST ) THEN
|
| USCAL1 = TEMP1 |
OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,
|
| USCAL2 = AAPP |
$ 2*N+LWRK_ZGEQRF, 2*N+N**2,
|
| * |
$ 2*N+N**2+LWRK_ZGESVJV,
|
| IF ( L2KILL ) THEN |
$ 2*N+N**2+N+LWRK_ZUNMQR,N+LWRK_ZUNMQRM )
|
| * L2KILL enforces computation of nonzero singular values in |
ELSE
|
| * the restricted range of condition number of the initial A, |
OPTWRK = MAX( N+LWRK_ZGEQP3, 2*N+LWRK_ZGEQRF,
|
| * sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN). |
$ 2*N+N**2, 2*N+N**2+LWRK_ZGESVJV,
|
| XSC = DSQRT( SFMIN ) |
$ 2*N+N**2+N+LWRK_ZUNMQR,
|
| ELSE |
$ N+LWRK_ZUNMQRM )
|
| XSC = SMALL |
END IF
|
| * |
END IF
|
| * Now, if the condition number of A is too big, |
END IF
|
| * sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN, |
IF ( L2TRAN .OR. ROWPIV ) THEN
|
| * as a precaution measure, the full SVD is computed using ZGESVJ |
MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON )
|
| * with accumulated Jacobi rotations. This provides numerically |
ELSE
|
| * more robust computation, at the cost of slightly increased run |
MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON )
|
| * time. Depending on the concrete implementation of BLAS and LAPACK |
END IF
|
| * (i.e. how they behave in presence of extreme ill-conditioning) the |
END IF
|
| * implementor may decide to remove this switch. |
MINWRK = MAX( 2, MINWRK )
|
| IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN |
OPTWRK = MAX( 2, OPTWRK )
|
| JRACC = .TRUE. |
IF ( LWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = - 17
|
| END IF |
IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19
|
| * |
END IF
|
| END IF |
*
|
| IF ( AAQQ .LT. XSC ) THEN |
IF ( INFO .NE. 0 ) THEN
|
| DO 700 p = 1, N |
* #:(
|
| IF ( SVA(p) .LT. XSC ) THEN |
CALL XERBLA( 'ZGEJSV', - INFO )
|
| CALL ZLASET( 'A', M, 1, CZERO, CZERO, A(1,p), LDA ) |
RETURN
|
| SVA(p) = ZERO |
ELSE IF ( LQUERY ) THEN
|
| END IF |
CWORK(1) = OPTWRK
|
| 700 CONTINUE |
CWORK(2) = MINWRK
|
| END IF |
RWORK(1) = MINRWRK
|
| * |
IWORK(1) = MAX( 4, MINIWRK )
|
| * Preconditioning using QR factorization with pivoting |
RETURN
|
| * |
END IF
|
| IF ( ROWPIV ) THEN |
*
|
| * Optional row permutation (Bjoerck row pivoting): |
* Quick return for void matrix (Y3K safe)
|
| * A result by Cox and Higham shows that the Bjoerck's |
* #:)
|
| * row pivoting combined with standard column pivoting |
IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
|
| * has similar effect as Powell-Reid complete pivoting. |
IWORK(1:4) = 0
|
| * The ell-infinity norms of A are made nonincreasing. |
RWORK(1:7) = 0
|
| DO 1952 p = 1, M - 1 |
RETURN
|
| q = IDAMAX( M-p+1, RWORK(M+N+p), 1 ) + p - 1 |
ENDIF
|
| IWORK(2*N+p) = q |
*
|
| IF ( p .NE. q ) THEN |
* Determine whether the matrix U should be M x N or M x M
|
| TEMP1 = RWORK(M+N+p) |
*
|
| RWORK(M+N+p) = RWORK(M+N+q) |
IF ( LSVEC ) THEN
|
| RWORK(M+N+q) = TEMP1 |
N1 = N
|
| END IF |
IF ( LSAME( JOBU, 'F' ) ) N1 = M
|
| 1952 CONTINUE |
END IF
|
| CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 ) |
*
|
| END IF |
* Set numerical parameters
|
| |
*
|
| * |
*! NOTE: Make sure DLAMCH() does not fail on the target architecture.
|
| * End of the preparation phase (scaling, optional sorting and |
*
|
| * transposing, optional flushing of small columns). |
EPSLN = DLAMCH('Epsilon')
|
| * |
SFMIN = DLAMCH('SafeMinimum')
|
| * Preconditioning |
SMALL = SFMIN / EPSLN
|
| * |
BIG = DLAMCH('O')
|
| * If the full SVD is needed, the right singular vectors are computed |
* BIG = ONE / SFMIN
|
| * from a matrix equation, and for that we need theoretical analysis |
*
|
| * of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF. |
* Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
|
| * In all other cases the first RR QRF can be chosen by other criteria |
*
|
| * (eg speed by replacing global with restricted window pivoting, such |
*(!) If necessary, scale SVA() to protect the largest norm from
|
| * as in xGEQPX from TOMS # 782). Good results will be obtained using |
* overflow. It is possible that this scaling pushes the smallest
|
| * xGEQPX with properly (!) chosen numerical parameters. |
* column norm left from the underflow threshold (extreme case).
|
| * Any improvement of ZGEQP3 improves overal performance of ZGEJSV. |
*
|
| * |
SCALEM = ONE / SQRT(DBLE(M)*DBLE(N))
|
| * A * P1 = Q1 * [ R1^* 0]^*: |
NOSCAL = .TRUE.
|
| DO 1963 p = 1, N |
GOSCAL = .TRUE.
|
| * .. all columns are free columns |
DO 1874 p = 1, N
|
| IWORK(p) = 0 |
AAPP = ZERO
|
| 1963 CONTINUE |
AAQQ = ONE
|
| CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LWORK-N, |
CALL ZLASSQ( M, A(1,p), 1, AAPP, AAQQ )
|
| $ RWORK, IERR ) |
IF ( AAPP .GT. BIG ) THEN
|
| * |
INFO = - 9
|
| * The upper triangular matrix R1 from the first QRF is inspected for |
CALL XERBLA( 'ZGEJSV', -INFO )
|
| * rank deficiency and possibilities for deflation, or possible |
RETURN
|
| * ill-conditioning. Depending on the user specified flag L2RANK, |
END IF
|
| * the procedure explores possibilities to reduce the numerical |
AAQQ = SQRT(AAQQ)
|
| * rank by inspecting the computed upper triangular factor. If |
IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN
|
| * L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of |
SVA(p) = AAPP * AAQQ
|
| * A + dA, where ||dA|| <= f(M,N)*EPSLN. |
ELSE
|
| * |
NOSCAL = .FALSE.
|
| NR = 1 |
SVA(p) = AAPP * ( AAQQ * SCALEM )
|
| IF ( L2ABER ) THEN |
IF ( GOSCAL ) THEN
|
| * Standard absolute error bound suffices. All sigma_i with |
GOSCAL = .FALSE.
|
| * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an |
CALL DSCAL( p-1, SCALEM, SVA, 1 )
|
| * agressive enforcement of lower numerical rank by introducing a |
END IF
|
| * backward error of the order of N*EPSLN*||A||. |
END IF
|
| TEMP1 = DSQRT(DBLE(N))*EPSLN |
1874 CONTINUE
|
| DO 3001 p = 2, N |
*
|
| IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN |
IF ( NOSCAL ) SCALEM = ONE
|
| NR = NR + 1 |
*
|
| ELSE |
AAPP = ZERO
|
| GO TO 3002 |
AAQQ = BIG
|
| END IF |
DO 4781 p = 1, N
|
| 3001 CONTINUE |
AAPP = MAX( AAPP, SVA(p) )
|
| 3002 CONTINUE |
IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
|
| ELSE IF ( L2RANK ) THEN |
4781 CONTINUE
|
| * .. similarly as above, only slightly more gentle (less agressive). |
*
|
| * Sudden drop on the diagonal of R1 is used as the criterion for |
* Quick return for zero M x N matrix
|
| * close-to-rank-defficient. |
* #:)
|
| TEMP1 = DSQRT(SFMIN) |
IF ( AAPP .EQ. ZERO ) THEN
|
| DO 3401 p = 2, N |
IF ( LSVEC ) CALL ZLASET( 'G', M, N1, CZERO, CONE, U, LDU )
|
| IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR. |
IF ( RSVEC ) CALL ZLASET( 'G', N, N, CZERO, CONE, V, LDV )
|
| $ ( ABS(A(p,p)) .LT. SMALL ) .OR. |
RWORK(1) = ONE
|
| $ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402 |
RWORK(2) = ONE
|
| NR = NR + 1 |
IF ( ERREST ) RWORK(3) = ONE
|
| 3401 CONTINUE |
IF ( LSVEC .AND. RSVEC ) THEN
|
| 3402 CONTINUE |
RWORK(4) = ONE
|
| * |
RWORK(5) = ONE
|
| ELSE |
END IF
|
| * The goal is high relative accuracy. However, if the matrix |
IF ( L2TRAN ) THEN
|
| * has high scaled condition number the relative accuracy is in |
RWORK(6) = ZERO
|
| * general not feasible. Later on, a condition number estimator |
RWORK(7) = ZERO
|
| * will be deployed to estimate the scaled condition number. |
END IF
|
| * Here we just remove the underflowed part of the triangular |
IWORK(1) = 0
|
| * factor. This prevents the situation in which the code is |
IWORK(2) = 0
|
| * working hard to get the accuracy not warranted by the data. |
IWORK(3) = 0
|
| TEMP1 = DSQRT(SFMIN) |
IWORK(4) = -1
|
| DO 3301 p = 2, N |
RETURN
|
| IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR. |
END IF
|
| $ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302 |
*
|
| NR = NR + 1 |
* Issue warning if denormalized column norms detected. Override the
|
| 3301 CONTINUE |
* high relative accuracy request. Issue licence to kill nonzero columns
|
| 3302 CONTINUE |
* (set them to zero) whose norm is less than sigma_max / BIG (roughly).
|
| * |
* #:(
|
| END IF |
WARNING = 0
|
| * |
IF ( AAQQ .LE. SFMIN ) THEN
|
| ALMORT = .FALSE. |
L2RANK = .TRUE.
|
| IF ( NR .EQ. N ) THEN |
L2KILL = .TRUE.
|
| MAXPRJ = ONE |
WARNING = 1
|
| DO 3051 p = 2, N |
END IF
|
| TEMP1 = ABS(A(p,p)) / SVA(IWORK(p)) |
*
|
| MAXPRJ = DMIN1( MAXPRJ, TEMP1 ) |
* Quick return for one-column matrix
|
| 3051 CONTINUE |
* #:)
|
| IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. |
IF ( N .EQ. 1 ) THEN
|
| END IF |
*
|
| * |
IF ( LSVEC ) THEN
|
| * |
CALL ZLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
|
| SCONDA = - ONE |
CALL ZLACPY( 'A', M, 1, A, LDA, U, LDU )
|
| CONDR1 = - ONE |
* computing all M left singular vectors of the M x 1 matrix
|
| CONDR2 = - ONE |
IF ( N1 .NE. N ) THEN
|
| * |
CALL ZGEQRF( M, N, U,LDU, CWORK, CWORK(N+1),LWORK-N,IERR )
|
| IF ( ERREST ) THEN |
CALL ZUNGQR( M,N1,1, U,LDU,CWORK,CWORK(N+1),LWORK-N,IERR )
|
| IF ( N .EQ. NR ) THEN |
CALL ZCOPY( M, A(1,1), 1, U(1,1), 1 )
|
| IF ( RSVEC ) THEN |
END IF
|
| * .. V is available as workspace |
END IF
|
| CALL ZLACPY( 'U', N, N, A, LDA, V, LDV ) |
IF ( RSVEC ) THEN
|
| DO 3053 p = 1, N |
V(1,1) = CONE
|
| TEMP1 = SVA(IWORK(p)) |
END IF
|
| CALL ZDSCAL( p, ONE/TEMP1, V(1,p), 1 ) |
IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
|
| 3053 CONTINUE |
SVA(1) = SVA(1) / SCALEM
|
| CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1, |
SCALEM = ONE
|
| $ CWORK(N+1), RWORK, IERR ) |
END IF
|
| * |
RWORK(1) = ONE / SCALEM
|
| ELSE IF ( LSVEC ) THEN |
RWORK(2) = ONE
|
| * .. U is available as workspace |
IF ( SVA(1) .NE. ZERO ) THEN
|
| CALL ZLACPY( 'U', N, N, A, LDA, U, LDU ) |
IWORK(1) = 1
|
| DO 3054 p = 1, N |
IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
|
| TEMP1 = SVA(IWORK(p)) |
IWORK(2) = 1
|
| CALL ZDSCAL( p, ONE/TEMP1, U(1,p), 1 ) |
ELSE
|
| 3054 CONTINUE |
IWORK(2) = 0
|
| CALL ZPOCON( 'U', N, U, LDU, ONE, TEMP1, |
END IF
|
| $ CWORK(N+1), RWORK, IERR ) |
ELSE
|
| ELSE |
IWORK(1) = 0
|
| CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N ) |
IWORK(2) = 0
|
| DO 3052 p = 1, N |
END IF
|
| TEMP1 = SVA(IWORK(p)) |
IWORK(3) = 0
|
| CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 ) |
IWORK(4) = -1
|
| 3052 CONTINUE |
IF ( ERREST ) RWORK(3) = ONE
|
| * .. the columns of R are scaled to have unit Euclidean lengths. |
IF ( LSVEC .AND. RSVEC ) THEN
|
| CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1, |
RWORK(4) = ONE
|
| $ CWORK(N+N*N+1), RWORK, IERR ) |
RWORK(5) = ONE
|
| * |
END IF
|
| END IF |
IF ( L2TRAN ) THEN
|
| SCONDA = ONE / DSQRT(TEMP1) |
RWORK(6) = ZERO
|
| * SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). |
RWORK(7) = ZERO
|
| * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA |
END IF
|
| ELSE |
RETURN
|
| SCONDA = - ONE |
*
|
| END IF |
END IF
|
| END IF |
*
|
| * |
TRANSP = .FALSE.
|
| L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) ) |
*
|
| * If there is no violent scaling, artificial perturbation is not needed. |
AATMAX = -ONE
|
| * |
AATMIN = BIG
|
| * Phase 3: |
IF ( ROWPIV .OR. L2TRAN ) THEN
|
| * |
*
|
| IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN |
* Compute the row norms, needed to determine row pivoting sequence
|
| * |
* (in the case of heavily row weighted A, row pivoting is strongly
|
| * Singular Values only |
* advised) and to collect information needed to compare the
|
| * |
* structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.).
|
| * .. transpose A(1:NR,1:N) |
*
|
| DO 1946 p = 1, MIN0( N-1, NR ) |
IF ( L2TRAN ) THEN
|
| CALL ZCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 ) |
DO 1950 p = 1, M
|
| CALL ZLACGV( N-p+1, A(p,p), 1 ) |
XSC = ZERO
|
| 1946 CONTINUE |
TEMP1 = ONE
|
| IF ( NR .EQ. N ) A(N,N) = DCONJG(A(N,N)) |
CALL ZLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
|
| * |
* ZLASSQ gets both the ell_2 and the ell_infinity norm
|
| * The following two DO-loops introduce small relative perturbation |
* in one pass through the vector
|
| * into the strict upper triangle of the lower triangular matrix. |
RWORK(M+p) = XSC * SCALEM
|
| * Small entries below the main diagonal are also changed. |
RWORK(p) = XSC * (SCALEM*SQRT(TEMP1))
|
| * This modification is useful if the computing environment does not |
AATMAX = MAX( AATMAX, RWORK(p) )
|
| * provide/allow FLUSH TO ZERO underflow, for it prevents many |
IF (RWORK(p) .NE. ZERO)
|
| * annoying denormalized numbers in case of strongly scaled matrices. |
$ AATMIN = MIN(AATMIN,RWORK(p))
|
| * The perturbation is structured so that it does not introduce any |
1950 CONTINUE
|
| * new perturbation of the singular values, and it does not destroy |
ELSE
|
| * the job done by the preconditioner. |
DO 1904 p = 1, M
|
| * The licence for this perturbation is in the variable L2PERT, which |
RWORK(M+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
|
| * should be .FALSE. if FLUSH TO ZERO underflow is active. |
AATMAX = MAX( AATMAX, RWORK(M+p) )
|
| * |
AATMIN = MIN( AATMIN, RWORK(M+p) )
|
| IF ( .NOT. ALMORT ) THEN |
1904 CONTINUE
|
| * |
END IF
|
| IF ( L2PERT ) THEN |
*
|
| * XSC = SQRT(SMALL) |
END IF
|
| XSC = EPSLN / DBLE(N) |
*
|
| DO 4947 q = 1, NR |
* For square matrix A try to determine whether A^* would be better
|
| CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO) |
* input for the preconditioned Jacobi SVD, with faster convergence.
|
| DO 4949 p = 1, N |
* The decision is based on an O(N) function of the vector of column
|
| IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) ) |
* and row norms of A, based on the Shannon entropy. This should give
|
| $ .OR. ( p .LT. q ) ) |
* the right choice in most cases when the difference actually matters.
|
| * $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) ) |
* It may fail and pick the slower converging side.
|
| $ A(p,q) = CTEMP |
*
|
| 4949 CONTINUE |
ENTRA = ZERO
|
| 4947 CONTINUE |
ENTRAT = ZERO
|
| ELSE |
IF ( L2TRAN ) THEN
|
| CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, A(1,2),LDA ) |
*
|
| END IF |
XSC = ZERO
|
| * |
TEMP1 = ONE
|
| * .. second preconditioning using the QR factorization |
CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
|
| * |
TEMP1 = ONE / TEMP1
|
| CALL ZGEQRF( N,NR, A,LDA, CWORK, CWORK(N+1),LWORK-N, IERR ) |
*
|
| * |
ENTRA = ZERO
|
| * .. and transpose upper to lower triangular |
DO 1113 p = 1, N
|
| DO 1948 p = 1, NR - 1 |
BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1
|
| CALL ZCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 ) |
IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
|
| CALL ZLACGV( NR-p+1, A(p,p), 1 ) |
1113 CONTINUE
|
| 1948 CONTINUE |
ENTRA = - ENTRA / DLOG(DBLE(N))
|
| * |
*
|
| END IF |
* Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.
|
| * |
* It is derived from the diagonal of A^* * A. Do the same with the
|
| * Row-cyclic Jacobi SVD algorithm with column pivoting |
* diagonal of A * A^*, compute the entropy of the corresponding
|
| * |
* probability distribution. Note that A * A^* and A^* * A have the
|
| * .. again some perturbation (a "background noise") is added |
* same trace.
|
| * to drown denormals |
*
|
| IF ( L2PERT ) THEN |
ENTRAT = ZERO
|
| * XSC = SQRT(SMALL) |
DO 1114 p = 1, M
|
| XSC = EPSLN / DBLE(N) |
BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1
|
| DO 1947 q = 1, NR |
IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
|
| CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO) |
1114 CONTINUE
|
| DO 1949 p = 1, NR |
ENTRAT = - ENTRAT / DLOG(DBLE(M))
|
| IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) ) |
*
|
| $ .OR. ( p .LT. q ) ) |
* Analyze the entropies and decide A or A^*. Smaller entropy
|
| * $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) ) |
* usually means better input for the algorithm.
|
| $ A(p,q) = CTEMP |
*
|
| 1949 CONTINUE |
TRANSP = ( ENTRAT .LT. ENTRA )
|
| 1947 CONTINUE |
*
|
| ELSE |
* If A^* is better than A, take the adjoint of A. This is allowed
|
| CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, A(1,2), LDA ) |
* only for square matrices, M=N.
|
| END IF |
IF ( TRANSP ) THEN
|
| * |
* In an optimal implementation, this trivial transpose
|
| * .. and one-sided Jacobi rotations are started on a lower |
* should be replaced with faster transpose.
|
| * triangular matrix (plus perturbation which is ignored in |
DO 1115 p = 1, N - 1
|
| * the part which destroys triangular form (confusing?!)) |
A(p,p) = CONJG(A(p,p))
|
| * |
DO 1116 q = p + 1, N
|
| CALL ZGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA, |
CTEMP = CONJG(A(q,p))
|
| $ N, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO ) |
A(q,p) = CONJG(A(p,q))
|
| * |
A(p,q) = CTEMP
|
| SCALEM = RWORK(1) |
1116 CONTINUE
|
| NUMRANK = NINT(RWORK(2)) |
1115 CONTINUE
|
| * |
A(N,N) = CONJG(A(N,N))
|
| * |
DO 1117 p = 1, N
|
| ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN |
RWORK(M+p) = SVA(p)
|
| * |
SVA(p) = RWORK(p)
|
| * -> Singular Values and Right Singular Vectors <- |
* previously computed row 2-norms are now column 2-norms
|
| * |
* of the transposed matrix
|
| IF ( ALMORT ) THEN |
1117 CONTINUE
|
| * |
TEMP1 = AAPP
|
| * .. in this case NR equals N |
AAPP = AATMAX
|
| DO 1998 p = 1, NR |
AATMAX = TEMP1
|
| CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) |
TEMP1 = AAQQ
|
| CALL ZLACGV( N-p+1, V(p,p), 1 ) |
AAQQ = AATMIN
|
| 1998 CONTINUE |
AATMIN = TEMP1
|
| CALL ZLASET( 'Upper', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV ) |
KILL = LSVEC
|
| * |
LSVEC = RSVEC
|
| CALL ZGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA, |
RSVEC = KILL
|
| $ CWORK, LWORK, RWORK, LRWORK, INFO ) |
IF ( LSVEC ) N1 = N
|
| SCALEM = RWORK(1) |
*
|
| NUMRANK = NINT(RWORK(2)) |
ROWPIV = .TRUE.
|
| |
END IF
|
| ELSE |
*
|
| * |
END IF
|
| * .. two more QR factorizations ( one QRF is not enough, two require |
* END IF L2TRAN
|
| * accumulated product of Jacobi rotations, three are perfect ) |
*
|
| * |
* Scale the matrix so that its maximal singular value remains less
|
| CALL ZLASET( 'Lower', NR-1,NR-1, CZERO, CZERO, A(2,1), LDA ) |
* than SQRT(BIG) -- the matrix is scaled so that its maximal column
|
| CALL ZGELQF( NR,N, A, LDA, CWORK, CWORK(N+1), LWORK-N, IERR) |
* has Euclidean norm equal to SQRT(BIG/N). The only reason to keep
|
| CALL ZLACPY( 'Lower', NR, NR, A, LDA, V, LDV ) |
* SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and
|
| CALL ZLASET( 'Upper', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV ) |
* BLAS routines that, in some implementations, are not capable of
|
| CALL ZGEQRF( NR, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), |
* working in the full interval [SFMIN,BIG] and that they may provoke
|
| $ LWORK-2*N, IERR ) |
* overflows in the intermediate results. If the singular values spread
|
| DO 8998 p = 1, NR |
* from SFMIN to BIG, then ZGESVJ will compute them. So, in that case,
|
| CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) |
* one should use ZGESVJ instead of ZGEJSV.
|
| CALL ZLACGV( NR-p+1, V(p,p), 1 ) |
* >> change in the April 2016 update: allow bigger range, i.e. the
|
| 8998 CONTINUE |
* largest column is allowed up to BIG/N and ZGESVJ will do the rest.
|
| CALL ZLASET('Upper', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV) |
BIG1 = SQRT( BIG )
|
| * |
TEMP1 = SQRT( BIG / DBLE(N) )
|
| CALL ZGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U, |
* TEMP1 = BIG/DBLE(N)
|
| $ LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO ) |
*
|
| SCALEM = RWORK(1) |
CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
|
| NUMRANK = NINT(RWORK(2)) |
IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
|
| IF ( NR .LT. N ) THEN |
AAQQ = ( AAQQ / AAPP ) * TEMP1
|
| CALL ZLASET( 'A',N-NR, NR, CZERO,CZERO, V(NR+1,1), LDV ) |
ELSE
|
| CALL ZLASET( 'A',NR, N-NR, CZERO,CZERO, V(1,NR+1), LDV ) |
AAQQ = ( AAQQ * TEMP1 ) / AAPP
|
| CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE, V(NR+1,NR+1),LDV ) |
END IF
|
| END IF |
TEMP1 = TEMP1 * SCALEM
|
| * |
CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
|
| CALL ZUNMLQ( 'Left', 'C', N, N, NR, A, LDA, CWORK, |
*
|
| $ V, LDV, CWORK(N+1), LWORK-N, IERR ) |
* To undo scaling at the end of this procedure, multiply the
|
| * |
* computed singular values with USCAL2 / USCAL1.
|
| END IF |
*
|
| * |
USCAL1 = TEMP1
|
| DO 8991 p = 1, N |
USCAL2 = AAPP
|
| CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA ) |
*
|
| 8991 CONTINUE |
IF ( L2KILL ) THEN
|
| CALL ZLACPY( 'All', N, N, A, LDA, V, LDV ) |
* L2KILL enforces computation of nonzero singular values in
|
| * |
* the restricted range of condition number of the initial A,
|
| IF ( TRANSP ) THEN |
* sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN).
|
| CALL ZLACPY( 'All', N, N, V, LDV, U, LDU ) |
XSC = SQRT( SFMIN )
|
| END IF |
ELSE
|
| * |
XSC = SMALL
|
| ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN |
*
|
| * |
* Now, if the condition number of A is too big,
|
| * .. Singular Values and Left Singular Vectors .. |
* sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN,
|
| * |
* as a precaution measure, the full SVD is computed using ZGESVJ
|
| * .. second preconditioning step to avoid need to accumulate |
* with accumulated Jacobi rotations. This provides numerically
|
| * Jacobi rotations in the Jacobi iterations. |
* more robust computation, at the cost of slightly increased run
|
| DO 1965 p = 1, NR |
* time. Depending on the concrete implementation of BLAS and LAPACK
|
| CALL ZCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 ) |
* (i.e. how they behave in presence of extreme ill-conditioning) the
|
| CALL ZLACGV( N-p+1, U(p,p), 1 ) |
* implementor may decide to remove this switch.
|
| 1965 CONTINUE |
IF ( ( AAQQ.LT.SQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
|
| CALL ZLASET( 'Upper', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU ) |
JRACC = .TRUE.
|
| * |
END IF
|
| CALL ZGEQRF( N, NR, U, LDU, CWORK(N+1), CWORK(2*N+1), |
*
|
| $ LWORK-2*N, IERR ) |
END IF
|
| * |
IF ( AAQQ .LT. XSC ) THEN
|
| DO 1967 p = 1, NR - 1 |
DO 700 p = 1, N
|
| CALL ZCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 ) |
IF ( SVA(p) .LT. XSC ) THEN
|
| CALL ZLACGV( N-p+1, U(p,p), 1 ) |
CALL ZLASET( 'A', M, 1, CZERO, CZERO, A(1,p), LDA )
|
| 1967 CONTINUE |
SVA(p) = ZERO
|
| CALL ZLASET( 'Upper', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU ) |
END IF
|
| * |
700 CONTINUE
|
| CALL ZGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A, |
END IF
|
| $ LDA, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO ) |
*
|
| SCALEM = RWORK(1) |
* Preconditioning using QR factorization with pivoting
|
| NUMRANK = NINT(RWORK(2)) |
*
|
| * |
IF ( ROWPIV ) THEN
|
| IF ( NR .LT. M ) THEN |
* Optional row permutation (Bjoerck row pivoting):
|
| CALL ZLASET( 'A', M-NR, NR,CZERO, CZERO, U(NR+1,1), LDU ) |
* A result by Cox and Higham shows that the Bjoerck's
|
| IF ( NR .LT. N1 ) THEN |
* row pivoting combined with standard column pivoting
|
| CALL ZLASET( 'A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU ) |
* has similar effect as Powell-Reid complete pivoting.
|
| CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,U(NR+1,NR+1),LDU ) |
* The ell-infinity norms of A are made nonincreasing.
|
| END IF |
IF ( ( LSVEC .AND. RSVEC ) .AND. .NOT.( JRACC ) ) THEN
|
| END IF |
IWOFF = 2*N
|
| * |
ELSE
|
| CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U, |
IWOFF = N
|
| $ LDU, CWORK(N+1), LWORK-N, IERR ) |
END IF
|
| * |
DO 1952 p = 1, M - 1
|
| IF ( ROWPIV ) |
q = IDAMAX( M-p+1, RWORK(M+p), 1 ) + p - 1
|
| $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
IWORK(IWOFF+p) = q
|
| * |
IF ( p .NE. q ) THEN
|
| DO 1974 p = 1, N1 |
TEMP1 = RWORK(M+p)
|
| XSC = ONE / DZNRM2( M, U(1,p), 1 ) |
RWORK(M+p) = RWORK(M+q)
|
| CALL ZDSCAL( M, XSC, U(1,p), 1 ) |
RWORK(M+q) = TEMP1
|
| 1974 CONTINUE |
END IF
|
| * |
1952 CONTINUE
|
| IF ( TRANSP ) THEN |
CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(IWOFF+1), 1 )
|
| CALL ZLACPY( 'All', N, N, U, LDU, V, LDV ) |
END IF
|
| END IF |
*
|
| * |
* End of the preparation phase (scaling, optional sorting and
|
| ELSE |
* transposing, optional flushing of small columns).
|
| * |
*
|
| * .. Full SVD .. |
* Preconditioning
|
| * |
*
|
| IF ( .NOT. JRACC ) THEN |
* If the full SVD is needed, the right singular vectors are computed
|
| * |
* from a matrix equation, and for that we need theoretical analysis
|
| IF ( .NOT. ALMORT ) THEN |
* of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF.
|
| * |
* In all other cases the first RR QRF can be chosen by other criteria
|
| * Second Preconditioning Step (QRF [with pivoting]) |
* (eg speed by replacing global with restricted window pivoting, such
|
| * Note that the composition of TRANSPOSE, QRF and TRANSPOSE is |
* as in xGEQPX from TOMS # 782). Good results will be obtained using
|
| * equivalent to an LQF CALL. Since in many libraries the QRF |
* xGEQPX with properly (!) chosen numerical parameters.
|
| * seems to be better optimized than the LQF, we do explicit |
* Any improvement of ZGEQP3 improves overal performance of ZGEJSV.
|
| * transpose and use the QRF. This is subject to changes in an |
*
|
| * optimized implementation of ZGEJSV. |
* A * P1 = Q1 * [ R1^* 0]^*:
|
| * |
DO 1963 p = 1, N
|
| DO 1968 p = 1, NR |
* .. all columns are free columns
|
| CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) |
IWORK(p) = 0
|
| CALL ZLACGV( N-p+1, V(p,p), 1 ) |
1963 CONTINUE
|
| 1968 CONTINUE |
CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LWORK-N,
|
| * |
$ RWORK, IERR )
|
| * .. the following two loops perturb small entries to avoid |
*
|
| * denormals in the second QR factorization, where they are |
* The upper triangular matrix R1 from the first QRF is inspected for
|
| * as good as zeros. This is done to avoid painfully slow |
* rank deficiency and possibilities for deflation, or possible
|
| * computation with denormals. The relative size of the perturbation |
* ill-conditioning. Depending on the user specified flag L2RANK,
|
| * is a parameter that can be changed by the implementer. |
* the procedure explores possibilities to reduce the numerical
|
| * This perturbation device will be obsolete on machines with |
* rank by inspecting the computed upper triangular factor. If
|
| * properly implemented arithmetic. |
* L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of
|
| * To switch it off, set L2PERT=.FALSE. To remove it from the |
* A + dA, where ||dA|| <= f(M,N)*EPSLN.
|
| * code, remove the action under L2PERT=.TRUE., leave the ELSE part. |
*
|
| * The following two loops should be blocked and fused with the |
NR = 1
|
| * transposed copy above. |
IF ( L2ABER ) THEN
|
| * |
* Standard absolute error bound suffices. All sigma_i with
|
| IF ( L2PERT ) THEN |
* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
|
| XSC = DSQRT(SMALL) |
* agressive enforcement of lower numerical rank by introducing a
|
| DO 2969 q = 1, NR |
* backward error of the order of N*EPSLN*||A||.
|
| CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO) |
TEMP1 = SQRT(DBLE(N))*EPSLN
|
| DO 2968 p = 1, N |
DO 3001 p = 2, N
|
| IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 ) |
IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
|
| $ .OR. ( p .LT. q ) ) |
NR = NR + 1
|
| * $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) ) |
ELSE
|
| $ V(p,q) = CTEMP |
GO TO 3002
|
| IF ( p .LT. q ) V(p,q) = - V(p,q) |
END IF
|
| 2968 CONTINUE |
3001 CONTINUE
|
| 2969 CONTINUE |
3002 CONTINUE
|
| ELSE |
ELSE IF ( L2RANK ) THEN
|
| CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV ) |
* .. similarly as above, only slightly more gentle (less agressive).
|
| END IF |
* Sudden drop on the diagonal of R1 is used as the criterion for
|
| * |
* close-to-rank-deficient.
|
| * Estimate the row scaled condition number of R1 |
TEMP1 = SQRT(SFMIN)
|
| * (If R1 is rectangular, N > NR, then the condition number |
DO 3401 p = 2, N
|
| * of the leading NR x NR submatrix is estimated.) |
IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
|
| * |
$ ( ABS(A(p,p)) .LT. SMALL ) .OR.
|
| CALL ZLACPY( 'L', NR, NR, V, LDV, CWORK(2*N+1), NR ) |
$ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
|
| DO 3950 p = 1, NR |
NR = NR + 1
|
| TEMP1 = DZNRM2(NR-p+1,CWORK(2*N+(p-1)*NR+p),1) |
3401 CONTINUE
|
| CALL ZDSCAL(NR-p+1,ONE/TEMP1,CWORK(2*N+(p-1)*NR+p),1) |
3402 CONTINUE
|
| 3950 CONTINUE |
*
|
| CALL ZPOCON('Lower',NR,CWORK(2*N+1),NR,ONE,TEMP1, |
ELSE
|
| $ CWORK(2*N+NR*NR+1),RWORK,IERR) |
* The goal is high relative accuracy. However, if the matrix
|
| CONDR1 = ONE / DSQRT(TEMP1) |
* has high scaled condition number the relative accuracy is in
|
| * .. here need a second oppinion on the condition number |
* general not feasible. Later on, a condition number estimator
|
| * .. then assume worst case scenario |
* will be deployed to estimate the scaled condition number.
|
| * R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) |
* Here we just remove the underflowed part of the triangular
|
| * more conservative <=> CONDR1 .LT. SQRT(DBLE(N)) |
* factor. This prevents the situation in which the code is
|
| * |
* working hard to get the accuracy not warranted by the data.
|
| COND_OK = DSQRT(DSQRT(DBLE(NR))) |
TEMP1 = SQRT(SFMIN)
|
| *[TP] COND_OK is a tuning parameter. |
DO 3301 p = 2, N
|
| * |
IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR.
|
| IF ( CONDR1 .LT. COND_OK ) THEN |
$ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
|
| * .. the second QRF without pivoting. Note: in an optimized |
NR = NR + 1
|
| * implementation, this QRF should be implemented as the QRF |
3301 CONTINUE
|
| * of a lower triangular matrix. |
3302 CONTINUE
|
| * R1^* = Q2 * R2 |
*
|
| CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), |
END IF
|
| $ LWORK-2*N, IERR ) |
*
|
| * |
ALMORT = .FALSE.
|
| IF ( L2PERT ) THEN |
IF ( NR .EQ. N ) THEN
|
| XSC = DSQRT(SMALL)/EPSLN |
MAXPRJ = ONE
|
| DO 3959 p = 2, NR |
DO 3051 p = 2, N
|
| DO 3958 q = 1, p - 1 |
TEMP1 = ABS(A(p,p)) / SVA(IWORK(p))
|
| CTEMP=DCMPLX(XSC*DMIN1(ABS(V(p,p)),ABS(V(q,q))), |
MAXPRJ = MIN( MAXPRJ, TEMP1 )
|
| $ ZERO) |
3051 CONTINUE
|
| IF ( ABS(V(q,p)) .LE. TEMP1 ) |
IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
|
| * $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) ) |
END IF
|
| $ V(q,p) = CTEMP |
*
|
| 3958 CONTINUE |
*
|
| 3959 CONTINUE |
SCONDA = - ONE
|
| END IF |
CONDR1 = - ONE
|
| * |
CONDR2 = - ONE
|
| IF ( NR .NE. N ) |
*
|
| $ CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N ) |
IF ( ERREST ) THEN
|
| * .. save ... |
IF ( N .EQ. NR ) THEN
|
| * |
IF ( RSVEC ) THEN
|
| * .. this transposed copy should be better than naive |
* .. V is available as workspace
|
| DO 1969 p = 1, NR - 1 |
CALL ZLACPY( 'U', N, N, A, LDA, V, LDV )
|
| CALL ZCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 ) |
DO 3053 p = 1, N
|
| CALL ZLACGV(NR-p+1, V(p,p), 1 ) |
TEMP1 = SVA(IWORK(p))
|
| 1969 CONTINUE |
CALL ZDSCAL( p, ONE/TEMP1, V(1,p), 1 )
|
| V(NR,NR)=DCONJG(V(NR,NR)) |
3053 CONTINUE
|
| * |
IF ( LSVEC )THEN
|
| CONDR2 = CONDR1 |
CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
|
| * |
$ CWORK(N+1), RWORK, IERR )
|
| ELSE |
ELSE
|
| * |
CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1,
|
| * .. ill-conditioned case: second QRF with pivoting |
$ CWORK, RWORK, IERR )
|
| * Note that windowed pivoting would be equaly good |
END IF
|
| * numerically, and more run-time efficient. So, in |
*
|
| * an optimal implementation, the next call to ZGEQP3 |
ELSE IF ( LSVEC ) THEN
|
| * should be replaced with eg. CALL ZGEQPX (ACM TOMS #782) |
* .. U is available as workspace
|
| * with properly (carefully) chosen parameters. |
CALL ZLACPY( 'U', N, N, A, LDA, U, LDU )
|
| * |
DO 3054 p = 1, N
|
| * R1^* * P2 = Q2 * R2 |
TEMP1 = SVA(IWORK(p))
|
| DO 3003 p = 1, NR |
CALL ZDSCAL( p, ONE/TEMP1, U(1,p), 1 )
|
| IWORK(N+p) = 0 |
3054 CONTINUE
|
| 3003 CONTINUE |
CALL ZPOCON( 'U', N, U, LDU, ONE, TEMP1,
|
| CALL ZGEQP3( N, NR, V, LDV, IWORK(N+1), CWORK(N+1), |
$ CWORK(N+1), RWORK, IERR )
|
| $ CWORK(2*N+1), LWORK-2*N, RWORK, IERR ) |
ELSE
|
| ** CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), |
CALL ZLACPY( 'U', N, N, A, LDA, CWORK, N )
|
| ** $ LWORK-2*N, IERR ) |
*[] CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
|
| IF ( L2PERT ) THEN |
* Change: here index shifted by N to the left, CWORK(1:N)
|
| XSC = DSQRT(SMALL) |
* not needed for SIGMA only computation
|
| DO 3969 p = 2, NR |
DO 3052 p = 1, N
|
| DO 3968 q = 1, p - 1 |
TEMP1 = SVA(IWORK(p))
|
| CTEMP=DCMPLX(XSC*DMIN1(ABS(V(p,p)),ABS(V(q,q))), |
*[] CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 )
|
| $ ZERO) |
CALL ZDSCAL( p, ONE/TEMP1, CWORK((p-1)*N+1), 1 )
|
| IF ( ABS(V(q,p)) .LE. TEMP1 ) |
3052 CONTINUE
|
| * $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) ) |
* .. the columns of R are scaled to have unit Euclidean lengths.
|
| $ V(q,p) = CTEMP |
*[] CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1,
|
| 3968 CONTINUE |
*[] $ CWORK(N+N*N+1), RWORK, IERR )
|
| 3969 CONTINUE |
CALL ZPOCON( 'U', N, CWORK, N, ONE, TEMP1,
|
| END IF |
$ CWORK(N*N+1), RWORK, IERR )
|
| * |
*
|
| CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N ) |
END IF
|
| * |
IF ( TEMP1 .NE. ZERO ) THEN
|
| IF ( L2PERT ) THEN |
SCONDA = ONE / SQRT(TEMP1)
|
| XSC = DSQRT(SMALL) |
ELSE
|
| DO 8970 p = 2, NR |
SCONDA = - ONE
|
| DO 8971 q = 1, p - 1 |
END IF
|
| CTEMP=DCMPLX(XSC*DMIN1(ABS(V(p,p)),ABS(V(q,q))), |
* SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
|
| $ ZERO) |
* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
|
| * V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) ) |
ELSE
|
| V(p,q) = - CTEMP |
SCONDA = - ONE
|
| 8971 CONTINUE |
END IF
|
| 8970 CONTINUE |
END IF
|
| ELSE |
*
|
| CALL ZLASET( 'L',NR-1,NR-1,CZERO,CZERO,V(2,1),LDV ) |
L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. SQRT(BIG1) )
|
| END IF |
* If there is no violent scaling, artificial perturbation is not needed.
|
| * Now, compute R2 = L3 * Q3, the LQ factorization. |
*
|
| CALL ZGELQF( NR, NR, V, LDV, CWORK(2*N+N*NR+1), |
* Phase 3:
|
| $ CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR ) |
*
|
| * .. and estimate the condition number |
IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
|
| CALL ZLACPY( 'L',NR,NR,V,LDV,CWORK(2*N+N*NR+NR+1),NR ) |
*
|
| DO 4950 p = 1, NR |
* Singular Values only
|
| TEMP1 = DZNRM2( p, CWORK(2*N+N*NR+NR+p), NR ) |
*
|
| CALL ZDSCAL( p, ONE/TEMP1, CWORK(2*N+N*NR+NR+p), NR ) |
* .. transpose A(1:NR,1:N)
|
| 4950 CONTINUE |
DO 1946 p = 1, MIN( N-1, NR )
|
| CALL ZPOCON( 'L',NR,CWORK(2*N+N*NR+NR+1),NR,ONE,TEMP1, |
CALL ZCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
|
| $ CWORK(2*N+N*NR+NR+NR*NR+1),RWORK,IERR ) |
CALL ZLACGV( N-p+1, A(p,p), 1 )
|
| CONDR2 = ONE / DSQRT(TEMP1) |
1946 CONTINUE
|
| * |
IF ( NR .EQ. N ) A(N,N) = CONJG(A(N,N))
|
| * |
*
|
| IF ( CONDR2 .GE. COND_OK ) THEN |
* The following two DO-loops introduce small relative perturbation
|
| * .. save the Householder vectors used for Q3 |
* into the strict upper triangle of the lower triangular matrix.
|
| * (this overwrittes the copy of R2, as it will not be |
* Small entries below the main diagonal are also changed.
|
| * needed in this branch, but it does not overwritte the |
* This modification is useful if the computing environment does not
|
| * Huseholder vectors of Q2.). |
* provide/allow FLUSH TO ZERO underflow, for it prevents many
|
| CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N ) |
* annoying denormalized numbers in case of strongly scaled matrices.
|
| * .. and the rest of the information on Q3 is in |
* The perturbation is structured so that it does not introduce any
|
| * WORK(2*N+N*NR+1:2*N+N*NR+N) |
* new perturbation of the singular values, and it does not destroy
|
| END IF |
* the job done by the preconditioner.
|
| * |
* The licence for this perturbation is in the variable L2PERT, which
|
| END IF |
* should be .FALSE. if FLUSH TO ZERO underflow is active.
|
| * |
*
|
| IF ( L2PERT ) THEN |
IF ( .NOT. ALMORT ) THEN
|
| XSC = DSQRT(SMALL) |
*
|
| DO 4968 q = 2, NR |
IF ( L2PERT ) THEN
|
| CTEMP = XSC * V(q,q) |
* XSC = SQRT(SMALL)
|
| DO 4969 p = 1, q - 1 |
XSC = EPSLN / DBLE(N)
|
| * V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) ) |
DO 4947 q = 1, NR
|
| V(p,q) = - CTEMP |
CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
|
| 4969 CONTINUE |
DO 4949 p = 1, N
|
| 4968 CONTINUE |
IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
|
| ELSE |
$ .OR. ( p .LT. q ) )
|
| CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV ) |
* $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
|
| END IF |
$ A(p,q) = CTEMP
|
| * |
4949 CONTINUE
|
| * Second preconditioning finished; continue with Jacobi SVD |
4947 CONTINUE
|
| * The input matrix is lower trinagular. |
ELSE
|
| * |
CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, A(1,2),LDA )
|
| * Recover the right singular vectors as solution of a well |
END IF
|
| * conditioned triangular matrix equation. |
*
|
| * |
* .. second preconditioning using the QR factorization
|
| IF ( CONDR1 .LT. COND_OK ) THEN |
*
|
| * |
CALL ZGEQRF( N,NR, A,LDA, CWORK, CWORK(N+1),LWORK-N, IERR )
|
| CALL ZGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, LDU, |
*
|
| $ CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,RWORK, |
* .. and transpose upper to lower triangular
|
| $ LRWORK, INFO ) |
DO 1948 p = 1, NR - 1
|
| SCALEM = RWORK(1) |
CALL ZCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
|
| NUMRANK = NINT(RWORK(2)) |
CALL ZLACGV( NR-p+1, A(p,p), 1 )
|
| DO 3970 p = 1, NR |
1948 CONTINUE
|
| CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 ) |
*
|
| CALL ZDSCAL( NR, SVA(p), V(1,p), 1 ) |
END IF
|
| 3970 CONTINUE |
*
|
| |
* Row-cyclic Jacobi SVD algorithm with column pivoting
|
| * .. pick the right matrix equation and solve it |
*
|
| * |
* .. again some perturbation (a "background noise") is added
|
| IF ( NR .EQ. N ) THEN |
* to drown denormals
|
| * :)) .. best case, R1 is inverted. The solution of this matrix |
IF ( L2PERT ) THEN
|
| * equation is Q2*V2 = the product of the Jacobi rotations |
* XSC = SQRT(SMALL)
|
| * used in ZGESVJ, premultiplied with the orthogonal matrix |
XSC = EPSLN / DBLE(N)
|
| * from the second QR factorization. |
DO 1947 q = 1, NR
|
| CALL ZTRSM('L','U','N','N', NR,NR,CONE, A,LDA, V,LDV) |
CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
|
| ELSE |
DO 1949 p = 1, NR
|
| * .. R1 is well conditioned, but non-square. Adjoint of R2 |
IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) )
|
| * is inverted to get the product of the Jacobi rotations |
$ .OR. ( p .LT. q ) )
|
| * used in ZGESVJ. The Q-factor from the second QR |
* $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) )
|
| * factorization is then built in explicitly. |
$ A(p,q) = CTEMP
|
| CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1), |
1949 CONTINUE
|
| $ N,V,LDV) |
1947 CONTINUE
|
| IF ( NR .LT. N ) THEN |
ELSE
|
| CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV) |
CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, A(1,2), LDA )
|
| CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV) |
END IF
|
| CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) |
*
|
| END IF |
* .. and one-sided Jacobi rotations are started on a lower
|
| CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), |
* triangular matrix (plus perturbation which is ignored in
|
| $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) |
* the part which destroys triangular form (confusing?!))
|
| END IF |
*
|
| * |
CALL ZGESVJ( 'L', 'N', 'N', NR, NR, A, LDA, SVA,
|
| ELSE IF ( CONDR2 .LT. COND_OK ) THEN |
$ N, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
|
| * |
*
|
| * The matrix R2 is inverted. The solution of the matrix equation |
SCALEM = RWORK(1)
|
| * is Q3^* * V3 = the product of the Jacobi rotations (appplied to |
NUMRANK = NINT(RWORK(2))
|
| * the lower triangular L3 from the LQ factorization of |
*
|
| * R2=L3*Q3), pre-multiplied with the transposed Q3. |
*
|
| CALL ZGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U, |
ELSE IF ( ( RSVEC .AND. ( .NOT. LSVEC ) .AND. ( .NOT. JRACC ) )
|
| $ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, |
$ .OR.
|
| $ RWORK, LRWORK, INFO ) |
$ ( JRACC .AND. ( .NOT. LSVEC ) .AND. ( NR .NE. N ) ) ) THEN
|
| SCALEM = RWORK(1) |
*
|
| NUMRANK = NINT(RWORK(2)) |
* -> Singular Values and Right Singular Vectors <-
|
| DO 3870 p = 1, NR |
*
|
| CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 ) |
IF ( ALMORT ) THEN
|
| CALL ZDSCAL( NR, SVA(p), U(1,p), 1 ) |
*
|
| 3870 CONTINUE |
* .. in this case NR equals N
|
| CALL ZTRSM('L','U','N','N',NR,NR,CONE,CWORK(2*N+1),N, |
DO 1998 p = 1, NR
|
| $ U,LDU) |
CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
|
| * .. apply the permutation from the second QR factorization |
CALL ZLACGV( N-p+1, V(p,p), 1 )
|
| DO 873 q = 1, NR |
1998 CONTINUE
|
| DO 872 p = 1, NR |
CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
|
| CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) |
*
|
| 872 CONTINUE |
CALL ZGESVJ( 'L','U','N', N, NR, V, LDV, SVA, NR, A, LDA,
|
| DO 874 p = 1, NR |
$ CWORK, LWORK, RWORK, LRWORK, INFO )
|
| U(p,q) = CWORK(2*N+N*NR+NR+p) |
SCALEM = RWORK(1)
|
| 874 CONTINUE |
NUMRANK = NINT(RWORK(2))
|
| 873 CONTINUE |
|
| IF ( NR .LT. N ) THEN |
ELSE
|
| CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) |
*
|
| CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) |
* .. two more QR factorizations ( one QRF is not enough, two require
|
| CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) |
* accumulated product of Jacobi rotations, three are perfect )
|
| END IF |
*
|
| CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), |
CALL ZLASET( 'L', NR-1,NR-1, CZERO, CZERO, A(2,1), LDA )
|
| $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
CALL ZGELQF( NR,N, A, LDA, CWORK, CWORK(N+1), LWORK-N, IERR)
|
| ELSE |
CALL ZLACPY( 'L', NR, NR, A, LDA, V, LDV )
|
| * Last line of defense. |
CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV )
|
| * #:( This is a rather pathological case: no scaled condition |
CALL ZGEQRF( NR, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
|
| * improvement after two pivoted QR factorizations. Other |
$ LWORK-2*N, IERR )
|
| * possibility is that the rank revealing QR factorization |
DO 8998 p = 1, NR
|
| * or the condition estimator has failed, or the COND_OK |
CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
|
| * is set very close to ONE (which is unnecessary). Normally, |
CALL ZLACGV( NR-p+1, V(p,p), 1 )
|
| * this branch should never be executed, but in rare cases of |
8998 CONTINUE
|
| * failure of the RRQR or condition estimator, the last line of |
CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV)
|
| * defense ensures that ZGEJSV completes the task. |
*
|
| * Compute the full SVD of L3 using ZGESVJ with explicit |
CALL ZGESVJ( 'L', 'U','N', NR, NR, V,LDV, SVA, NR, U,
|
| * accumulation of Jacobi rotations. |
$ LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
|
| CALL ZGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U, |
SCALEM = RWORK(1)
|
| $ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, |
NUMRANK = NINT(RWORK(2))
|
| $ RWORK, LRWORK, INFO ) |
IF ( NR .LT. N ) THEN
|
| SCALEM = RWORK(1) |
CALL ZLASET( 'A',N-NR, NR, CZERO,CZERO, V(NR+1,1), LDV )
|
| NUMRANK = NINT(RWORK(2)) |
CALL ZLASET( 'A',NR, N-NR, CZERO,CZERO, V(1,NR+1), LDV )
|
| IF ( NR .LT. N ) THEN |
CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE, V(NR+1,NR+1),LDV )
|
| CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) |
END IF
|
| CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) |
*
|
| CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) |
CALL ZUNMLQ( 'L', 'C', N, N, NR, A, LDA, CWORK,
|
| END IF |
$ V, LDV, CWORK(N+1), LWORK-N, IERR )
|
| CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), |
*
|
| $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
END IF
|
| * |
* .. permute the rows of V
|
| CALL ZUNMLQ( 'L', 'C', NR, NR, NR, CWORK(2*N+1), N, |
* DO 8991 p = 1, N
|
| $ CWORK(2*N+N*NR+1), U, LDU, CWORK(2*N+N*NR+NR+1), |
* CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
|
| $ LWORK-2*N-N*NR-NR, IERR ) |
* 8991 CONTINUE
|
| DO 773 q = 1, NR |
* CALL ZLACPY( 'All', N, N, A, LDA, V, LDV )
|
| DO 772 p = 1, NR |
CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
|
| CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) |
*
|
| 772 CONTINUE |
IF ( TRANSP ) THEN
|
| DO 774 p = 1, NR |
CALL ZLACPY( 'A', N, N, V, LDV, U, LDU )
|
| U(p,q) = CWORK(2*N+N*NR+NR+p) |
END IF
|
| 774 CONTINUE |
*
|
| 773 CONTINUE |
ELSE IF ( JRACC .AND. (.NOT. LSVEC) .AND. ( NR.EQ. N ) ) THEN
|
| * |
*
|
| END IF |
CALL ZLASET( 'L', N-1,N-1, CZERO, CZERO, A(2,1), LDA )
|
| * |
*
|
| * Permute the rows of V using the (column) permutation from the |
CALL ZGESVJ( 'U','N','V', N, N, A, LDA, SVA, N, V, LDV,
|
| * first QRF. Also, scale the columns to make them unit in |
$ CWORK, LWORK, RWORK, LRWORK, INFO )
|
| * Euclidean norm. This applies to all cases. |
SCALEM = RWORK(1)
|
| * |
NUMRANK = NINT(RWORK(2))
|
| TEMP1 = DSQRT(DBLE(N)) * EPSLN |
CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK )
|
| DO 1972 q = 1, N |
*
|
| DO 972 p = 1, N |
ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
|
| CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) |
*
|
| 972 CONTINUE |
* .. Singular Values and Left Singular Vectors ..
|
| DO 973 p = 1, N |
*
|
| V(p,q) = CWORK(2*N+N*NR+NR+p) |
* .. second preconditioning step to avoid need to accumulate
|
| 973 CONTINUE |
* Jacobi rotations in the Jacobi iterations.
|
| XSC = ONE / DZNRM2( N, V(1,q), 1 ) |
DO 1965 p = 1, NR
|
| IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
CALL ZCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
|
| $ CALL ZDSCAL( N, XSC, V(1,q), 1 ) |
CALL ZLACGV( N-p+1, U(p,p), 1 )
|
| 1972 CONTINUE |
1965 CONTINUE
|
| * At this moment, V contains the right singular vectors of A. |
CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
|
| * Next, assemble the left singular vector matrix U (M x N). |
*
|
| IF ( NR .LT. M ) THEN |
CALL ZGEQRF( N, NR, U, LDU, CWORK(N+1), CWORK(2*N+1),
|
| CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU) |
$ LWORK-2*N, IERR )
|
| IF ( NR .LT. N1 ) THEN |
*
|
| CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU) |
DO 1967 p = 1, NR - 1
|
| CALL ZLASET('A',M-NR,N1-NR,CZERO,CONE, |
CALL ZCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
|
| $ U(NR+1,NR+1),LDU) |
CALL ZLACGV( N-p+1, U(p,p), 1 )
|
| END IF |
1967 CONTINUE
|
| END IF |
CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
|
| * |
*
|
| * The Q matrix from the first QRF is built into the left singular |
CALL ZGESVJ( 'L', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
|
| * matrix U. This applies to all cases. |
$ LDA, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
|
| * |
SCALEM = RWORK(1)
|
| CALL ZUNMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, CWORK, U, |
NUMRANK = NINT(RWORK(2))
|
| $ LDU, CWORK(N+1), LWORK-N, IERR ) |
*
|
| |
IF ( NR .LT. M ) THEN
|
| * The columns of U are normalized. The cost is O(M*N) flops. |
CALL ZLASET( 'A', M-NR, NR,CZERO, CZERO, U(NR+1,1), LDU )
|
| TEMP1 = DSQRT(DBLE(M)) * EPSLN |
IF ( NR .LT. N1 ) THEN
|
| DO 1973 p = 1, NR |
CALL ZLASET( 'A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU )
|
| XSC = ONE / DZNRM2( M, U(1,p), 1 ) |
CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,U(NR+1,NR+1),LDU )
|
| IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
END IF
|
| $ CALL ZDSCAL( M, XSC, U(1,p), 1 ) |
END IF
|
| 1973 CONTINUE |
*
|
| * |
CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
|
| * If the initial QRF is computed with row pivoting, the left |
$ LDU, CWORK(N+1), LWORK-N, IERR )
|
| * singular vectors must be adjusted. |
*
|
| * |
IF ( ROWPIV )
|
| IF ( ROWPIV ) |
$ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
|
| $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
*
|
| * |
DO 1974 p = 1, N1
|
| ELSE |
XSC = ONE / DZNRM2( M, U(1,p), 1 )
|
| * |
CALL ZDSCAL( M, XSC, U(1,p), 1 )
|
| * .. the initial matrix A has almost orthogonal columns and |
1974 CONTINUE
|
| * the second QRF is not needed |
*
|
| * |
IF ( TRANSP ) THEN
|
| CALL ZLACPY( 'Upper', N, N, A, LDA, CWORK(N+1), N ) |
CALL ZLACPY( 'A', N, N, U, LDU, V, LDV )
|
| IF ( L2PERT ) THEN |
END IF
|
| XSC = DSQRT(SMALL) |
*
|
| DO 5970 p = 2, N |
ELSE
|
| CTEMP = XSC * CWORK( N + (p-1)*N + p ) |
*
|
| DO 5971 q = 1, p - 1 |
* .. Full SVD ..
|
| * CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) / |
*
|
| * $ ABS(CWORK(N+(p-1)*N+q)) ) |
IF ( .NOT. JRACC ) THEN
|
| CWORK(N+(q-1)*N+p)=-CTEMP |
*
|
| 5971 CONTINUE |
IF ( .NOT. ALMORT ) THEN
|
| 5970 CONTINUE |
*
|
| ELSE |
* Second Preconditioning Step (QRF [with pivoting])
|
| CALL ZLASET( 'Lower',N-1,N-1,CZERO,CZERO,CWORK(N+2),N ) |
* Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
|
| END IF |
* equivalent to an LQF CALL. Since in many libraries the QRF
|
| * |
* seems to be better optimized than the LQF, we do explicit
|
| CALL ZGESVJ( 'Upper', 'U', 'N', N, N, CWORK(N+1), N, SVA, |
* transpose and use the QRF. This is subject to changes in an
|
| $ N, U, LDU, CWORK(N+N*N+1), LWORK-N-N*N, RWORK, LRWORK, |
* optimized implementation of ZGEJSV.
|
| $ INFO ) |
*
|
| * |
DO 1968 p = 1, NR
|
| SCALEM = RWORK(1) |
CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
|
| NUMRANK = NINT(RWORK(2)) |
CALL ZLACGV( N-p+1, V(p,p), 1 )
|
| DO 6970 p = 1, N |
1968 CONTINUE
|
| CALL ZCOPY( N, CWORK(N+(p-1)*N+1), 1, U(1,p), 1 ) |
*
|
| CALL ZDSCAL( N, SVA(p), CWORK(N+(p-1)*N+1), 1 ) |
* .. the following two loops perturb small entries to avoid
|
| 6970 CONTINUE |
* denormals in the second QR factorization, where they are
|
| * |
* as good as zeros. This is done to avoid painfully slow
|
| CALL ZTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N, |
* computation with denormals. The relative size of the perturbation
|
| $ CONE, A, LDA, CWORK(N+1), N ) |
* is a parameter that can be changed by the implementer.
|
| DO 6972 p = 1, N |
* This perturbation device will be obsolete on machines with
|
| CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV ) |
* properly implemented arithmetic.
|
| 6972 CONTINUE |
* To switch it off, set L2PERT=.FALSE. To remove it from the
|
| TEMP1 = DSQRT(DBLE(N))*EPSLN |
* code, remove the action under L2PERT=.TRUE., leave the ELSE part.
|
| DO 6971 p = 1, N |
* The following two loops should be blocked and fused with the
|
| XSC = ONE / DZNRM2( N, V(1,p), 1 ) |
* transposed copy above.
|
| IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
*
|
| $ CALL ZDSCAL( N, XSC, V(1,p), 1 ) |
IF ( L2PERT ) THEN
|
| 6971 CONTINUE |
XSC = SQRT(SMALL)
|
| * |
DO 2969 q = 1, NR
|
| * Assemble the left singular vector matrix U (M x N). |
CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
|
| * |
DO 2968 p = 1, N
|
| IF ( N .LT. M ) THEN |
IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
|
| CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U(N+1,1), LDU ) |
$ .OR. ( p .LT. q ) )
|
| IF ( N .LT. N1 ) THEN |
* $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
|
| CALL ZLASET('A',N, N1-N, CZERO, CZERO, U(1,N+1),LDU) |
$ V(p,q) = CTEMP
|
| CALL ZLASET( 'A',M-N,N1-N, CZERO, CONE,U(N+1,N+1),LDU) |
IF ( p .LT. q ) V(p,q) = - V(p,q)
|
| END IF |
2968 CONTINUE
|
| END IF |
2969 CONTINUE
|
| CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U, |
ELSE
|
| $ LDU, CWORK(N+1), LWORK-N, IERR ) |
CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
|
| TEMP1 = DSQRT(DBLE(M))*EPSLN |
END IF
|
| DO 6973 p = 1, N1 |
*
|
| XSC = ONE / DZNRM2( M, U(1,p), 1 ) |
* Estimate the row scaled condition number of R1
|
| IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
* (If R1 is rectangular, N > NR, then the condition number
|
| $ CALL ZDSCAL( M, XSC, U(1,p), 1 ) |
* of the leading NR x NR submatrix is estimated.)
|
| 6973 CONTINUE |
*
|
| * |
CALL ZLACPY( 'L', NR, NR, V, LDV, CWORK(2*N+1), NR )
|
| IF ( ROWPIV ) |
DO 3950 p = 1, NR
|
| $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
TEMP1 = DZNRM2(NR-p+1,CWORK(2*N+(p-1)*NR+p),1)
|
| * |
CALL ZDSCAL(NR-p+1,ONE/TEMP1,CWORK(2*N+(p-1)*NR+p),1)
|
| END IF |
3950 CONTINUE
|
| * |
CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1,
|
| * end of the >> almost orthogonal case << in the full SVD |
$ CWORK(2*N+NR*NR+1),RWORK,IERR)
|
| * |
CONDR1 = ONE / SQRT(TEMP1)
|
| ELSE |
* .. here need a second oppinion on the condition number
|
| * |
* .. then assume worst case scenario
|
| * This branch deploys a preconditioned Jacobi SVD with explicitly |
* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
|
| * accumulated rotations. It is included as optional, mainly for |
* more conservative <=> CONDR1 .LT. SQRT(DBLE(N))
|
| * experimental purposes. It does perfom well, and can also be used. |
*
|
| * In this implementation, this branch will be automatically activated |
COND_OK = SQRT(SQRT(DBLE(NR)))
|
| * if the condition number sigma_max(A) / sigma_min(A) is predicted |
*[TP] COND_OK is a tuning parameter.
|
| * to be greater than the overflow threshold. This is because the |
*
|
| * a posteriori computation of the singular vectors assumes robust |
IF ( CONDR1 .LT. COND_OK ) THEN
|
| * implementation of BLAS and some LAPACK procedures, capable of working |
* .. the second QRF without pivoting. Note: in an optimized
|
| * in presence of extreme values. Since that is not always the case, ... |
* implementation, this QRF should be implemented as the QRF
|
| * |
* of a lower triangular matrix.
|
| DO 7968 p = 1, NR |
* R1^* = Q2 * R2
|
| CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) |
CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
|
| CALL ZLACGV( N-p+1, V(p,p), 1 ) |
$ LWORK-2*N, IERR )
|
| 7968 CONTINUE |
*
|
| * |
IF ( L2PERT ) THEN
|
| IF ( L2PERT ) THEN |
XSC = SQRT(SMALL)/EPSLN
|
| XSC = DSQRT(SMALL/EPSLN) |
DO 3959 p = 2, NR
|
| DO 5969 q = 1, NR |
DO 3958 q = 1, p - 1
|
| CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO) |
CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
|
| DO 5968 p = 1, N |
$ ZERO)
|
| IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 ) |
IF ( ABS(V(q,p)) .LE. TEMP1 )
|
| $ .OR. ( p .LT. q ) ) |
* $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
|
| * $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) ) |
$ V(q,p) = CTEMP
|
| $ V(p,q) = CTEMP |
3958 CONTINUE
|
| IF ( p .LT. q ) V(p,q) = - V(p,q) |
3959 CONTINUE
|
| 5968 CONTINUE |
END IF
|
| 5969 CONTINUE |
*
|
| ELSE |
IF ( NR .NE. N )
|
| CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV ) |
$ CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
|
| END IF |
* .. save ...
|
| |
*
|
| CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), |
* .. this transposed copy should be better than naive
|
| $ LWORK-2*N, IERR ) |
DO 1969 p = 1, NR - 1
|
| CALL ZLACPY( 'L', N, NR, V, LDV, CWORK(2*N+1), N ) |
CALL ZCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
|
| * |
CALL ZLACGV(NR-p+1, V(p,p), 1 )
|
| DO 7969 p = 1, NR |
1969 CONTINUE
|
| CALL ZCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 ) |
V(NR,NR)=CONJG(V(NR,NR))
|
| CALL ZLACGV( NR-p+1, U(p,p), 1 ) |
*
|
| 7969 CONTINUE |
CONDR2 = CONDR1
|
| |
*
|
| IF ( L2PERT ) THEN |
ELSE
|
| XSC = DSQRT(SMALL/EPSLN) |
*
|
| DO 9970 q = 2, NR |
* .. ill-conditioned case: second QRF with pivoting
|
| DO 9971 p = 1, q - 1 |
* Note that windowed pivoting would be equaly good
|
| CTEMP = DCMPLX(XSC * DMIN1(ABS(U(p,p)),ABS(U(q,q))), |
* numerically, and more run-time efficient. So, in
|
| $ ZERO) |
* an optimal implementation, the next call to ZGEQP3
|
| * U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) ) |
* should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)
|
| U(p,q) = - CTEMP |
* with properly (carefully) chosen parameters.
|
| 9971 CONTINUE |
*
|
| 9970 CONTINUE |
* R1^* * P2 = Q2 * R2
|
| ELSE |
DO 3003 p = 1, NR
|
| CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU ) |
IWORK(N+p) = 0
|
| END IF |
3003 CONTINUE
|
| |
CALL ZGEQP3( N, NR, V, LDV, IWORK(N+1), CWORK(N+1),
|
| CALL ZGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA, |
$ CWORK(2*N+1), LWORK-2*N, RWORK, IERR )
|
| $ N, V, LDV, CWORK(2*N+N*NR+1), LWORK-2*N-N*NR, |
** CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
|
| $ RWORK, LRWORK, INFO ) |
** $ LWORK-2*N, IERR )
|
| SCALEM = RWORK(1) |
IF ( L2PERT ) THEN
|
| NUMRANK = NINT(RWORK(2)) |
XSC = SQRT(SMALL)
|
| |
DO 3969 p = 2, NR
|
| IF ( NR .LT. N ) THEN |
DO 3968 q = 1, p - 1
|
| CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) |
CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
|
| CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) |
$ ZERO)
|
| CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV ) |
IF ( ABS(V(q,p)) .LE. TEMP1 )
|
| END IF |
* $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) )
|
| |
$ V(q,p) = CTEMP
|
| CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), |
3968 CONTINUE
|
| $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
3969 CONTINUE
|
| * |
END IF
|
| * Permute the rows of V using the (column) permutation from the |
*
|
| * first QRF. Also, scale the columns to make them unit in |
CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N )
|
| * Euclidean norm. This applies to all cases. |
*
|
| * |
IF ( L2PERT ) THEN
|
| TEMP1 = DSQRT(DBLE(N)) * EPSLN |
XSC = SQRT(SMALL)
|
| DO 7972 q = 1, N |
DO 8970 p = 2, NR
|
| DO 8972 p = 1, N |
DO 8971 q = 1, p - 1
|
| CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) |
CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))),
|
| 8972 CONTINUE |
$ ZERO)
|
| DO 8973 p = 1, N |
* V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) )
|
| V(p,q) = CWORK(2*N+N*NR+NR+p) |
V(p,q) = - CTEMP
|
| 8973 CONTINUE |
8971 CONTINUE
|
| XSC = ONE / DZNRM2( N, V(1,q), 1 ) |
8970 CONTINUE
|
| IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
ELSE
|
| $ CALL ZDSCAL( N, XSC, V(1,q), 1 ) |
CALL ZLASET( 'L',NR-1,NR-1,CZERO,CZERO,V(2,1),LDV )
|
| 7972 CONTINUE |
END IF
|
| * |
* Now, compute R2 = L3 * Q3, the LQ factorization.
|
| * At this moment, V contains the right singular vectors of A. |
CALL ZGELQF( NR, NR, V, LDV, CWORK(2*N+N*NR+1),
|
| * Next, assemble the left singular vector matrix U (M x N). |
$ CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
|
| * |
* .. and estimate the condition number
|
| IF ( NR .LT. M ) THEN |
CALL ZLACPY( 'L',NR,NR,V,LDV,CWORK(2*N+N*NR+NR+1),NR )
|
| CALL ZLASET( 'A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU ) |
DO 4950 p = 1, NR
|
| IF ( NR .LT. N1 ) THEN |
TEMP1 = DZNRM2( p, CWORK(2*N+N*NR+NR+p), NR )
|
| CALL ZLASET('A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU) |
CALL ZDSCAL( p, ONE/TEMP1, CWORK(2*N+N*NR+NR+p), NR )
|
| CALL ZLASET('A',M-NR,N1-NR, CZERO, CONE,U(NR+1,NR+1),LDU) |
4950 CONTINUE
|
| END IF |
CALL ZPOCON( 'L',NR,CWORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
|
| END IF |
$ CWORK(2*N+N*NR+NR+NR*NR+1),RWORK,IERR )
|
| * |
CONDR2 = ONE / SQRT(TEMP1)
|
| CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U, |
*
|
| $ LDU, CWORK(N+1), LWORK-N, IERR ) |
*
|
| * |
IF ( CONDR2 .GE. COND_OK ) THEN
|
| IF ( ROWPIV ) |
* .. save the Householder vectors used for Q3
|
| $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
* (this overwrittes the copy of R2, as it will not be
|
| * |
* needed in this branch, but it does not overwritte the
|
| * |
* Huseholder vectors of Q2.).
|
| END IF |
CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N )
|
| IF ( TRANSP ) THEN |
* .. and the rest of the information on Q3 is in
|
| * .. swap U and V because the procedure worked on A^* |
* WORK(2*N+N*NR+1:2*N+N*NR+N)
|
| DO 6974 p = 1, N |
END IF
|
| CALL ZSWAP( N, U(1,p), 1, V(1,p), 1 ) |
*
|
| 6974 CONTINUE |
END IF
|
| END IF |
*
|
| * |
IF ( L2PERT ) THEN
|
| END IF |
XSC = SQRT(SMALL)
|
| * end of the full SVD |
DO 4968 q = 2, NR
|
| * |
CTEMP = XSC * V(q,q)
|
| * Undo scaling, if necessary (and possible) |
DO 4969 p = 1, q - 1
|
| * |
* V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) )
|
| IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN |
V(p,q) = - CTEMP
|
| CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR ) |
4969 CONTINUE
|
| USCAL1 = ONE |
4968 CONTINUE
|
| USCAL2 = ONE |
ELSE
|
| END IF |
CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
|
| * |
END IF
|
| IF ( NR .LT. N ) THEN |
*
|
| DO 3004 p = NR+1, N |
* Second preconditioning finished; continue with Jacobi SVD
|
| SVA(p) = ZERO |
* The input matrix is lower trinagular.
|
| 3004 CONTINUE |
*
|
| END IF |
* Recover the right singular vectors as solution of a well
|
| * |
* conditioned triangular matrix equation.
|
| RWORK(1) = USCAL2 * SCALEM |
*
|
| RWORK(2) = USCAL1 |
IF ( CONDR1 .LT. COND_OK ) THEN
|
| IF ( ERREST ) RWORK(3) = SCONDA |
*
|
| IF ( LSVEC .AND. RSVEC ) THEN |
CALL ZGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, LDU,
|
| RWORK(4) = CONDR1 |
$ CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,RWORK,
|
| RWORK(5) = CONDR2 |
$ LRWORK, INFO )
|
| END IF |
SCALEM = RWORK(1)
|
| IF ( L2TRAN ) THEN |
NUMRANK = NINT(RWORK(2))
|
| RWORK(6) = ENTRA |
DO 3970 p = 1, NR
|
| RWORK(7) = ENTRAT |
CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 )
|
| END IF |
CALL ZDSCAL( NR, SVA(p), V(1,p), 1 )
|
| * |
3970 CONTINUE
|
| IWORK(1) = NR |
|
| IWORK(2) = NUMRANK |
* .. pick the right matrix equation and solve it
|
| IWORK(3) = WARNING |
*
|
| * |
IF ( NR .EQ. N ) THEN
|
| RETURN |
* :)) .. best case, R1 is inverted. The solution of this matrix
|
| * .. |
* equation is Q2*V2 = the product of the Jacobi rotations
|
| * .. END OF ZGEJSV |
* used in ZGESVJ, premultiplied with the orthogonal matrix
|
| * .. |
* from the second QR factorization.
|
| END |
CALL ZTRSM('L','U','N','N', NR,NR,CONE, A,LDA, V,LDV)
|
| * |
ELSE
|
| |
* .. R1 is well conditioned, but non-square. Adjoint of R2
|
| |
* is inverted to get the product of the Jacobi rotations
|
| |
* used in ZGESVJ. The Q-factor from the second QR
|
| |
* factorization is then built in explicitly.
|
| |
CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
|
| |
$ N,V,LDV)
|
| |
IF ( NR .LT. N ) THEN
|
| |
CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
|
| |
CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
|
| |
CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
|
| |
END IF
|
| |
CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
|
| |
$ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
|
| |
END IF
|
| |
*
|
| |
ELSE IF ( CONDR2 .LT. COND_OK ) THEN
|
| |
*
|
| |
* The matrix R2 is inverted. The solution of the matrix equation
|
| |
* is Q3^* * V3 = the product of the Jacobi rotations (appplied to
|
| |
* the lower triangular L3 from the LQ factorization of
|
| |
* R2=L3*Q3), pre-multiplied with the transposed Q3.
|
| |
CALL ZGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
|
| |
$ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
|
| |
$ RWORK, LRWORK, INFO )
|
| |
SCALEM = RWORK(1)
|
| |
NUMRANK = NINT(RWORK(2))
|
| |
DO 3870 p = 1, NR
|
| |
CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 )
|
| |
CALL ZDSCAL( NR, SVA(p), U(1,p), 1 )
|
| |
3870 CONTINUE
|
| |
CALL ZTRSM('L','U','N','N',NR,NR,CONE,CWORK(2*N+1),N,
|
| |
$ U,LDU)
|
| |
* .. apply the permutation from the second QR factorization
|
| |
DO 873 q = 1, NR
|
| |
DO 872 p = 1, NR
|
| |
CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
|
| |
872 CONTINUE
|
| |
DO 874 p = 1, NR
|
| |
U(p,q) = CWORK(2*N+N*NR+NR+p)
|
| |
874 CONTINUE
|
| |
873 CONTINUE
|
| |
IF ( NR .LT. N ) THEN
|
| |
CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
|
| |
CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
|
| |
CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
|
| |
END IF
|
| |
CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
|
| |
$ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
|
| |
ELSE
|
| |
* Last line of defense.
|
| |
* #:( This is a rather pathological case: no scaled condition
|
| |
* improvement after two pivoted QR factorizations. Other
|
| |
* possibility is that the rank revealing QR factorization
|
| |
* or the condition estimator has failed, or the COND_OK
|
| |
* is set very close to ONE (which is unnecessary). Normally,
|
| |
* this branch should never be executed, but in rare cases of
|
| |
* failure of the RRQR or condition estimator, the last line of
|
| |
* defense ensures that ZGEJSV completes the task.
|
| |
* Compute the full SVD of L3 using ZGESVJ with explicit
|
| |
* accumulation of Jacobi rotations.
|
| |
CALL ZGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
|
| |
$ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR,
|
| |
$ RWORK, LRWORK, INFO )
|
| |
SCALEM = RWORK(1)
|
| |
NUMRANK = NINT(RWORK(2))
|
| |
IF ( NR .LT. N ) THEN
|
| |
CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
|
| |
CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
|
| |
CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
|
| |
END IF
|
| |
CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
|
| |
$ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
|
| |
*
|
| |
CALL ZUNMLQ( 'L', 'C', NR, NR, NR, CWORK(2*N+1), N,
|
| |
$ CWORK(2*N+N*NR+1), U, LDU, CWORK(2*N+N*NR+NR+1),
|
| |
$ LWORK-2*N-N*NR-NR, IERR )
|
| |
DO 773 q = 1, NR
|
| |
DO 772 p = 1, NR
|
| |
CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
|
| |
772 CONTINUE
|
| |
DO 774 p = 1, NR
|
| |
U(p,q) = CWORK(2*N+N*NR+NR+p)
|
| |
774 CONTINUE
|
| |
773 CONTINUE
|
| |
*
|
| |
END IF
|
| |
*
|
| |
* Permute the rows of V using the (column) permutation from the
|
| |
* first QRF. Also, scale the columns to make them unit in
|
| |
* Euclidean norm. This applies to all cases.
|
| |
*
|
| |
TEMP1 = SQRT(DBLE(N)) * EPSLN
|
| |
DO 1972 q = 1, N
|
| |
DO 972 p = 1, N
|
| |
CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
|
| |
972 CONTINUE
|
| |
DO 973 p = 1, N
|
| |
V(p,q) = CWORK(2*N+N*NR+NR+p)
|
| |
973 CONTINUE
|
| |
XSC = ONE / DZNRM2( N, V(1,q), 1 )
|
| |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
|
| |
$ CALL ZDSCAL( N, XSC, V(1,q), 1 )
|
| |
1972 CONTINUE
|
| |
* At this moment, V contains the right singular vectors of A.
|
| |
* Next, assemble the left singular vector matrix U (M x N).
|
| |
IF ( NR .LT. M ) THEN
|
| |
CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU)
|
| |
IF ( NR .LT. N1 ) THEN
|
| |
CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
|
| |
CALL ZLASET('A',M-NR,N1-NR,CZERO,CONE,
|
| |
$ U(NR+1,NR+1),LDU)
|
| |
END IF
|
| |
END IF
|
| |
*
|
| |
* The Q matrix from the first QRF is built into the left singular
|
| |
* matrix U. This applies to all cases.
|
| |
*
|
| |
CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
|
| |
$ LDU, CWORK(N+1), LWORK-N, IERR )
|
| |
|
| |
* The columns of U are normalized. The cost is O(M*N) flops.
|
| |
TEMP1 = SQRT(DBLE(M)) * EPSLN
|
| |
DO 1973 p = 1, NR
|
| |
XSC = ONE / DZNRM2( M, U(1,p), 1 )
|
| |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
|
| |
$ CALL ZDSCAL( M, XSC, U(1,p), 1 )
|
| |
1973 CONTINUE
|
| |
*
|
| |
* If the initial QRF is computed with row pivoting, the left
|
| |
* singular vectors must be adjusted.
|
| |
*
|
| |
IF ( ROWPIV )
|
| |
$ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
|
| |
*
|
| |
ELSE
|
| |
*
|
| |
* .. the initial matrix A has almost orthogonal columns and
|
| |
* the second QRF is not needed
|
| |
*
|
| |
CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N )
|
| |
IF ( L2PERT ) THEN
|
| |
XSC = SQRT(SMALL)
|
| |
DO 5970 p = 2, N
|
| |
CTEMP = XSC * CWORK( N + (p-1)*N + p )
|
| |
DO 5971 q = 1, p - 1
|
| |
* CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) /
|
| |
* $ ABS(CWORK(N+(p-1)*N+q)) )
|
| |
CWORK(N+(q-1)*N+p)=-CTEMP
|
| |
5971 CONTINUE
|
| |
5970 CONTINUE
|
| |
ELSE
|
| |
CALL ZLASET( 'L',N-1,N-1,CZERO,CZERO,CWORK(N+2),N )
|
| |
END IF
|
| |
*
|
| |
CALL ZGESVJ( 'U', 'U', 'N', N, N, CWORK(N+1), N, SVA,
|
| |
$ N, U, LDU, CWORK(N+N*N+1), LWORK-N-N*N, RWORK, LRWORK,
|
| |
$ INFO )
|
| |
*
|
| |
SCALEM = RWORK(1)
|
| |
NUMRANK = NINT(RWORK(2))
|
| |
DO 6970 p = 1, N
|
| |
CALL ZCOPY( N, CWORK(N+(p-1)*N+1), 1, U(1,p), 1 )
|
| |
CALL ZDSCAL( N, SVA(p), CWORK(N+(p-1)*N+1), 1 )
|
| |
6970 CONTINUE
|
| |
*
|
| |
CALL ZTRSM( 'L', 'U', 'N', 'N', N, N,
|
| |
$ CONE, A, LDA, CWORK(N+1), N )
|
| |
DO 6972 p = 1, N
|
| |
CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )
|
| |
6972 CONTINUE
|
| |
TEMP1 = SQRT(DBLE(N))*EPSLN
|
| |
DO 6971 p = 1, N
|
| |
XSC = ONE / DZNRM2( N, V(1,p), 1 )
|
| |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
|
| |
$ CALL ZDSCAL( N, XSC, V(1,p), 1 )
|
| |
6971 CONTINUE
|
| |
*
|
| |
* Assemble the left singular vector matrix U (M x N).
|
| |
*
|
| |
IF ( N .LT. M ) THEN
|
| |
CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U(N+1,1), LDU )
|
| |
IF ( N .LT. N1 ) THEN
|
| |
CALL ZLASET('A',N, N1-N, CZERO, CZERO, U(1,N+1),LDU)
|
| |
CALL ZLASET( 'A',M-N,N1-N, CZERO, CONE,U(N+1,N+1),LDU)
|
| |
END IF
|
| |
END IF
|
| |
CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
|
| |
$ LDU, CWORK(N+1), LWORK-N, IERR )
|
| |
TEMP1 = SQRT(DBLE(M))*EPSLN
|
| |
DO 6973 p = 1, N1
|
| |
XSC = ONE / DZNRM2( M, U(1,p), 1 )
|
| |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
|
| |
$ CALL ZDSCAL( M, XSC, U(1,p), 1 )
|
| |
6973 CONTINUE
|
| |
*
|
| |
IF ( ROWPIV )
|
| |
$ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
|
| |
*
|
| |
END IF
|
| |
*
|
| |
* end of the >> almost orthogonal case << in the full SVD
|
| |
*
|
| |
ELSE
|
| |
*
|
| |
* This branch deploys a preconditioned Jacobi SVD with explicitly
|
| |
* accumulated rotations. It is included as optional, mainly for
|
| |
* experimental purposes. It does perfom well, and can also be used.
|
| |
* In this implementation, this branch will be automatically activated
|
| |
* if the condition number sigma_max(A) / sigma_min(A) is predicted
|
| |
* to be greater than the overflow threshold. This is because the
|
| |
* a posteriori computation of the singular vectors assumes robust
|
| |
* implementation of BLAS and some LAPACK procedures, capable of working
|
| |
* in presence of extreme values, e.g. when the singular values spread from
|
| |
* the underflow to the overflow threshold.
|
| |
*
|
| |
DO 7968 p = 1, NR
|
| |
CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
|
| |
CALL ZLACGV( N-p+1, V(p,p), 1 )
|
| |
7968 CONTINUE
|
| |
*
|
| |
IF ( L2PERT ) THEN
|
| |
XSC = SQRT(SMALL/EPSLN)
|
| |
DO 5969 q = 1, NR
|
| |
CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO)
|
| |
DO 5968 p = 1, N
|
| |
IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 )
|
| |
$ .OR. ( p .LT. q ) )
|
| |
* $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) )
|
| |
$ V(p,q) = CTEMP
|
| |
IF ( p .LT. q ) V(p,q) = - V(p,q)
|
| |
5968 CONTINUE
|
| |
5969 CONTINUE
|
| |
ELSE
|
| |
CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV )
|
| |
END IF
|
| |
|
| |
CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1),
|
| |
$ LWORK-2*N, IERR )
|
| |
CALL ZLACPY( 'L', N, NR, V, LDV, CWORK(2*N+1), N )
|
| |
*
|
| |
DO 7969 p = 1, NR
|
| |
CALL ZCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
|
| |
CALL ZLACGV( NR-p+1, U(p,p), 1 )
|
| |
7969 CONTINUE
|
| |
|
| |
IF ( L2PERT ) THEN
|
| |
XSC = SQRT(SMALL/EPSLN)
|
| |
DO 9970 q = 2, NR
|
| |
DO 9971 p = 1, q - 1
|
| |
CTEMP = DCMPLX(XSC * MIN(ABS(U(p,p)),ABS(U(q,q))),
|
| |
$ ZERO)
|
| |
* U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) )
|
| |
U(p,q) = - CTEMP
|
| |
9971 CONTINUE
|
| |
9970 CONTINUE
|
| |
ELSE
|
| |
CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU )
|
| |
END IF
|
| |
|
| |
CALL ZGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA,
|
| |
$ N, V, LDV, CWORK(2*N+N*NR+1), LWORK-2*N-N*NR,
|
| |
$ RWORK, LRWORK, INFO )
|
| |
SCALEM = RWORK(1)
|
| |
NUMRANK = NINT(RWORK(2))
|
| |
|
| |
IF ( NR .LT. N ) THEN
|
| |
CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
|
| |
CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
|
| |
CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
|
| |
END IF
|
| |
|
| |
CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
|
| |
$ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
|
| |
*
|
| |
* Permute the rows of V using the (column) permutation from the
|
| |
* first QRF. Also, scale the columns to make them unit in
|
| |
* Euclidean norm. This applies to all cases.
|
| |
*
|
| |
TEMP1 = SQRT(DBLE(N)) * EPSLN
|
| |
DO 7972 q = 1, N
|
| |
DO 8972 p = 1, N
|
| |
CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
|
| |
8972 CONTINUE
|
| |
DO 8973 p = 1, N
|
| |
V(p,q) = CWORK(2*N+N*NR+NR+p)
|
| |
8973 CONTINUE
|
| |
XSC = ONE / DZNRM2( N, V(1,q), 1 )
|
| |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
|
| |
$ CALL ZDSCAL( N, XSC, V(1,q), 1 )
|
| |
7972 CONTINUE
|
| |
*
|
| |
* At this moment, V contains the right singular vectors of A.
|
| |
* Next, assemble the left singular vector matrix U (M x N).
|
| |
*
|
| |
IF ( NR .LT. M ) THEN
|
| |
CALL ZLASET( 'A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU )
|
| |
IF ( NR .LT. N1 ) THEN
|
| |
CALL ZLASET('A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU)
|
| |
CALL ZLASET('A',M-NR,N1-NR, CZERO, CONE,U(NR+1,NR+1),LDU)
|
| |
END IF
|
| |
END IF
|
| |
*
|
| |
CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
|
| |
$ LDU, CWORK(N+1), LWORK-N, IERR )
|
| |
*
|
| |
IF ( ROWPIV )
|
| |
$ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 )
|
| |
*
|
| |
*
|
| |
END IF
|
| |
IF ( TRANSP ) THEN
|
| |
* .. swap U and V because the procedure worked on A^*
|
| |
DO 6974 p = 1, N
|
| |
CALL ZSWAP( N, U(1,p), 1, V(1,p), 1 )
|
| |
6974 CONTINUE
|
| |
END IF
|
| |
*
|
| |
END IF
|
| |
* end of the full SVD
|
| |
*
|
| |
* Undo scaling, if necessary (and possible)
|
| |
*
|
| |
IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
|
| |
CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
|
| |
USCAL1 = ONE
|
| |
USCAL2 = ONE
|
| |
END IF
|
| |
*
|
| |
IF ( NR .LT. N ) THEN
|
| |
DO 3004 p = NR+1, N
|
| |
SVA(p) = ZERO
|
| |
3004 CONTINUE
|
| |
END IF
|
| |
*
|
| |
RWORK(1) = USCAL2 * SCALEM
|
| |
RWORK(2) = USCAL1
|
| |
IF ( ERREST ) RWORK(3) = SCONDA
|
| |
IF ( LSVEC .AND. RSVEC ) THEN
|
| |
RWORK(4) = CONDR1
|
| |
RWORK(5) = CONDR2
|
| |
END IF
|
| |
IF ( L2TRAN ) THEN
|
| |
RWORK(6) = ENTRA
|
| |
RWORK(7) = ENTRAT
|
| |
END IF
|
| |
*
|
| |
IWORK(1) = NR
|
| |
IWORK(2) = NUMRANK
|
| |
IWORK(3) = WARNING
|
| |
IF ( TRANSP ) THEN
|
| |
IWORK(4) = 1
|
| |
ELSE
|
| |
IWORK(4) = -1
|
| |
END IF
|
| |
|
| |
*
|
| |
RETURN
|
| |
* ..
|
| |
* .. END OF ZGEJSV
|
| |
* ..
|
| |
END
|
| |
*
|
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgejsv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack