LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgejsv.f
Go to the documentation of this file.
00001 *> \brief \b DGEJSV
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DGEJSV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
00022 *                          M, N, A, LDA, SVA, U, LDU, V, LDV,
00023 *                          WORK, LWORK, IWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       IMPLICIT    NONE
00027 *       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
00031 *      $            WORK( LWORK )
00032 *       INTEGER     IWORK( * )
00033 *       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> DGEJSV computes the singular value decomposition (SVD) of a real M-by-N
00043 *> matrix [A], where M >= N. The SVD of [A] is written as
00044 *>
00045 *>              [A] = [U] * [SIGMA] * [V]^t,
00046 *>
00047 *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
00048 *> diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and
00049 *> [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
00050 *> the singular values of [A]. The columns of [U] and [V] are the left and
00051 *> the right singular vectors of [A], respectively. The matrices [U] and [V]
00052 *> are computed and stored in the arrays U and V, respectively. The diagonal
00053 *> of [SIGMA] is computed and stored in the array SVA.
00054 *> \endverbatim
00055 *
00056 *  Arguments:
00057 *  ==========
00058 *
00059 *> \param[in] JOBA
00060 *> \verbatim
00061 *>          JOBA is CHARACTER*1
00062 *>        Specifies the level of accuracy:
00063 *>       = 'C': This option works well (high relative accuracy) if A = B * D,
00064 *>             with well-conditioned B and arbitrary diagonal matrix D.
00065 *>             The accuracy cannot be spoiled by COLUMN scaling. The
00066 *>             accuracy of the computed output depends on the condition of
00067 *>             B, and the procedure aims at the best theoretical accuracy.
00068 *>             The relative error max_{i=1:N}|d sigma_i| / sigma_i is
00069 *>             bounded by f(M,N)*epsilon* cond(B), independent of D.
00070 *>             The input matrix is preprocessed with the QRF with column
00071 *>             pivoting. This initial preprocessing and preconditioning by
00072 *>             a rank revealing QR factorization is common for all values of
00073 *>             JOBA. Additional actions are specified as follows:
00074 *>       = 'E': Computation as with 'C' with an additional estimate of the
00075 *>             condition number of B. It provides a realistic error bound.
00076 *>       = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
00077 *>             D1, D2, and well-conditioned matrix C, this option gives
00078 *>             higher accuracy than the 'C' option. If the structure of the
00079 *>             input matrix is not known, and relative accuracy is
00080 *>             desirable, then this option is advisable. The input matrix A
00081 *>             is preprocessed with QR factorization with FULL (row and
00082 *>             column) pivoting.
00083 *>       = 'G'  Computation as with 'F' with an additional estimate of the
00084 *>             condition number of B, where A=D*B. If A has heavily weighted
00085 *>             rows, then using this condition number gives too pessimistic
00086 *>             error bound.
00087 *>       = 'A': Small singular values are the noise and the matrix is treated
00088 *>             as numerically rank defficient. The error in the computed
00089 *>             singular values is bounded by f(m,n)*epsilon*||A||.
00090 *>             The computed SVD A = U * S * V^t restores A up to
00091 *>             f(m,n)*epsilon*||A||.
00092 *>             This gives the procedure the licence to discard (set to zero)
00093 *>             all singular values below N*epsilon*||A||.
00094 *>       = 'R': Similar as in 'A'. Rank revealing property of the initial
00095 *>             QR factorization is used do reveal (using triangular factor)
00096 *>             a gap sigma_{r+1} < epsilon * sigma_r in which case the
00097 *>             numerical RANK is declared to be r. The SVD is computed with
00098 *>             absolute error bounds, but more accurately than with 'A'.
00099 *> \endverbatim
00100 *>
00101 *> \param[in] JOBU
00102 *> \verbatim
00103 *>          JOBU is CHARACTER*1
00104 *>        Specifies whether to compute the columns of U:
00105 *>       = 'U': N columns of U are returned in the array U.
00106 *>       = 'F': full set of M left sing. vectors is returned in the array U.
00107 *>       = 'W': U may be used as workspace of length M*N. See the description
00108 *>             of U.
00109 *>       = 'N': U is not computed.
00110 *> \endverbatim
00111 *>
00112 *> \param[in] JOBV
00113 *> \verbatim
00114 *>          JOBV is CHARACTER*1
00115 *>        Specifies whether to compute the matrix V:
00116 *>       = 'V': N columns of V are returned in the array V; Jacobi rotations
00117 *>             are not explicitly accumulated.
00118 *>       = 'J': N columns of V are returned in the array V, but they are
00119 *>             computed as the product of Jacobi rotations. This option is
00120 *>             allowed only if JOBU .NE. 'N', i.e. in computing the full SVD.
00121 *>       = 'W': V may be used as workspace of length N*N. See the description
00122 *>             of V.
00123 *>       = 'N': V is not computed.
00124 *> \endverbatim
00125 *>
00126 *> \param[in] JOBR
00127 *> \verbatim
00128 *>          JOBR is CHARACTER*1
00129 *>        Specifies the RANGE for the singular values. Issues the licence to
00130 *>        set to zero small positive singular values if they are outside
00131 *>        specified range. If A .NE. 0 is scaled so that the largest singular
00132 *>        value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
00133 *>        the licence to kill columns of A whose norm in c*A is less than
00134 *>        DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
00135 *>        where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
00136 *>       = 'N': Do not kill small columns of c*A. This option assumes that
00137 *>             BLAS and QR factorizations and triangular solvers are
00138 *>             implemented to work in that range. If the condition of A
00139 *>             is greater than BIG, use DGESVJ.
00140 *>       = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)]
00141 *>             (roughly, as described above). This option is recommended.
00142 *>                                            ~~~~~~~~~~~~~~~~~~~~~~~~~~~
00143 *>        For computing the singular values in the FULL range [SFMIN,BIG]
00144 *>        use DGESVJ.
00145 *> \endverbatim
00146 *>
00147 *> \param[in] JOBT
00148 *> \verbatim
00149 *>          JOBT is CHARACTER*1
00150 *>        If the matrix is square then the procedure may determine to use
00151 *>        transposed A if A^t seems to be better with respect to convergence.
00152 *>        If the matrix is not square, JOBT is ignored. This is subject to
00153 *>        changes in the future.
00154 *>        The decision is based on two values of entropy over the adjoint
00155 *>        orbit of A^t * A. See the descriptions of WORK(6) and WORK(7).
00156 *>       = 'T': transpose if entropy test indicates possibly faster
00157 *>        convergence of Jacobi process if A^t is taken as input. If A is
00158 *>        replaced with A^t, then the row pivoting is included automatically.
00159 *>       = 'N': do not speculate.
00160 *>        This option can be used to compute only the singular values, or the
00161 *>        full SVD (U, SIGMA and V). For only one set of singular vectors
00162 *>        (U or V), the caller should provide both U and V, as one of the
00163 *>        matrices is used as workspace if the matrix A is transposed.
00164 *>        The implementer can easily remove this constraint and make the
00165 *>        code more complicated. See the descriptions of U and V.
00166 *> \endverbatim
00167 *>
00168 *> \param[in] JOBP
00169 *> \verbatim
00170 *>          JOBP is CHARACTER*1
00171 *>        Issues the licence to introduce structured perturbations to drown
00172 *>        denormalized numbers. This licence should be active if the
00173 *>        denormals are poorly implemented, causing slow computation,
00174 *>        especially in cases of fast convergence (!). For details see [1,2].
00175 *>        For the sake of simplicity, this perturbations are included only
00176 *>        when the full SVD or only the singular values are requested. The
00177 *>        implementer/user can easily add the perturbation for the cases of
00178 *>        computing one set of singular vectors.
00179 *>       = 'P': introduce perturbation
00180 *>       = 'N': do not perturb
00181 *> \endverbatim
00182 *>
00183 *> \param[in] M
00184 *> \verbatim
00185 *>          M is INTEGER
00186 *>         The number of rows of the input matrix A.  M >= 0.
00187 *> \endverbatim
00188 *>
00189 *> \param[in] N
00190 *> \verbatim
00191 *>          N is INTEGER
00192 *>         The number of columns of the input matrix A. M >= N >= 0.
00193 *> \endverbatim
00194 *>
00195 *> \param[in,out] A
00196 *> \verbatim
00197 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00198 *>          On entry, the M-by-N matrix A.
00199 *> \endverbatim
00200 *>
00201 *> \param[in] LDA
00202 *> \verbatim
00203 *>          LDA is INTEGER
00204 *>          The leading dimension of the array A.  LDA >= max(1,M).
00205 *> \endverbatim
00206 *>
00207 *> \param[out] SVA
00208 *> \verbatim
00209 *>          SVA is DOUBLE PRECISION array, dimension (N)
00210 *>          On exit,
00211 *>          - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
00212 *>            computation SVA contains Euclidean column norms of the
00213 *>            iterated matrices in the array A.
00214 *>          - For WORK(1) .NE. WORK(2): The singular values of A are
00215 *>            (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
00216 *>            sigma_max(A) overflows or if small singular values have been
00217 *>            saved from underflow by scaling the input matrix A.
00218 *>          - If JOBR='R' then some of the singular values may be returned
00219 *>            as exact zeros obtained by "set to zero" because they are
00220 *>            below the numerical rank threshold or are denormalized numbers.
00221 *> \endverbatim
00222 *>
00223 *> \param[out] U
00224 *> \verbatim
00225 *>          U is DOUBLE PRECISION array, dimension ( LDU, N )
00226 *>          If JOBU = 'U', then U contains on exit the M-by-N matrix of
00227 *>                         the left singular vectors.
00228 *>          If JOBU = 'F', then U contains on exit the M-by-M matrix of
00229 *>                         the left singular vectors, including an ONB
00230 *>                         of the orthogonal complement of the Range(A).
00231 *>          If JOBU = 'W'  .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
00232 *>                         then U is used as workspace if the procedure
00233 *>                         replaces A with A^t. In that case, [V] is computed
00234 *>                         in U as left singular vectors of A^t and then
00235 *>                         copied back to the V array. This 'W' option is just
00236 *>                         a reminder to the caller that in this case U is
00237 *>                         reserved as workspace of length N*N.
00238 *>          If JOBU = 'N'  U is not referenced.
00239 *> \endverbatim
00240 *>
00241 *> \param[in] LDU
00242 *> \verbatim
00243 *>          LDU is INTEGER
00244 *>          The leading dimension of the array U,  LDU >= 1.
00245 *>          IF  JOBU = 'U' or 'F' or 'W',  then LDU >= M.
00246 *> \endverbatim
00247 *>
00248 *> \param[out] V
00249 *> \verbatim
00250 *>          V is DOUBLE PRECISION array, dimension ( LDV, N )
00251 *>          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
00252 *>                         the right singular vectors;
00253 *>          If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
00254 *>                         then V is used as workspace if the pprocedure
00255 *>                         replaces A with A^t. In that case, [U] is computed
00256 *>                         in V as right singular vectors of A^t and then
00257 *>                         copied back to the U array. This 'W' option is just
00258 *>                         a reminder to the caller that in this case V is
00259 *>                         reserved as workspace of length N*N.
00260 *>          If JOBV = 'N'  V is not referenced.
00261 *> \endverbatim
00262 *>
00263 *> \param[in] LDV
00264 *> \verbatim
00265 *>          LDV is INTEGER
00266 *>          The leading dimension of the array V,  LDV >= 1.
00267 *>          If JOBV = 'V' or 'J' or 'W', then LDV >= N.
00268 *> \endverbatim
00269 *>
00270 *> \param[out] WORK
00271 *> \verbatim
00272 *>          WORK is DOUBLE PRECISION array, dimension at least LWORK.
00273 *>          On exit, if N.GT.0 .AND. M.GT.0 (else not referenced), 
00274 *>          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
00275 *>                    that SCALE*SVA(1:N) are the computed singular values
00276 *>                    of A. (See the description of SVA().)
00277 *>          WORK(2) = See the description of WORK(1).
00278 *>          WORK(3) = SCONDA is an estimate for the condition number of
00279 *>                    column equilibrated A. (If JOBA .EQ. 'E' or 'G')
00280 *>                    SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
00281 *>                    It is computed using DPOCON. It holds
00282 *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
00283 *>                    where R is the triangular factor from the QRF of A.
00284 *>                    However, if R is truncated and the numerical rank is
00285 *>                    determined to be strictly smaller than N, SCONDA is
00286 *>                    returned as -1, thus indicating that the smallest
00287 *>                    singular values might be lost.
00288 *>
00289 *>          If full SVD is needed, the following two condition numbers are
00290 *>          useful for the analysis of the algorithm. They are provied for
00291 *>          a developer/implementer who is familiar with the details of
00292 *>          the method.
00293 *>
00294 *>          WORK(4) = an estimate of the scaled condition number of the
00295 *>                    triangular factor in the first QR factorization.
00296 *>          WORK(5) = an estimate of the scaled condition number of the
00297 *>                    triangular factor in the second QR factorization.
00298 *>          The following two parameters are computed if JOBT .EQ. 'T'.
00299 *>          They are provided for a developer/implementer who is familiar
00300 *>          with the details of the method.
00301 *>
00302 *>          WORK(6) = the entropy of A^t*A :: this is the Shannon entropy
00303 *>                    of diag(A^t*A) / Trace(A^t*A) taken as point in the
00304 *>                    probability simplex.
00305 *>          WORK(7) = the entropy of A*A^t.
00306 *> \endverbatim
00307 *>
00308 *> \param[in] LWORK
00309 *> \verbatim
00310 *>          LWORK is INTEGER
00311 *>          Length of WORK to confirm proper allocation of work space.
00312 *>          LWORK depends on the job:
00313 *>
00314 *>          If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
00315 *>            -> .. no scaled condition estimate required (JOBE.EQ.'N'):
00316 *>               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
00317 *>               ->> For optimal performance (blocked code) the optimal value
00318 *>               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
00319 *>               block size for DGEQP3 and DGEQRF.
00320 *>               In general, optimal LWORK is computed as 
00321 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).        
00322 *>            -> .. an estimate of the scaled condition number of A is
00323 *>               required (JOBA='E', 'G'). In this case, LWORK is the maximum
00324 *>               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7).
00325 *>               ->> For optimal performance (blocked code) the optimal value 
00326 *>               is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7).
00327 *>               In general, the optimal length LWORK is computed as
00328 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 
00329 *>                                                     N+N*N+LWORK(DPOCON),7).
00330 *>
00331 *>          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
00332 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
00333 *>            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
00334 *>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ,
00335 *>               DORMLQ. In general, the optimal length LWORK is computed as
00336 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), 
00337 *>                       N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
00338 *>
00339 *>          If SIGMA and the left singular vectors are needed
00340 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
00341 *>            -> For optimal performance:
00342 *>               if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
00343 *>               if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
00344 *>               where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
00345 *>               In general, the optimal length LWORK is computed as
00346 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
00347 *>                        2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). 
00348 *>               Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or 
00349 *>               M*NB (for JOBU.EQ.'F').
00350 *>               
00351 *>          If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and 
00352 *>            -> if JOBV.EQ.'V'  
00353 *>               the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). 
00354 *>            -> if JOBV.EQ.'J' the minimal requirement is 
00355 *>               LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
00356 *>            -> For optimal performance, LWORK should be additionally
00357 *>               larger than N+M*NB, where NB is the optimal block size
00358 *>               for DORMQR.
00359 *> \endverbatim
00360 *>
00361 *> \param[out] IWORK
00362 *> \verbatim
00363 *>          IWORK is INTEGER array, dimension M+3*N.
00364 *>          On exit,
00365 *>          IWORK(1) = the numerical rank determined after the initial
00366 *>                     QR factorization with pivoting. See the descriptions
00367 *>                     of JOBA and JOBR.
00368 *>          IWORK(2) = the number of the computed nonzero singular values
00369 *>          IWORK(3) = if nonzero, a warning message:
00370 *>                     If IWORK(3).EQ.1 then some of the column norms of A
00371 *>                     were denormalized floats. The requested high accuracy
00372 *>                     is not warranted by the data.
00373 *> \endverbatim
00374 *>
00375 *> \param[out] INFO
00376 *> \verbatim
00377 *>          INFO is INTEGER
00378 *>           < 0  : if INFO = -i, then the i-th argument had an illegal value.
00379 *>           = 0 :  successfull exit;
00380 *>           > 0 :  DGEJSV  did not converge in the maximal allowed number
00381 *>                  of sweeps. The computed values may be inaccurate.
00382 *> \endverbatim
00383 *
00384 *  Authors:
00385 *  ========
00386 *
00387 *> \author Univ. of Tennessee 
00388 *> \author Univ. of California Berkeley 
00389 *> \author Univ. of Colorado Denver 
00390 *> \author NAG Ltd. 
00391 *
00392 *> \date November 2011
00393 *
00394 *> \ingroup doubleGEcomputational
00395 *
00396 *> \par Further Details:
00397 *  =====================
00398 *>
00399 *> \verbatim
00400 *>
00401 *>  DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses DGEQP3,
00402 *>  DGEQRF, and DGELQF as preprocessors and preconditioners. Optionally, an
00403 *>  additional row pivoting can be used as a preprocessor, which in some
00404 *>  cases results in much higher accuracy. An example is matrix A with the
00405 *>  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
00406 *>  diagonal matrices and C is well-conditioned matrix. In that case, complete
00407 *>  pivoting in the first QR factorizations provides accuracy dependent on the
00408 *>  condition number of C, and independent of D1, D2. Such higher accuracy is
00409 *>  not completely understood theoretically, but it works well in practice.
00410 *>  Further, if A can be written as A = B*D, with well-conditioned B and some
00411 *>  diagonal D, then the high accuracy is guaranteed, both theoretically and
00412 *>  in software, independent of D. For more details see [1], [2].
00413 *>     The computational range for the singular values can be the full range
00414 *>  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
00415 *>  & LAPACK routines called by DGEJSV are implemented to work in that range.
00416 *>  If that is not the case, then the restriction for safe computation with
00417 *>  the singular values in the range of normalized IEEE numbers is that the
00418 *>  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
00419 *>  overflow. This code (DGEJSV) is best used in this restricted range,
00420 *>  meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
00421 *>  returned as zeros. See JOBR for details on this.
00422 *>     Further, this implementation is somewhat slower than the one described
00423 *>  in [1,2] due to replacement of some non-LAPACK components, and because
00424 *>  the choice of some tuning parameters in the iterative part (DGESVJ) is
00425 *>  left to the implementer on a particular machine.
00426 *>     The rank revealing QR factorization (in this code: DGEQP3) should be
00427 *>  implemented as in [3]. We have a new version of DGEQP3 under development
00428 *>  that is more robust than the current one in LAPACK, with a cleaner cut in
00429 *>  rank defficient cases. It will be available in the SIGMA library [4].
00430 *>  If M is much larger than N, it is obvious that the inital QRF with
00431 *>  column pivoting can be preprocessed by the QRF without pivoting. That
00432 *>  well known trick is not used in DGEJSV because in some cases heavy row
00433 *>  weighting can be treated with complete pivoting. The overhead in cases
00434 *>  M much larger than N is then only due to pivoting, but the benefits in
00435 *>  terms of accuracy have prevailed. The implementer/user can incorporate
00436 *>  this extra QRF step easily. The implementer can also improve data movement
00437 *>  (matrix transpose, matrix copy, matrix transposed copy) - this
00438 *>  implementation of DGEJSV uses only the simplest, naive data movement.
00439 *> \endverbatim
00440 *
00441 *> \par Contributors:
00442 *  ==================
00443 *>
00444 *>  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
00445 *
00446 *> \par References:
00447 *  ================
00448 *>
00449 *> \verbatim
00450 *>
00451 *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
00452 *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
00453 *>     LAPACK Working note 169.
00454 *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
00455 *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
00456 *>     LAPACK Working note 170.
00457 *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
00458 *>     factorization software - a case study.
00459 *>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
00460 *>     LAPACK Working note 176.
00461 *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
00462 *>     QSVD, (H,K)-SVD computations.
00463 *>     Department of Mathematics, University of Zagreb, 2008.
00464 *> \endverbatim
00465 *
00466 *>  \par Bugs, examples and comments:
00467 *   =================================
00468 *>
00469 *>  Please report all bugs and send interesting examples and/or comments to
00470 *>  drmac@math.hr. Thank you.
00471 *>
00472 *  =====================================================================
00473       SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
00474      $                   M, N, A, LDA, SVA, U, LDU, V, LDV,
00475      $                   WORK, LWORK, IWORK, INFO )
00476 *
00477 *  -- LAPACK computational routine (version 3.4.0) --
00478 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00479 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00480 *     November 2011
00481 *
00482 *     .. Scalar Arguments ..
00483       IMPLICIT    NONE
00484       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
00485 *     ..
00486 *     .. Array Arguments ..
00487       DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
00488      $            WORK( LWORK )
00489       INTEGER     IWORK( * )
00490       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
00491 *     ..
00492 *
00493 *  ===========================================================================
00494 *
00495 *     .. Local Parameters ..
00496       DOUBLE PRECISION   ZERO,  ONE
00497       PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
00498 *     ..
00499 *     .. Local Scalars ..
00500       DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
00501      $        CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,  MAXPRJ, SCALEM,
00502      $        SCONDA, SFMIN,  SMALL,  TEMP1,  USCAL1, USCAL2, XSC
00503       INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
00504       LOGICAL ALMORT, DEFR,   ERREST, GOSCAL, JRACC,  KILL,   LSVEC,
00505      $        L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
00506      $        NOSCAL, ROWPIV, RSVEC,  TRANSP
00507 *     ..
00508 *     .. Intrinsic Functions ..
00509       INTRINSIC DABS,  DLOG, DMAX1, DMIN1, DBLE,
00510      $          MAX0, MIN0, IDNINT,  DSIGN,  DSQRT
00511 *     ..
00512 *     .. External Functions ..
00513       DOUBLE PRECISION  DLAMCH, DNRM2
00514       INTEGER   IDAMAX
00515       LOGICAL   LSAME
00516       EXTERNAL  IDAMAX, LSAME, DLAMCH, DNRM2
00517 *     ..
00518 *     .. External Subroutines ..
00519       EXTERNAL  DCOPY,  DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL,
00520      $          DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ,
00521      $          DORMQR, DPOCON, DSCAL,  DSWAP,  DTRSM,  XERBLA
00522 *
00523       EXTERNAL  DGESVJ
00524 *     ..
00525 *
00526 *     Test the input arguments
00527 *
00528       LSVEC  = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
00529       JRACC  = LSAME( JOBV, 'J' )
00530       RSVEC  = LSAME( JOBV, 'V' ) .OR. JRACC
00531       ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
00532       L2RANK = LSAME( JOBA, 'R' )
00533       L2ABER = LSAME( JOBA, 'A' )
00534       ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
00535       L2TRAN = LSAME( JOBT, 'T' )
00536       L2KILL = LSAME( JOBR, 'R' )
00537       DEFR   = LSAME( JOBR, 'N' )
00538       L2PERT = LSAME( JOBP, 'P' )
00539 *
00540       IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
00541      $     ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
00542          INFO = - 1
00543       ELSE IF ( .NOT.( LSVEC  .OR. LSAME( JOBU, 'N' ) .OR.
00544      $                             LSAME( JOBU, 'W' )) ) THEN
00545          INFO = - 2
00546       ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
00547      $   LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN
00548          INFO = - 3
00549       ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) )    THEN
00550          INFO = - 4
00551       ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN
00552          INFO = - 5
00553       ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
00554          INFO = - 6
00555       ELSE IF ( M .LT. 0 ) THEN
00556          INFO = - 7
00557       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
00558          INFO = - 8
00559       ELSE IF ( LDA .LT. M ) THEN
00560          INFO = - 10
00561       ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
00562          INFO = - 13
00563       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
00564          INFO = - 14
00565       ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
00566      &                           (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR.
00567      & (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND.
00568      &                         (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR.
00569      & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1)))
00570      & .OR.
00571      & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1)))
00572      & .OR.
00573      & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. 
00574      &                          (LWORK.LT.MAX0(2*M+N,6*N+2*N*N)))
00575      & .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND.
00576      &                          LWORK.LT.MAX0(2*M+N,4*N+N*N,2*N+N*N+6)))
00577      &   THEN
00578          INFO = - 17
00579       ELSE
00580 *        #:)
00581          INFO = 0
00582       END IF
00583 *
00584       IF ( INFO .NE. 0 ) THEN
00585 *       #:(
00586          CALL XERBLA( 'DGEJSV', - INFO )
00587          RETURN
00588       END IF
00589 *
00590 *     Quick return for void matrix (Y3K safe)
00591 * #:)
00592       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
00593 *
00594 *     Determine whether the matrix U should be M x N or M x M
00595 *
00596       IF ( LSVEC ) THEN
00597          N1 = N
00598          IF ( LSAME( JOBU, 'F' ) ) N1 = M
00599       END IF
00600 *
00601 *     Set numerical parameters
00602 *
00603 *!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
00604 *
00605       EPSLN = DLAMCH('Epsilon')
00606       SFMIN = DLAMCH('SafeMinimum')
00607       SMALL = SFMIN / EPSLN
00608       BIG   = DLAMCH('O')
00609 *     BIG   = ONE / SFMIN
00610 *
00611 *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
00612 *
00613 *(!)  If necessary, scale SVA() to protect the largest norm from
00614 *     overflow. It is possible that this scaling pushes the smallest
00615 *     column norm left from the underflow threshold (extreme case).
00616 *
00617       SCALEM  = ONE / DSQRT(DBLE(M)*DBLE(N))
00618       NOSCAL  = .TRUE.
00619       GOSCAL  = .TRUE.
00620       DO 1874 p = 1, N
00621          AAPP = ZERO
00622          AAQQ = ONE
00623          CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ )
00624          IF ( AAPP .GT. BIG ) THEN
00625             INFO = - 9
00626             CALL XERBLA( 'DGEJSV', -INFO )
00627             RETURN
00628          END IF
00629          AAQQ = DSQRT(AAQQ)
00630          IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL  ) THEN
00631             SVA(p)  = AAPP * AAQQ
00632          ELSE
00633             NOSCAL  = .FALSE.
00634             SVA(p)  = AAPP * ( AAQQ * SCALEM )
00635             IF ( GOSCAL ) THEN
00636                GOSCAL = .FALSE.
00637                CALL DSCAL( p-1, SCALEM, SVA, 1 )
00638             END IF
00639          END IF
00640  1874 CONTINUE
00641 *
00642       IF ( NOSCAL ) SCALEM = ONE
00643 *
00644       AAPP = ZERO
00645       AAQQ = BIG
00646       DO 4781 p = 1, N
00647          AAPP = DMAX1( AAPP, SVA(p) )
00648          IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) )
00649  4781 CONTINUE
00650 *
00651 *     Quick return for zero M x N matrix
00652 * #:)
00653       IF ( AAPP .EQ. ZERO ) THEN
00654          IF ( LSVEC ) CALL DLASET( 'G', M, N1, ZERO, ONE, U, LDU )
00655          IF ( RSVEC ) CALL DLASET( 'G', N, N,  ZERO, ONE, V, LDV )
00656          WORK(1) = ONE
00657          WORK(2) = ONE
00658          IF ( ERREST ) WORK(3) = ONE
00659          IF ( LSVEC .AND. RSVEC ) THEN
00660             WORK(4) = ONE
00661             WORK(5) = ONE
00662          END IF
00663          IF ( L2TRAN ) THEN
00664             WORK(6) = ZERO
00665             WORK(7) = ZERO
00666          END IF
00667          IWORK(1) = 0
00668          IWORK(2) = 0
00669          IWORK(3) = 0
00670          RETURN
00671       END IF
00672 *
00673 *     Issue warning if denormalized column norms detected. Override the
00674 *     high relative accuracy request. Issue licence to kill columns
00675 *     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
00676 * #:(
00677       WARNING = 0
00678       IF ( AAQQ .LE. SFMIN ) THEN
00679          L2RANK = .TRUE.
00680          L2KILL = .TRUE.
00681          WARNING = 1
00682       END IF
00683 *
00684 *     Quick return for one-column matrix
00685 * #:)
00686       IF ( N .EQ. 1 ) THEN
00687 *
00688          IF ( LSVEC ) THEN
00689             CALL DLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
00690             CALL DLACPY( 'A', M, 1, A, LDA, U, LDU )
00691 *           computing all M left singular vectors of the M x 1 matrix
00692             IF ( N1 .NE. N  ) THEN
00693                CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
00694                CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
00695                CALL DCOPY( M, A(1,1), 1, U(1,1), 1 )
00696             END IF
00697          END IF
00698          IF ( RSVEC ) THEN
00699              V(1,1) = ONE
00700          END IF
00701          IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
00702             SVA(1)  = SVA(1) / SCALEM
00703             SCALEM  = ONE
00704          END IF
00705          WORK(1) = ONE / SCALEM
00706          WORK(2) = ONE
00707          IF ( SVA(1) .NE. ZERO ) THEN
00708             IWORK(1) = 1
00709             IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
00710                IWORK(2) = 1
00711             ELSE
00712                IWORK(2) = 0
00713             END IF
00714          ELSE
00715             IWORK(1) = 0
00716             IWORK(2) = 0
00717          END IF
00718          IF ( ERREST ) WORK(3) = ONE
00719          IF ( LSVEC .AND. RSVEC ) THEN
00720             WORK(4) = ONE
00721             WORK(5) = ONE
00722          END IF
00723          IF ( L2TRAN ) THEN
00724             WORK(6) = ZERO
00725             WORK(7) = ZERO
00726          END IF
00727          RETURN
00728 *
00729       END IF
00730 *
00731       TRANSP = .FALSE.
00732       L2TRAN = L2TRAN .AND. ( M .EQ. N )
00733 *
00734       AATMAX = -ONE
00735       AATMIN =  BIG
00736       IF ( ROWPIV .OR. L2TRAN ) THEN
00737 *
00738 *     Compute the row norms, needed to determine row pivoting sequence
00739 *     (in the case of heavily row weighted A, row pivoting is strongly
00740 *     advised) and to collect information needed to compare the
00741 *     structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.).
00742 *
00743          IF ( L2TRAN ) THEN
00744             DO 1950 p = 1, M
00745                XSC   = ZERO
00746                TEMP1 = ONE
00747                CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
00748 *              DLASSQ gets both the ell_2 and the ell_infinity norm
00749 *              in one pass through the vector
00750                WORK(M+N+p)  = XSC * SCALEM
00751                WORK(N+p)    = XSC * (SCALEM*DSQRT(TEMP1))
00752                AATMAX = DMAX1( AATMAX, WORK(N+p) )
00753                IF (WORK(N+p) .NE. ZERO) AATMIN = DMIN1(AATMIN,WORK(N+p))
00754  1950       CONTINUE
00755          ELSE
00756             DO 1904 p = 1, M
00757                WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) )
00758                AATMAX = DMAX1( AATMAX, WORK(M+N+p) )
00759                AATMIN = DMIN1( AATMIN, WORK(M+N+p) )
00760  1904       CONTINUE
00761          END IF
00762 *
00763       END IF
00764 *
00765 *     For square matrix A try to determine whether A^t  would be  better
00766 *     input for the preconditioned Jacobi SVD, with faster convergence.
00767 *     The decision is based on an O(N) function of the vector of column
00768 *     and row norms of A, based on the Shannon entropy. This should give
00769 *     the right choice in most cases when the difference actually matters.
00770 *     It may fail and pick the slower converging side.
00771 *
00772       ENTRA  = ZERO
00773       ENTRAT = ZERO
00774       IF ( L2TRAN ) THEN
00775 *
00776          XSC   = ZERO
00777          TEMP1 = ONE
00778          CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
00779          TEMP1 = ONE / TEMP1
00780 *
00781          ENTRA = ZERO
00782          DO 1113 p = 1, N
00783             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
00784             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
00785  1113    CONTINUE
00786          ENTRA = - ENTRA / DLOG(DBLE(N))
00787 *
00788 *        Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex.
00789 *        It is derived from the diagonal of  A^t * A.  Do the same with the
00790 *        diagonal of A * A^t, compute the entropy of the corresponding
00791 *        probability distribution. Note that A * A^t and A^t * A have the
00792 *        same trace.
00793 *
00794          ENTRAT = ZERO
00795          DO 1114 p = N+1, N+M
00796             BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1
00797             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
00798  1114    CONTINUE
00799          ENTRAT = - ENTRAT / DLOG(DBLE(M))
00800 *
00801 *        Analyze the entropies and decide A or A^t. Smaller entropy
00802 *        usually means better input for the algorithm.
00803 *
00804          TRANSP = ( ENTRAT .LT. ENTRA )
00805 *
00806 *        If A^t is better than A, transpose A.
00807 *
00808          IF ( TRANSP ) THEN
00809 *           In an optimal implementation, this trivial transpose
00810 *           should be replaced with faster transpose.
00811             DO 1115 p = 1, N - 1
00812                DO 1116 q = p + 1, N
00813                    TEMP1 = A(q,p)
00814                   A(q,p) = A(p,q)
00815                   A(p,q) = TEMP1
00816  1116          CONTINUE
00817  1115       CONTINUE
00818             DO 1117 p = 1, N
00819                WORK(M+N+p) = SVA(p)
00820                SVA(p)      = WORK(N+p)
00821  1117       CONTINUE
00822             TEMP1  = AAPP
00823             AAPP   = AATMAX
00824             AATMAX = TEMP1
00825             TEMP1  = AAQQ
00826             AAQQ   = AATMIN
00827             AATMIN = TEMP1
00828             KILL   = LSVEC
00829             LSVEC  = RSVEC
00830             RSVEC  = KILL
00831             IF ( LSVEC ) N1 = N
00832 *
00833             ROWPIV = .TRUE.
00834          END IF
00835 *
00836       END IF
00837 *     END IF L2TRAN
00838 *
00839 *     Scale the matrix so that its maximal singular value remains less
00840 *     than DSQRT(BIG) -- the matrix is scaled so that its maximal column
00841 *     has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep
00842 *     DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and
00843 *     BLAS routines that, in some implementations, are not capable of
00844 *     working in the full interval [SFMIN,BIG] and that they may provoke
00845 *     overflows in the intermediate results. If the singular values spread
00846 *     from SFMIN to BIG, then DGESVJ will compute them. So, in that case,
00847 *     one should use DGESVJ instead of DGEJSV.
00848 *
00849       BIG1   = DSQRT( BIG )
00850       TEMP1  = DSQRT( BIG / DBLE(N) )
00851 *
00852       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
00853       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
00854           AAQQ = ( AAQQ / AAPP ) * TEMP1
00855       ELSE
00856           AAQQ = ( AAQQ * TEMP1 ) / AAPP
00857       END IF
00858       TEMP1 = TEMP1 * SCALEM
00859       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
00860 *
00861 *     To undo scaling at the end of this procedure, multiply the
00862 *     computed singular values with USCAL2 / USCAL1.
00863 *
00864       USCAL1 = TEMP1
00865       USCAL2 = AAPP
00866 *
00867       IF ( L2KILL ) THEN
00868 *        L2KILL enforces computation of nonzero singular values in
00869 *        the restricted range of condition number of the initial A,
00870 *        sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN).
00871          XSC = DSQRT( SFMIN )
00872       ELSE
00873          XSC = SMALL
00874 *
00875 *        Now, if the condition number of A is too big,
00876 *        sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN,
00877 *        as a precaution measure, the full SVD is computed using DGESVJ
00878 *        with accumulated Jacobi rotations. This provides numerically
00879 *        more robust computation, at the cost of slightly increased run
00880 *        time. Depending on the concrete implementation of BLAS and LAPACK
00881 *        (i.e. how they behave in presence of extreme ill-conditioning) the
00882 *        implementor may decide to remove this switch.
00883          IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
00884             JRACC = .TRUE.
00885          END IF
00886 *
00887       END IF
00888       IF ( AAQQ .LT. XSC ) THEN
00889          DO 700 p = 1, N
00890             IF ( SVA(p) .LT. XSC ) THEN
00891                CALL DLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA )
00892                SVA(p) = ZERO
00893             END IF
00894  700     CONTINUE
00895       END IF
00896 *
00897 *     Preconditioning using QR factorization with pivoting
00898 *
00899       IF ( ROWPIV ) THEN
00900 *        Optional row permutation (Bjoerck row pivoting):
00901 *        A result by Cox and Higham shows that the Bjoerck's
00902 *        row pivoting combined with standard column pivoting
00903 *        has similar effect as Powell-Reid complete pivoting.
00904 *        The ell-infinity norms of A are made nonincreasing.
00905          DO 1952 p = 1, M - 1
00906             q = IDAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1
00907             IWORK(2*N+p) = q
00908             IF ( p .NE. q ) THEN
00909                TEMP1       = WORK(M+N+p)
00910                WORK(M+N+p) = WORK(M+N+q)
00911                WORK(M+N+q) = TEMP1
00912             END IF
00913  1952    CONTINUE
00914          CALL DLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 )
00915       END IF
00916 *
00917 *     End of the preparation phase (scaling, optional sorting and
00918 *     transposing, optional flushing of small columns).
00919 *
00920 *     Preconditioning
00921 *
00922 *     If the full SVD is needed, the right singular vectors are computed
00923 *     from a matrix equation, and for that we need theoretical analysis
00924 *     of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF.
00925 *     In all other cases the first RR QRF can be chosen by other criteria
00926 *     (eg speed by replacing global with restricted window pivoting, such
00927 *     as in SGEQPX from TOMS # 782). Good results will be obtained using
00928 *     SGEQPX with properly (!) chosen numerical parameters.
00929 *     Any improvement of DGEQP3 improves overal performance of DGEJSV.
00930 *
00931 *     A * P1 = Q1 * [ R1^t 0]^t:
00932       DO 1963 p = 1, N
00933 *        .. all columns are free columns
00934          IWORK(p) = 0
00935  1963 CONTINUE
00936       CALL DGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR )
00937 *
00938 *     The upper triangular matrix R1 from the first QRF is inspected for
00939 *     rank deficiency and possibilities for deflation, or possible
00940 *     ill-conditioning. Depending on the user specified flag L2RANK,
00941 *     the procedure explores possibilities to reduce the numerical
00942 *     rank by inspecting the computed upper triangular factor. If
00943 *     L2RANK or L2ABER are up, then DGEJSV will compute the SVD of
00944 *     A + dA, where ||dA|| <= f(M,N)*EPSLN.
00945 *
00946       NR = 1
00947       IF ( L2ABER ) THEN
00948 *        Standard absolute error bound suffices. All sigma_i with
00949 *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
00950 *        agressive enforcement of lower numerical rank by introducing a
00951 *        backward error of the order of N*EPSLN*||A||.
00952          TEMP1 = DSQRT(DBLE(N))*EPSLN
00953          DO 3001 p = 2, N
00954             IF ( DABS(A(p,p)) .GE. (TEMP1*DABS(A(1,1))) ) THEN
00955                NR = NR + 1
00956             ELSE
00957                GO TO 3002
00958             END IF
00959  3001    CONTINUE
00960  3002    CONTINUE
00961       ELSE IF ( L2RANK ) THEN
00962 *        .. similarly as above, only slightly more gentle (less agressive).
00963 *        Sudden drop on the diagonal of R1 is used as the criterion for
00964 *        close-to-rank-defficient.
00965          TEMP1 = DSQRT(SFMIN)
00966          DO 3401 p = 2, N
00967             IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR.
00968      $           ( DABS(A(p,p)) .LT. SMALL ) .OR.
00969      $           ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
00970             NR = NR + 1
00971  3401    CONTINUE
00972  3402    CONTINUE
00973 *
00974       ELSE
00975 *        The goal is high relative accuracy. However, if the matrix
00976 *        has high scaled condition number the relative accuracy is in
00977 *        general not feasible. Later on, a condition number estimator
00978 *        will be deployed to estimate the scaled condition number.
00979 *        Here we just remove the underflowed part of the triangular
00980 *        factor. This prevents the situation in which the code is
00981 *        working hard to get the accuracy not warranted by the data.
00982          TEMP1  = DSQRT(SFMIN)
00983          DO 3301 p = 2, N
00984             IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR.
00985      $          ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
00986             NR = NR + 1
00987  3301    CONTINUE
00988  3302    CONTINUE
00989 *
00990       END IF
00991 *
00992       ALMORT = .FALSE.
00993       IF ( NR .EQ. N ) THEN
00994          MAXPRJ = ONE
00995          DO 3051 p = 2, N
00996             TEMP1  = DABS(A(p,p)) / SVA(IWORK(p))
00997             MAXPRJ = DMIN1( MAXPRJ, TEMP1 )
00998  3051    CONTINUE
00999          IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
01000       END IF
01001 *
01002 *
01003       SCONDA = - ONE
01004       CONDR1 = - ONE
01005       CONDR2 = - ONE
01006 *
01007       IF ( ERREST ) THEN
01008          IF ( N .EQ. NR ) THEN
01009             IF ( RSVEC ) THEN
01010 *              .. V is available as workspace
01011                CALL DLACPY( 'U', N, N, A, LDA, V, LDV )
01012                DO 3053 p = 1, N
01013                   TEMP1 = SVA(IWORK(p))
01014                   CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 )
01015  3053          CONTINUE
01016                CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1,
01017      $              WORK(N+1), IWORK(2*N+M+1), IERR )
01018             ELSE IF ( LSVEC ) THEN
01019 *              .. U is available as workspace
01020                CALL DLACPY( 'U', N, N, A, LDA, U, LDU )
01021                DO 3054 p = 1, N
01022                   TEMP1 = SVA(IWORK(p))
01023                   CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 )
01024  3054          CONTINUE
01025                CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1,
01026      $              WORK(N+1), IWORK(2*N+M+1), IERR )
01027             ELSE
01028                CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N )
01029                DO 3052 p = 1, N
01030                   TEMP1 = SVA(IWORK(p))
01031                   CALL DSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 )
01032  3052          CONTINUE
01033 *           .. the columns of R are scaled to have unit Euclidean lengths.
01034                CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1,
01035      $              WORK(N+N*N+1), IWORK(2*N+M+1), IERR )
01036             END IF
01037             SCONDA = ONE / DSQRT(TEMP1)
01038 *           SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
01039 *           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
01040          ELSE
01041             SCONDA = - ONE
01042          END IF
01043       END IF
01044 *
01045       L2PERT = L2PERT .AND. ( DABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) )
01046 *     If there is no violent scaling, artificial perturbation is not needed.
01047 *
01048 *     Phase 3:
01049 *
01050       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
01051 *
01052 *         Singular Values only
01053 *
01054 *         .. transpose A(1:NR,1:N)
01055          DO 1946 p = 1, MIN0( N-1, NR )
01056             CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
01057  1946    CONTINUE
01058 *
01059 *        The following two DO-loops introduce small relative perturbation
01060 *        into the strict upper triangle of the lower triangular matrix.
01061 *        Small entries below the main diagonal are also changed.
01062 *        This modification is useful if the computing environment does not
01063 *        provide/allow FLUSH TO ZERO underflow, for it prevents many
01064 *        annoying denormalized numbers in case of strongly scaled matrices.
01065 *        The perturbation is structured so that it does not introduce any
01066 *        new perturbation of the singular values, and it does not destroy
01067 *        the job done by the preconditioner.
01068 *        The licence for this perturbation is in the variable L2PERT, which
01069 *        should be .FALSE. if FLUSH TO ZERO underflow is active.
01070 *
01071          IF ( .NOT. ALMORT ) THEN
01072 *
01073             IF ( L2PERT ) THEN
01074 *              XSC = DSQRT(SMALL)
01075                XSC = EPSLN / DBLE(N)
01076                DO 4947 q = 1, NR
01077                   TEMP1 = XSC*DABS(A(q,q))
01078                   DO 4949 p = 1, N
01079                      IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
01080      $                    .OR. ( p .LT. q ) )
01081      $                     A(p,q) = DSIGN( TEMP1, A(p,q) )
01082  4949             CONTINUE
01083  4947          CONTINUE
01084             ELSE
01085                CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA )
01086             END IF
01087 *
01088 *            .. second preconditioning using the QR factorization
01089 *
01090             CALL DGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR )
01091 *
01092 *           .. and transpose upper to lower triangular
01093             DO 1948 p = 1, NR - 1
01094                CALL DCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
01095  1948       CONTINUE
01096 *
01097          END IF
01098 *
01099 *           Row-cyclic Jacobi SVD algorithm with column pivoting
01100 *
01101 *           .. again some perturbation (a "background noise") is added
01102 *           to drown denormals
01103             IF ( L2PERT ) THEN
01104 *              XSC = DSQRT(SMALL)
01105                XSC = EPSLN / DBLE(N)
01106                DO 1947 q = 1, NR
01107                   TEMP1 = XSC*DABS(A(q,q))
01108                   DO 1949 p = 1, NR
01109                      IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
01110      $                       .OR. ( p .LT. q ) )
01111      $                   A(p,q) = DSIGN( TEMP1, A(p,q) )
01112  1949             CONTINUE
01113  1947          CONTINUE
01114             ELSE
01115                CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
01116             END IF
01117 *
01118 *           .. and one-sided Jacobi rotations are started on a lower
01119 *           triangular matrix (plus perturbation which is ignored in
01120 *           the part which destroys triangular form (confusing?!))
01121 *
01122             CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA,
01123      $                      N, V, LDV, WORK, LWORK, INFO )
01124 *
01125             SCALEM  = WORK(1)
01126             NUMRANK = IDNINT(WORK(2))
01127 *
01128 *
01129       ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
01130 *
01131 *        -> Singular Values and Right Singular Vectors <-
01132 *
01133          IF ( ALMORT ) THEN
01134 *
01135 *           .. in this case NR equals N
01136             DO 1998 p = 1, NR
01137                CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
01138  1998       CONTINUE
01139             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01140 *
01141             CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
01142      $                  WORK, LWORK, INFO )
01143             SCALEM  = WORK(1)
01144             NUMRANK = IDNINT(WORK(2))
01145 
01146          ELSE
01147 *
01148 *        .. two more QR factorizations ( one QRF is not enough, two require
01149 *        accumulated product of Jacobi rotations, three are perfect )
01150 *
01151             CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
01152             CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
01153             CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
01154             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01155             CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01156      $                   LWORK-2*N, IERR )
01157             DO 8998 p = 1, NR
01158                CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
01159  8998       CONTINUE
01160             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01161 *
01162             CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
01163      $                  LDU, WORK(N+1), LWORK, INFO )
01164             SCALEM  = WORK(N+1)
01165             NUMRANK = IDNINT(WORK(N+2))
01166             IF ( NR .LT. N ) THEN
01167                CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),   LDV )
01168                CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),   LDV )
01169                CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
01170             END IF
01171 *
01172          CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
01173      $               V, LDV, WORK(N+1), LWORK-N, IERR )
01174 *
01175          END IF
01176 *
01177          DO 8991 p = 1, N
01178             CALL DCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
01179  8991    CONTINUE
01180          CALL DLACPY( 'All', N, N, A, LDA, V, LDV )
01181 *
01182          IF ( TRANSP ) THEN
01183             CALL DLACPY( 'All', N, N, V, LDV, U, LDU )
01184          END IF
01185 *
01186       ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
01187 *
01188 *        .. Singular Values and Left Singular Vectors                 ..
01189 *
01190 *        .. second preconditioning step to avoid need to accumulate
01191 *        Jacobi rotations in the Jacobi iterations.
01192          DO 1965 p = 1, NR
01193             CALL DCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
01194  1965    CONTINUE
01195          CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
01196 *
01197          CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1),
01198      $              LWORK-2*N, IERR )
01199 *
01200          DO 1967 p = 1, NR - 1
01201             CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
01202  1967    CONTINUE
01203          CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
01204 *
01205          CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
01206      $        LDA, WORK(N+1), LWORK-N, INFO )
01207          SCALEM  = WORK(N+1)
01208          NUMRANK = IDNINT(WORK(N+2))
01209 *
01210          IF ( NR .LT. M ) THEN
01211             CALL DLASET( 'A',  M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
01212             IF ( NR .LT. N1 ) THEN
01213                CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
01214                CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
01215             END IF
01216          END IF
01217 *
01218          CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
01219      $               LDU, WORK(N+1), LWORK-N, IERR )
01220 *
01221          IF ( ROWPIV )
01222      $       CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01223 *
01224          DO 1974 p = 1, N1
01225             XSC = ONE / DNRM2( M, U(1,p), 1 )
01226             CALL DSCAL( M, XSC, U(1,p), 1 )
01227  1974    CONTINUE
01228 *
01229          IF ( TRANSP ) THEN
01230             CALL DLACPY( 'All', N, N, U, LDU, V, LDV )
01231          END IF
01232 *
01233       ELSE
01234 *
01235 *        .. Full SVD ..
01236 *
01237          IF ( .NOT. JRACC ) THEN
01238 *
01239          IF ( .NOT. ALMORT ) THEN
01240 *
01241 *           Second Preconditioning Step (QRF [with pivoting])
01242 *           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
01243 *           equivalent to an LQF CALL. Since in many libraries the QRF
01244 *           seems to be better optimized than the LQF, we do explicit
01245 *           transpose and use the QRF. This is subject to changes in an
01246 *           optimized implementation of DGEJSV.
01247 *
01248             DO 1968 p = 1, NR
01249                CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
01250  1968       CONTINUE
01251 *
01252 *           .. the following two loops perturb small entries to avoid
01253 *           denormals in the second QR factorization, where they are
01254 *           as good as zeros. This is done to avoid painfully slow
01255 *           computation with denormals. The relative size of the perturbation
01256 *           is a parameter that can be changed by the implementer.
01257 *           This perturbation device will be obsolete on machines with
01258 *           properly implemented arithmetic.
01259 *           To switch it off, set L2PERT=.FALSE. To remove it from  the
01260 *           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
01261 *           The following two loops should be blocked and fused with the
01262 *           transposed copy above.
01263 *
01264             IF ( L2PERT ) THEN
01265                XSC = DSQRT(SMALL)
01266                DO 2969 q = 1, NR
01267                   TEMP1 = XSC*DABS( V(q,q) )
01268                   DO 2968 p = 1, N
01269                      IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
01270      $                   .OR. ( p .LT. q ) )
01271      $                   V(p,q) = DSIGN( TEMP1, V(p,q) )
01272                      IF ( p .LT. q ) V(p,q) = - V(p,q)
01273  2968             CONTINUE
01274  2969          CONTINUE
01275             ELSE
01276                CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01277             END IF
01278 *
01279 *           Estimate the row scaled condition number of R1
01280 *           (If R1 is rectangular, N > NR, then the condition number
01281 *           of the leading NR x NR submatrix is estimated.)
01282 *
01283             CALL DLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR )
01284             DO 3950 p = 1, NR
01285                TEMP1 = DNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1)
01286                CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1)
01287  3950       CONTINUE
01288             CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
01289      $                   WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
01290             CONDR1 = ONE / DSQRT(TEMP1)
01291 *           .. here need a second oppinion on the condition number
01292 *           .. then assume worst case scenario
01293 *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
01294 *           more conservative    <=> CONDR1 .LT. DSQRT(DBLE(N))
01295 *
01296             COND_OK = DSQRT(DBLE(NR))
01297 *[TP]       COND_OK is a tuning parameter.
01298 
01299             IF ( CONDR1 .LT. COND_OK ) THEN
01300 *              .. the second QRF without pivoting. Note: in an optimized
01301 *              implementation, this QRF should be implemented as the QRF
01302 *              of a lower triangular matrix.
01303 *              R1^t = Q2 * R2
01304                CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01305      $              LWORK-2*N, IERR )
01306 *
01307                IF ( L2PERT ) THEN
01308                   XSC = DSQRT(SMALL)/EPSLN
01309                   DO 3959 p = 2, NR
01310                      DO 3958 q = 1, p - 1
01311                         TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
01312                         IF ( DABS(V(q,p)) .LE. TEMP1 )
01313      $                     V(q,p) = DSIGN( TEMP1, V(q,p) )
01314  3958                CONTINUE
01315  3959             CONTINUE
01316                END IF
01317 *
01318                IF ( NR .NE. N )
01319      $         CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
01320 *              .. save ...
01321 *
01322 *           .. this transposed copy should be better than naive
01323                DO 1969 p = 1, NR - 1
01324                   CALL DCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
01325  1969          CONTINUE
01326 *
01327                CONDR2 = CONDR1
01328 *
01329             ELSE
01330 *
01331 *              .. ill-conditioned case: second QRF with pivoting
01332 *              Note that windowed pivoting would be equaly good
01333 *              numerically, and more run-time efficient. So, in
01334 *              an optimal implementation, the next call to DGEQP3
01335 *              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
01336 *              with properly (carefully) chosen parameters.
01337 *
01338 *              R1^t * P2 = Q2 * R2
01339                DO 3003 p = 1, NR
01340                   IWORK(N+p) = 0
01341  3003          CONTINUE
01342                CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1),
01343      $                  WORK(2*N+1), LWORK-2*N, IERR )
01344 **               CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01345 **     $              LWORK-2*N, IERR )
01346                IF ( L2PERT ) THEN
01347                   XSC = DSQRT(SMALL)
01348                   DO 3969 p = 2, NR
01349                      DO 3968 q = 1, p - 1
01350                         TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
01351                         IF ( DABS(V(q,p)) .LE. TEMP1 )
01352      $                     V(q,p) = DSIGN( TEMP1, V(q,p) )
01353  3968                CONTINUE
01354  3969             CONTINUE
01355                END IF
01356 *
01357                CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
01358 *
01359                IF ( L2PERT ) THEN
01360                   XSC = DSQRT(SMALL)
01361                   DO 8970 p = 2, NR
01362                      DO 8971 q = 1, p - 1
01363                         TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
01364                         V(p,q) = - DSIGN( TEMP1, V(q,p) )
01365  8971                CONTINUE
01366  8970             CONTINUE
01367                ELSE
01368                   CALL DLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV )
01369                END IF
01370 *              Now, compute R2 = L3 * Q3, the LQ factorization.
01371                CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1),
01372      $               WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
01373 *              .. and estimate the condition number
01374                CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR )
01375                DO 4950 p = 1, NR
01376                   TEMP1 = DNRM2( p, WORK(2*N+N*NR+NR+p), NR )
01377                   CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR )
01378  4950          CONTINUE
01379                CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
01380      $              WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR )
01381                CONDR2 = ONE / DSQRT(TEMP1)
01382 *
01383                IF ( CONDR2 .GE. COND_OK ) THEN
01384 *                 .. save the Householder vectors used for Q3
01385 *                 (this overwrittes the copy of R2, as it will not be
01386 *                 needed in this branch, but it does not overwritte the
01387 *                 Huseholder vectors of Q2.).
01388                   CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
01389 *                 .. and the rest of the information on Q3 is in
01390 *                 WORK(2*N+N*NR+1:2*N+N*NR+N)
01391                END IF
01392 *
01393             END IF
01394 *
01395             IF ( L2PERT ) THEN
01396                XSC = DSQRT(SMALL)
01397                DO 4968 q = 2, NR
01398                   TEMP1 = XSC * V(q,q)
01399                   DO 4969 p = 1, q - 1
01400 *                    V(p,q) = - DSIGN( TEMP1, V(q,p) )
01401                      V(p,q) = - DSIGN( TEMP1, V(p,q) )
01402  4969             CONTINUE
01403  4968          CONTINUE
01404             ELSE
01405                CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
01406             END IF
01407 *
01408 *        Second preconditioning finished; continue with Jacobi SVD
01409 *        The input matrix is lower trinagular.
01410 *
01411 *        Recover the right singular vectors as solution of a well
01412 *        conditioned triangular matrix equation.
01413 *
01414             IF ( CONDR1 .LT. COND_OK ) THEN
01415 *
01416                CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U,
01417      $              LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO )
01418                SCALEM  = WORK(2*N+N*NR+NR+1)
01419                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
01420                DO 3970 p = 1, NR
01421                   CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
01422                   CALL DSCAL( NR, SVA(p),    V(1,p), 1 )
01423  3970          CONTINUE
01424 
01425 *        .. pick the right matrix equation and solve it
01426 *
01427                IF ( NR .EQ. N ) THEN
01428 * :))             .. best case, R1 is inverted. The solution of this matrix
01429 *                 equation is Q2*V2 = the product of the Jacobi rotations
01430 *                 used in DGESVJ, premultiplied with the orthogonal matrix
01431 *                 from the second QR factorization.
01432                   CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
01433                ELSE
01434 *                 .. R1 is well conditioned, but non-square. Transpose(R2)
01435 *                 is inverted to get the product of the Jacobi rotations
01436 *                 used in DGESVJ. The Q-factor from the second QR
01437 *                 factorization is then built in explicitly.
01438                   CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1),
01439      $                 N,V,LDV)
01440                   IF ( NR .LT. N ) THEN
01441                     CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
01442                     CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
01443                     CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
01444                   END IF
01445                   CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01446      $                 V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
01447                END IF
01448 *
01449             ELSE IF ( CONDR2 .LT. COND_OK ) THEN
01450 *
01451 * :)           .. the input matrix A is very likely a relative of
01452 *              the Kahan matrix :)
01453 *              The matrix R2 is inverted. The solution of the matrix equation
01454 *              is Q3^T*V3 = the product of the Jacobi rotations (appplied to
01455 *              the lower triangular L3 from the LQ factorization of
01456 *              R2=L3*Q3), pre-multiplied with the transposed Q3.
01457                CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
01458      $              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
01459                SCALEM  = WORK(2*N+N*NR+NR+1)
01460                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
01461                DO 3870 p = 1, NR
01462                   CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
01463                   CALL DSCAL( NR, SVA(p),    U(1,p), 1 )
01464  3870          CONTINUE
01465                CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
01466 *              .. apply the permutation from the second QR factorization
01467                DO 873 q = 1, NR
01468                   DO 872 p = 1, NR
01469                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
01470  872              CONTINUE
01471                   DO 874 p = 1, NR
01472                      U(p,q) = WORK(2*N+N*NR+NR+p)
01473  874              CONTINUE
01474  873           CONTINUE
01475                IF ( NR .LT. N ) THEN
01476                   CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
01477                   CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
01478                   CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
01479                END IF
01480                CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01481      $              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
01482             ELSE
01483 *              Last line of defense.
01484 * #:(          This is a rather pathological case: no scaled condition
01485 *              improvement after two pivoted QR factorizations. Other
01486 *              possibility is that the rank revealing QR factorization
01487 *              or the condition estimator has failed, or the COND_OK
01488 *              is set very close to ONE (which is unnecessary). Normally,
01489 *              this branch should never be executed, but in rare cases of
01490 *              failure of the RRQR or condition estimator, the last line of
01491 *              defense ensures that DGEJSV completes the task.
01492 *              Compute the full SVD of L3 using DGESVJ with explicit
01493 *              accumulation of Jacobi rotations.
01494                CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
01495      $              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
01496                SCALEM  = WORK(2*N+N*NR+NR+1)
01497                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
01498                IF ( NR .LT. N ) THEN
01499                   CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
01500                   CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
01501                   CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
01502                END IF
01503                CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01504      $              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
01505 *
01506                CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N,
01507      $              WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1),
01508      $              LWORK-2*N-N*NR-NR, IERR )
01509                DO 773 q = 1, NR
01510                   DO 772 p = 1, NR
01511                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
01512  772              CONTINUE
01513                   DO 774 p = 1, NR
01514                      U(p,q) = WORK(2*N+N*NR+NR+p)
01515  774              CONTINUE
01516  773           CONTINUE
01517 *
01518             END IF
01519 *
01520 *           Permute the rows of V using the (column) permutation from the
01521 *           first QRF. Also, scale the columns to make them unit in
01522 *           Euclidean norm. This applies to all cases.
01523 *
01524             TEMP1 = DSQRT(DBLE(N)) * EPSLN
01525             DO 1972 q = 1, N
01526                DO 972 p = 1, N
01527                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
01528   972          CONTINUE
01529                DO 973 p = 1, N
01530                   V(p,q) = WORK(2*N+N*NR+NR+p)
01531   973          CONTINUE
01532                XSC = ONE / DNRM2( N, V(1,q), 1 )
01533                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01534      $           CALL DSCAL( N, XSC, V(1,q), 1 )
01535  1972       CONTINUE
01536 *           At this moment, V contains the right singular vectors of A.
01537 *           Next, assemble the left singular vector matrix U (M x N).
01538             IF ( NR .LT. M ) THEN
01539                CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
01540                IF ( NR .LT. N1 ) THEN
01541                   CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
01542                   CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
01543                END IF
01544             END IF
01545 *
01546 *           The Q matrix from the first QRF is built into the left singular
01547 *           matrix U. This applies to all cases.
01548 *
01549             CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U,
01550      $           LDU, WORK(N+1), LWORK-N, IERR )
01551 
01552 *           The columns of U are normalized. The cost is O(M*N) flops.
01553             TEMP1 = DSQRT(DBLE(M)) * EPSLN
01554             DO 1973 p = 1, NR
01555                XSC = ONE / DNRM2( M, U(1,p), 1 )
01556                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01557      $          CALL DSCAL( M, XSC, U(1,p), 1 )
01558  1973       CONTINUE
01559 *
01560 *           If the initial QRF is computed with row pivoting, the left
01561 *           singular vectors must be adjusted.
01562 *
01563             IF ( ROWPIV )
01564      $          CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01565 *
01566          ELSE
01567 *
01568 *        .. the initial matrix A has almost orthogonal columns and
01569 *        the second QRF is not needed
01570 *
01571             CALL DLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N )
01572             IF ( L2PERT ) THEN
01573                XSC = DSQRT(SMALL)
01574                DO 5970 p = 2, N
01575                   TEMP1 = XSC * WORK( N + (p-1)*N + p )
01576                   DO 5971 q = 1, p - 1
01577                      WORK(N+(q-1)*N+p)=-DSIGN(TEMP1,WORK(N+(p-1)*N+q))
01578  5971             CONTINUE
01579  5970          CONTINUE
01580             ELSE
01581                CALL DLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N )
01582             END IF
01583 *
01584             CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA,
01585      $           N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO )
01586 *
01587             SCALEM  = WORK(N+N*N+1)
01588             NUMRANK = IDNINT(WORK(N+N*N+2))
01589             DO 6970 p = 1, N
01590                CALL DCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 )
01591                CALL DSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 )
01592  6970       CONTINUE
01593 *
01594             CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N,
01595      $           ONE, A, LDA, WORK(N+1), N )
01596             DO 6972 p = 1, N
01597                CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV )
01598  6972       CONTINUE
01599             TEMP1 = DSQRT(DBLE(N))*EPSLN
01600             DO 6971 p = 1, N
01601                XSC = ONE / DNRM2( N, V(1,p), 1 )
01602                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01603      $            CALL DSCAL( N, XSC, V(1,p), 1 )
01604  6971       CONTINUE
01605 *
01606 *           Assemble the left singular vector matrix U (M x N).
01607 *
01608             IF ( N .LT. M ) THEN
01609                CALL DLASET( 'A',  M-N, N, ZERO, ZERO, U(N+1,1), LDU )
01610                IF ( N .LT. N1 ) THEN
01611                   CALL DLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
01612                   CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
01613                END IF
01614             END IF
01615             CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
01616      $           LDU, WORK(N+1), LWORK-N, IERR )
01617             TEMP1 = DSQRT(DBLE(M))*EPSLN
01618             DO 6973 p = 1, N1
01619                XSC = ONE / DNRM2( M, U(1,p), 1 )
01620                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01621      $            CALL DSCAL( M, XSC, U(1,p), 1 )
01622  6973       CONTINUE
01623 *
01624             IF ( ROWPIV )
01625      $         CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01626 *
01627          END IF
01628 *
01629 *        end of the  >> almost orthogonal case <<  in the full SVD
01630 *
01631          ELSE
01632 *
01633 *        This branch deploys a preconditioned Jacobi SVD with explicitly
01634 *        accumulated rotations. It is included as optional, mainly for
01635 *        experimental purposes. It does perfom well, and can also be used.
01636 *        In this implementation, this branch will be automatically activated
01637 *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
01638 *        to be greater than the overflow threshold. This is because the
01639 *        a posteriori computation of the singular vectors assumes robust
01640 *        implementation of BLAS and some LAPACK procedures, capable of working
01641 *        in presence of extreme values. Since that is not always the case, ...
01642 *
01643          DO 7968 p = 1, NR
01644             CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
01645  7968    CONTINUE
01646 *
01647          IF ( L2PERT ) THEN
01648             XSC = DSQRT(SMALL/EPSLN)
01649             DO 5969 q = 1, NR
01650                TEMP1 = XSC*DABS( V(q,q) )
01651                DO 5968 p = 1, N
01652                   IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
01653      $                .OR. ( p .LT. q ) )
01654      $                V(p,q) = DSIGN( TEMP1, V(p,q) )
01655                   IF ( p .LT. q ) V(p,q) = - V(p,q)
01656  5968          CONTINUE
01657  5969       CONTINUE
01658          ELSE
01659             CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
01660          END IF
01661 
01662          CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
01663      $        LWORK-2*N, IERR )
01664          CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N )
01665 *
01666          DO 7969 p = 1, NR
01667             CALL DCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
01668  7969    CONTINUE
01669 
01670          IF ( L2PERT ) THEN
01671             XSC = DSQRT(SMALL/EPSLN)
01672             DO 9970 q = 2, NR
01673                DO 9971 p = 1, q - 1
01674                   TEMP1 = XSC * DMIN1(DABS(U(p,p)),DABS(U(q,q)))
01675                   U(p,q) = - DSIGN( TEMP1, U(q,p) )
01676  9971          CONTINUE
01677  9970       CONTINUE
01678          ELSE
01679             CALL DLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
01680          END IF
01681 
01682          CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA,
01683      $        N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO )
01684          SCALEM  = WORK(2*N+N*NR+1)
01685          NUMRANK = IDNINT(WORK(2*N+N*NR+2))
01686 
01687          IF ( NR .LT. N ) THEN
01688             CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
01689             CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
01690             CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
01691          END IF
01692 
01693          CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
01694      $        V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
01695 *
01696 *           Permute the rows of V using the (column) permutation from the
01697 *           first QRF. Also, scale the columns to make them unit in
01698 *           Euclidean norm. This applies to all cases.
01699 *
01700             TEMP1 = DSQRT(DBLE(N)) * EPSLN
01701             DO 7972 q = 1, N
01702                DO 8972 p = 1, N
01703                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
01704  8972          CONTINUE
01705                DO 8973 p = 1, N
01706                   V(p,q) = WORK(2*N+N*NR+NR+p)
01707  8973          CONTINUE
01708                XSC = ONE / DNRM2( N, V(1,q), 1 )
01709                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
01710      $           CALL DSCAL( N, XSC, V(1,q), 1 )
01711  7972       CONTINUE
01712 *
01713 *           At this moment, V contains the right singular vectors of A.
01714 *           Next, assemble the left singular vector matrix U (M x N).
01715 *
01716          IF ( NR .LT. M ) THEN
01717             CALL DLASET( 'A',  M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
01718             IF ( NR .LT. N1 ) THEN
01719                CALL DLASET( 'A',NR,  N1-NR, ZERO, ZERO,  U(1,NR+1),LDU )
01720                CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
01721             END IF
01722          END IF
01723 *
01724          CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
01725      $        LDU, WORK(N+1), LWORK-N, IERR )
01726 *
01727             IF ( ROWPIV )
01728      $         CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
01729 *
01730 *
01731          END IF
01732          IF ( TRANSP ) THEN
01733 *           .. swap U and V because the procedure worked on A^t
01734             DO 6974 p = 1, N
01735                CALL DSWAP( N, U(1,p), 1, V(1,p), 1 )
01736  6974       CONTINUE
01737          END IF
01738 *
01739       END IF
01740 *     end of the full SVD
01741 *
01742 *     Undo scaling, if necessary (and possible)
01743 *
01744       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
01745          CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
01746          USCAL1 = ONE
01747          USCAL2 = ONE
01748       END IF
01749 *
01750       IF ( NR .LT. N ) THEN
01751          DO 3004 p = NR+1, N
01752             SVA(p) = ZERO
01753  3004    CONTINUE
01754       END IF
01755 *
01756       WORK(1) = USCAL2 * SCALEM
01757       WORK(2) = USCAL1
01758       IF ( ERREST ) WORK(3) = SCONDA
01759       IF ( LSVEC .AND. RSVEC ) THEN
01760          WORK(4) = CONDR1
01761          WORK(5) = CONDR2
01762       END IF
01763       IF ( L2TRAN ) THEN
01764          WORK(6) = ENTRA
01765          WORK(7) = ENTRAT
01766       END IF
01767 *
01768       IWORK(1) = NR
01769       IWORK(2) = NUMRANK
01770       IWORK(3) = WARNING
01771 *
01772       RETURN
01773 *     ..
01774 *     .. END OF DGEJSV
01775 *     ..
01776       END
01777 *
 All Files Functions