SRC/dgejsv.f

   1 *> \brief \b DGEJSV
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGEJSV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
  22 *                          M, N, A, LDA, SVA, U, LDU, V, LDV,
  23 *                          WORK, LWORK, IWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       IMPLICIT    NONE
  27 *       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
  31 *      $            WORK( LWORK )
  32 *       INTEGER     IWORK( * )
  33 *       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> DGEJSV computes the singular value decomposition (SVD) of a real M-by-N
  43 *> matrix [A], where M >= N. The SVD of [A] is written as
  44 *>
  45 *>              [A] = [U] * [SIGMA] * [V]^t,
  46 *>
  47 *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
  48 *> diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and
  49 *> [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
  50 *> the singular values of [A]. The columns of [U] and [V] are the left and
  51 *> the right singular vectors of [A], respectively. The matrices [U] and [V]
  52 *> are computed and stored in the arrays U and V, respectively. The diagonal
  53 *> of [SIGMA] is computed and stored in the array SVA.
  54 *> DGEJSV can sometimes compute tiny singular values and their singular vectors much
  55 *> more accurately than other SVD routines, see below under Further Details.
  56 *> \endverbatim
  57 *
  58 *  Arguments:
  59 *  ==========
  60 *
  61 *> \param[in] JOBA
  62 *> \verbatim
  63 *>          JOBA is CHARACTER*1
  64 *>        Specifies the level of accuracy:
  65 *>       = 'C': This option works well (high relative accuracy) if A = B * D,
  66 *>             with well-conditioned B and arbitrary diagonal matrix D.
  67 *>             The accuracy cannot be spoiled by COLUMN scaling. The
  68 *>             accuracy of the computed output depends on the condition of
  69 *>             B, and the procedure aims at the best theoretical accuracy.
  70 *>             The relative error max_{i=1:N}|d sigma_i| / sigma_i is
  71 *>             bounded by f(M,N)*epsilon* cond(B), independent of D.
  72 *>             The input matrix is preprocessed with the QRF with column
  73 *>             pivoting. This initial preprocessing and preconditioning by
  74 *>             a rank revealing QR factorization is common for all values of
  75 *>             JOBA. Additional actions are specified as follows:
  76 *>       = 'E': Computation as with 'C' with an additional estimate of the
  77 *>             condition number of B. It provides a realistic error bound.
  78 *>       = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
  79 *>             D1, D2, and well-conditioned matrix C, this option gives
  80 *>             higher accuracy than the 'C' option. If the structure of the
  81 *>             input matrix is not known, and relative accuracy is
  82 *>             desirable, then this option is advisable. The input matrix A
  83 *>             is preprocessed with QR factorization with FULL (row and
  84 *>             column) pivoting.
  85 *>       = 'G'  Computation as with 'F' with an additional estimate of the
  86 *>             condition number of B, where A=D*B. If A has heavily weighted
  87 *>             rows, then using this condition number gives too pessimistic
  88 *>             error bound.
  89 *>       = 'A': Small singular values are the noise and the matrix is treated
  90 *>             as numerically rank defficient. The error in the computed
  91 *>             singular values is bounded by f(m,n)*epsilon*||A||.
  92 *>             The computed SVD A = U * S * V^t restores A up to
  93 *>             f(m,n)*epsilon*||A||.
  94 *>             This gives the procedure the licence to discard (set to zero)
  95 *>             all singular values below N*epsilon*||A||.
  96 *>       = 'R': Similar as in 'A'. Rank revealing property of the initial
  97 *>             QR factorization is used do reveal (using triangular factor)
  98 *>             a gap sigma_{r+1} < epsilon * sigma_r in which case the
  99 *>             numerical RANK is declared to be r. The SVD is computed with
 100 *>             absolute error bounds, but more accurately than with 'A'.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] JOBU
 104 *> \verbatim
 105 *>          JOBU is CHARACTER*1
 106 *>        Specifies whether to compute the columns of U:
 107 *>       = 'U': N columns of U are returned in the array U.
 108 *>       = 'F': full set of M left sing. vectors is returned in the array U.
 109 *>       = 'W': U may be used as workspace of length M*N. See the description
 110 *>             of U.
 111 *>       = 'N': U is not computed.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] JOBV
 115 *> \verbatim
 116 *>          JOBV is CHARACTER*1
 117 *>        Specifies whether to compute the matrix V:
 118 *>       = 'V': N columns of V are returned in the array V; Jacobi rotations
 119 *>             are not explicitly accumulated.
 120 *>       = 'J': N columns of V are returned in the array V, but they are
 121 *>             computed as the product of Jacobi rotations. This option is
 122 *>             allowed only if JOBU .NE. 'N', i.e. in computing the full SVD.
 123 *>       = 'W': V may be used as workspace of length N*N. See the description
 124 *>             of V.
 125 *>       = 'N': V is not computed.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] JOBR
 129 *> \verbatim
 130 *>          JOBR is CHARACTER*1
 131 *>        Specifies the RANGE for the singular values. Issues the licence to
 132 *>        set to zero small positive singular values if they are outside
 133 *>        specified range. If A .NE. 0 is scaled so that the largest singular
 134 *>        value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
 135 *>        the licence to kill columns of A whose norm in c*A is less than
 136 *>        DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
 137 *>        where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
 138 *>       = 'N': Do not kill small columns of c*A. This option assumes that
 139 *>             BLAS and QR factorizations and triangular solvers are
 140 *>             implemented to work in that range. If the condition of A
 141 *>             is greater than BIG, use DGESVJ.
 142 *>       = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)]
 143 *>             (roughly, as described above). This option is recommended.
 144 *>                                            ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 145 *>        For computing the singular values in the FULL range [SFMIN,BIG]
 146 *>        use DGESVJ.
 147 *> \endverbatim
 148 *>
 149 *> \param[in] JOBT
 150 *> \verbatim
 151 *>          JOBT is CHARACTER*1
 152 *>        If the matrix is square then the procedure may determine to use
 153 *>        transposed A if A^t seems to be better with respect to convergence.
 154 *>        If the matrix is not square, JOBT is ignored. This is subject to
 155 *>        changes in the future.
 156 *>        The decision is based on two values of entropy over the adjoint
 157 *>        orbit of A^t * A. See the descriptions of WORK(6) and WORK(7).
 158 *>       = 'T': transpose if entropy test indicates possibly faster
 159 *>        convergence of Jacobi process if A^t is taken as input. If A is
 160 *>        replaced with A^t, then the row pivoting is included automatically.
 161 *>       = 'N': do not speculate.
 162 *>        This option can be used to compute only the singular values, or the
 163 *>        full SVD (U, SIGMA and V). For only one set of singular vectors
 164 *>        (U or V), the caller should provide both U and V, as one of the
 165 *>        matrices is used as workspace if the matrix A is transposed.
 166 *>        The implementer can easily remove this constraint and make the
 167 *>        code more complicated. See the descriptions of U and V.
 168 *> \endverbatim
 169 *>
 170 *> \param[in] JOBP
 171 *> \verbatim
 172 *>          JOBP is CHARACTER*1
 173 *>        Issues the licence to introduce structured perturbations to drown
 174 *>        denormalized numbers. This licence should be active if the
 175 *>        denormals are poorly implemented, causing slow computation,
 176 *>        especially in cases of fast convergence (!). For details see [1,2].
 177 *>        For the sake of simplicity, this perturbations are included only
 178 *>        when the full SVD or only the singular values are requested. The
 179 *>        implementer/user can easily add the perturbation for the cases of
 180 *>        computing one set of singular vectors.
 181 *>       = 'P': introduce perturbation
 182 *>       = 'N': do not perturb
 183 *> \endverbatim
 184 *>
 185 *> \param[in] M
 186 *> \verbatim
 187 *>          M is INTEGER
 188 *>         The number of rows of the input matrix A.  M >= 0.
 189 *> \endverbatim
 190 *>
 191 *> \param[in] N
 192 *> \verbatim
 193 *>          N is INTEGER
 194 *>         The number of columns of the input matrix A. M >= N >= 0.
 195 *> \endverbatim
 196 *>
 197 *> \param[in,out] A
 198 *> \verbatim
 199 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 200 *>          On entry, the M-by-N matrix A.
 201 *> \endverbatim
 202 *>
 203 *> \param[in] LDA
 204 *> \verbatim
 205 *>          LDA is INTEGER
 206 *>          The leading dimension of the array A.  LDA >= max(1,M).
 207 *> \endverbatim
 208 *>
 209 *> \param[out] SVA
 210 *> \verbatim
 211 *>          SVA is DOUBLE PRECISION array, dimension (N)
 212 *>          On exit,
 213 *>          - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
 214 *>            computation SVA contains Euclidean column norms of the
 215 *>            iterated matrices in the array A.
 216 *>          - For WORK(1) .NE. WORK(2): The singular values of A are
 217 *>            (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
 218 *>            sigma_max(A) overflows or if small singular values have been
 219 *>            saved from underflow by scaling the input matrix A.
 220 *>          - If JOBR='R' then some of the singular values may be returned
 221 *>            as exact zeros obtained by "set to zero" because they are
 222 *>            below the numerical rank threshold or are denormalized numbers.
 223 *> \endverbatim
 224 *>
 225 *> \param[out] U
 226 *> \verbatim
 227 *>          U is DOUBLE PRECISION array, dimension ( LDU, N )
 228 *>          If JOBU = 'U', then U contains on exit the M-by-N matrix of
 229 *>                         the left singular vectors.
 230 *>          If JOBU = 'F', then U contains on exit the M-by-M matrix of
 231 *>                         the left singular vectors, including an ONB
 232 *>                         of the orthogonal complement of the Range(A).
 233 *>          If JOBU = 'W'  .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
 234 *>                         then U is used as workspace if the procedure
 235 *>                         replaces A with A^t. In that case, [V] is computed
 236 *>                         in U as left singular vectors of A^t and then
 237 *>                         copied back to the V array. This 'W' option is just
 238 *>                         a reminder to the caller that in this case U is
 239 *>                         reserved as workspace of length N*N.
 240 *>          If JOBU = 'N'  U is not referenced, unless JOBT='T'.
 241 *> \endverbatim
 242 *>
 243 *> \param[in] LDU
 244 *> \verbatim
 245 *>          LDU is INTEGER
 246 *>          The leading dimension of the array U,  LDU >= 1.
 247 *>          IF  JOBU = 'U' or 'F' or 'W',  then LDU >= M.
 248 *> \endverbatim
 249 *>
 250 *> \param[out] V
 251 *> \verbatim
 252 *>          V is DOUBLE PRECISION array, dimension ( LDV, N )
 253 *>          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
 254 *>                         the right singular vectors;
 255 *>          If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
 256 *>                         then V is used as workspace if the pprocedure
 257 *>                         replaces A with A^t. In that case, [U] is computed
 258 *>                         in V as right singular vectors of A^t and then
 259 *>                         copied back to the U array. This 'W' option is just
 260 *>                         a reminder to the caller that in this case V is
 261 *>                         reserved as workspace of length N*N.
 262 *>          If JOBV = 'N'  V is not referenced, unless JOBT='T'.
 263 *> \endverbatim
 264 *>
 265 *> \param[in] LDV
 266 *> \verbatim
 267 *>          LDV is INTEGER
 268 *>          The leading dimension of the array V,  LDV >= 1.
 269 *>          If JOBV = 'V' or 'J' or 'W', then LDV >= N.
 270 *> \endverbatim
 271 *>
 272 *> \param[out] WORK
 273 *> \verbatim
 274 *>          WORK is DOUBLE PRECISION array, dimension at least LWORK.
 275 *>          On exit, if N.GT.0 .AND. M.GT.0 (else not referenced),
 276 *>          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
 277 *>                    that SCALE*SVA(1:N) are the computed singular values
 278 *>                    of A. (See the description of SVA().)
 279 *>          WORK(2) = See the description of WORK(1).
 280 *>          WORK(3) = SCONDA is an estimate for the condition number of
 281 *>                    column equilibrated A. (If JOBA .EQ. 'E' or 'G')
 282 *>                    SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
 283 *>                    It is computed using DPOCON. It holds
 284 *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
 285 *>                    where R is the triangular factor from the QRF of A.
 286 *>                    However, if R is truncated and the numerical rank is
 287 *>                    determined to be strictly smaller than N, SCONDA is
 288 *>                    returned as -1, thus indicating that the smallest
 289 *>                    singular values might be lost.
 290 *>
 291 *>          If full SVD is needed, the following two condition numbers are
 292 *>          useful for the analysis of the algorithm. They are provied for
 293 *>          a developer/implementer who is familiar with the details of
 294 *>          the method.
 295 *>
 296 *>          WORK(4) = an estimate of the scaled condition number of the
 297 *>                    triangular factor in the first QR factorization.
 298 *>          WORK(5) = an estimate of the scaled condition number of the
 299 *>                    triangular factor in the second QR factorization.
 300 *>          The following two parameters are computed if JOBT .EQ. 'T'.
 301 *>          They are provided for a developer/implementer who is familiar
 302 *>          with the details of the method.
 303 *>
 304 *>          WORK(6) = the entropy of A^t*A :: this is the Shannon entropy
 305 *>                    of diag(A^t*A) / Trace(A^t*A) taken as point in the
 306 *>                    probability simplex.
 307 *>          WORK(7) = the entropy of A*A^t.
 308 *> \endverbatim
 309 *>
 310 *> \param[in] LWORK
 311 *> \verbatim
 312 *>          LWORK is INTEGER
 313 *>          Length of WORK to confirm proper allocation of work space.
 314 *>          LWORK depends on the job:
 315 *>
 316 *>          If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
 317 *>            -> .. no scaled condition estimate required (JOBE.EQ.'N'):
 318 *>               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
 319 *>               ->> For optimal performance (blocked code) the optimal value
 320 *>               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
 321 *>               block size for DGEQP3 and DGEQRF.
 322 *>               In general, optimal LWORK is computed as
 323 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7).
 324 *>            -> .. an estimate of the scaled condition number of A is
 325 *>               required (JOBA='E', 'G'). In this case, LWORK is the maximum
 326 *>               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7).
 327 *>               ->> For optimal performance (blocked code) the optimal value
 328 *>               is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7).
 329 *>               In general, the optimal length LWORK is computed as
 330 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
 331 *>                                                     N+N*N+LWORK(DPOCON),7).
 332 *>
 333 *>          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
 334 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
 335 *>            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
 336 *>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,
 337 *>               DORMLQ. In general, the optimal length LWORK is computed as
 338 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON),
 339 *>                       N+LWORK(DGELQF), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)).
 340 *>
 341 *>          If SIGMA and the left singular vectors are needed
 342 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
 343 *>            -> For optimal performance:
 344 *>               if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
 345 *>               if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
 346 *>               where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
 347 *>               In general, the optimal length LWORK is computed as
 348 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
 349 *>                        2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
 350 *>               Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or
 351 *>               M*NB (for JOBU.EQ.'F').
 352 *>
 353 *>          If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and
 354 *>            -> if JOBV.EQ.'V'
 355 *>               the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
 356 *>            -> if JOBV.EQ.'J' the minimal requirement is
 357 *>               LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
 358 *>            -> For optimal performance, LWORK should be additionally
 359 *>               larger than N+M*NB, where NB is the optimal block size
 360 *>               for DORMQR.
 361 *> \endverbatim
 362 *>
 363 *> \param[out] IWORK
 364 *> \verbatim
 365 *>          IWORK is INTEGER array, dimension M+3*N.
 366 *>          On exit,
 367 *>          IWORK(1) = the numerical rank determined after the initial
 368 *>                     QR factorization with pivoting. See the descriptions
 369 *>                     of JOBA and JOBR.
 370 *>          IWORK(2) = the number of the computed nonzero singular values
 371 *>          IWORK(3) = if nonzero, a warning message:
 372 *>                     If IWORK(3).EQ.1 then some of the column norms of A
 373 *>                     were denormalized floats. The requested high accuracy
 374 *>                     is not warranted by the data.
 375 *> \endverbatim
 376 *>
 377 *> \param[out] INFO
 378 *> \verbatim
 379 *>          INFO is INTEGER
 380 *>           < 0  : if INFO = -i, then the i-th argument had an illegal value.
 381 *>           = 0 :  successful exit;
 382 *>           > 0 :  DGEJSV  did not converge in the maximal allowed number
 383 *>                  of sweeps. The computed values may be inaccurate.
 384 *> \endverbatim
 385 *
 386 *  Authors:
 387 *  ========
 388 *
 389 *> \author Univ. of Tennessee
 390 *> \author Univ. of California Berkeley
 391 *> \author Univ. of Colorado Denver
 392 *> \author NAG Ltd.
 393 *
 394 *> \date June 2016
 395 *
 396 *> \ingroup doubleGEsing
 397 *
 398 *> \par Further Details:
 399 *  =====================
 400 *>
 401 *> \verbatim
 402 *>
 403 *>  DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses DGEQP3,
 404 *>  DGEQRF, and DGELQF as preprocessors and preconditioners. Optionally, an
 405 *>  additional row pivoting can be used as a preprocessor, which in some
 406 *>  cases results in much higher accuracy. An example is matrix A with the
 407 *>  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
 408 *>  diagonal matrices and C is well-conditioned matrix. In that case, complete
 409 *>  pivoting in the first QR factorizations provides accuracy dependent on the
 410 *>  condition number of C, and independent of D1, D2. Such higher accuracy is
 411 *>  not completely understood theoretically, but it works well in practice.
 412 *>  Further, if A can be written as A = B*D, with well-conditioned B and some
 413 *>  diagonal D, then the high accuracy is guaranteed, both theoretically and
 414 *>  in software, independent of D. For more details see [1], [2].
 415 *>     The computational range for the singular values can be the full range
 416 *>  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
 417 *>  & LAPACK routines called by DGEJSV are implemented to work in that range.
 418 *>  If that is not the case, then the restriction for safe computation with
 419 *>  the singular values in the range of normalized IEEE numbers is that the
 420 *>  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
 421 *>  overflow. This code (DGEJSV) is best used in this restricted range,
 422 *>  meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are
 423 *>  returned as zeros. See JOBR for details on this.
 424 *>     Further, this implementation is somewhat slower than the one described
 425 *>  in [1,2] due to replacement of some non-LAPACK components, and because
 426 *>  the choice of some tuning parameters in the iterative part (DGESVJ) is
 427 *>  left to the implementer on a particular machine.
 428 *>     The rank revealing QR factorization (in this code: DGEQP3) should be
 429 *>  implemented as in [3]. We have a new version of DGEQP3 under development
 430 *>  that is more robust than the current one in LAPACK, with a cleaner cut in
 431 *>  rank defficient cases. It will be available in the SIGMA library [4].
 432 *>  If M is much larger than N, it is obvious that the initial QRF with
 433 *>  column pivoting can be preprocessed by the QRF without pivoting. That
 434 *>  well known trick is not used in DGEJSV because in some cases heavy row
 435 *>  weighting can be treated with complete pivoting. The overhead in cases
 436 *>  M much larger than N is then only due to pivoting, but the benefits in
 437 *>  terms of accuracy have prevailed. The implementer/user can incorporate
 438 *>  this extra QRF step easily. The implementer can also improve data movement
 439 *>  (matrix transpose, matrix copy, matrix transposed copy) - this
 440 *>  implementation of DGEJSV uses only the simplest, naive data movement.
 441 *> \endverbatim
 442 *
 443 *> \par Contributors:
 444 *  ==================
 445 *>
 446 *>  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
 447 *
 448 *> \par References:
 449 *  ================
 450 *>
 451 *> \verbatim
 452 *>
 453 *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
 454 *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
 455 *>     LAPACK Working note 169.
 456 *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
 457 *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
 458 *>     LAPACK Working note 170.
 459 *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
 460 *>     factorization software - a case study.
 461 *>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
 462 *>     LAPACK Working note 176.
 463 *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
 464 *>     QSVD, (H,K)-SVD computations.
 465 *>     Department of Mathematics, University of Zagreb, 2008.
 466 *> \endverbatim
 467 *
 468 *>  \par Bugs, examples and comments:
 469 *   =================================
 470 *>
 471 *>  Please report all bugs and send interesting examples and/or comments to
 472 *>  drmac@math.hr. Thank you.
 473 *>
 474 *  =====================================================================
 475       SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
 476      $                   M, N, A, LDA, SVA, U, LDU, V, LDV,
 477      $                   WORK, LWORK, IWORK, INFO )
 478 *
 479 *  -- LAPACK computational routine (version 3.6.1) --
 480 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 481 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 482 *     June 2016
 483 *
 484 *     .. Scalar Arguments ..
 485       IMPLICIT    NONE
 486       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
 487 *     ..
 488 *     .. Array Arguments ..
 489       DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
 490      $            WORK( LWORK )
 491       INTEGER     IWORK( * )
 492       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
 493 *     ..
 494 *
 495 *  ===========================================================================
 496 *
 497 *     .. Local Parameters ..
 498       DOUBLE PRECISION   ZERO,  ONE
 499       PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
 500 *     ..
 501 *     .. Local Scalars ..
 502       DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
 503      $        CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,  MAXPRJ, SCALEM,
 504      $        SCONDA, SFMIN,  SMALL,  TEMP1,  USCAL1, USCAL2, XSC
 505       INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
 506       LOGICAL ALMORT, DEFR,   ERREST, GOSCAL, JRACC,  KILL,   LSVEC,
 507      $        L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
 508      $        NOSCAL, ROWPIV, RSVEC,  TRANSP
 509 *     ..
 510 *     .. Intrinsic Functions ..
 511       INTRINSIC DABS, DLOG, MAX, MIN, DBLE, IDNINT, DSIGN, DSQRT
 512 *     ..
 513 *     .. External Functions ..
 514       DOUBLE PRECISION  DLAMCH, DNRM2
 515       INTEGER   IDAMAX
 516       LOGICAL   LSAME
 517       EXTERNAL  IDAMAX, LSAME, DLAMCH, DNRM2
 518 *     ..
 519 *     .. External Subroutines ..
 520       EXTERNAL  DCOPY,  DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL,
 521      $          DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ,
 522      $          DORMQR, DPOCON, DSCAL,  DSWAP,  DTRSM,  XERBLA
 523 *
 524       EXTERNAL  DGESVJ
 525 *     ..
 526 *
 527 *     Test the input arguments
 528 *
 529       LSVEC  = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
 530       JRACC  = LSAME( JOBV, 'J' )
 531       RSVEC  = LSAME( JOBV, 'V' ) .OR. JRACC
 532       ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
 533       L2RANK = LSAME( JOBA, 'R' )
 534       L2ABER = LSAME( JOBA, 'A' )
 535       ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
 536       L2TRAN = LSAME( JOBT, 'T' )
 537       L2KILL = LSAME( JOBR, 'R' )
 538       DEFR   = LSAME( JOBR, 'N' )
 539       L2PERT = LSAME( JOBP, 'P' )
 540 *
 541       IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
 542      $     ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
 543          INFO = - 1
 544       ELSE IF ( .NOT.( LSVEC  .OR. LSAME( JOBU, 'N' ) .OR.
 545      $                             LSAME( JOBU, 'W' )) ) THEN
 546          INFO = - 2
 547       ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
 548      $   LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN
 549          INFO = - 3
 550       ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) )    THEN
 551          INFO = - 4
 552       ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN
 553          INFO = - 5
 554       ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
 555          INFO = - 6
 556       ELSE IF ( M .LT. 0 ) THEN
 557          INFO = - 7
 558       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
 559          INFO = - 8
 560       ELSE IF ( LDA .LT. M ) THEN
 561          INFO = - 10
 562       ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
 563          INFO = - 13
 564       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
 565          INFO = - 14
 566       ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
 567      &                           (LWORK .LT. MAX(7,4*N+1,2*M+N))) .OR.
 568      & (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND.
 569      &                         (LWORK .LT. MAX(7,4*N+N*N,2*M+N))) .OR.
 570      & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1)))
 571      & .OR.
 572      & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1)))
 573      & .OR.
 574      & (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND.
 575      &                          (LWORK.LT.MAX(2*M+N,6*N+2*N*N)))
 576      & .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND.
 577      &                          LWORK.LT.MAX(2*M+N,4*N+N*N,2*N+N*N+6)))
 578      &   THEN
 579          INFO = - 17
 580       ELSE
 581 *        #:)
 582          INFO = 0
 583       END IF
 584 *
 585       IF ( INFO .NE. 0 ) THEN
 586 *       #:(
 587          CALL XERBLA( 'DGEJSV', - INFO )
 588          RETURN
 589       END IF
 590 *
 591 *     Quick return for void matrix (Y3K safe)
 592 * #:)
 593       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
 594          IWORK(1:3) = 0
 595          WORK(1:7) = 0
 596          RETURN
 597       ENDIF
 598 *
 599 *     Determine whether the matrix U should be M x N or M x M
 600 *
 601       IF ( LSVEC ) THEN
 602          N1 = N
 603          IF ( LSAME( JOBU, 'F' ) ) N1 = M
 604       END IF
 605 *
 606 *     Set numerical parameters
 607 *
 608 *!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
 609 *
 610       EPSLN = DLAMCH('Epsilon')
 611       SFMIN = DLAMCH('SafeMinimum')
 612       SMALL = SFMIN / EPSLN
 613       BIG   = DLAMCH('O')
 614 *     BIG   = ONE / SFMIN
 615 *
 616 *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
 617 *
 618 *(!)  If necessary, scale SVA() to protect the largest norm from
 619 *     overflow. It is possible that this scaling pushes the smallest
 620 *     column norm left from the underflow threshold (extreme case).
 621 *
 622       SCALEM  = ONE / DSQRT(DBLE(M)*DBLE(N))
 623       NOSCAL  = .TRUE.
 624       GOSCAL  = .TRUE.
 625       DO 1874 p = 1, N
 626          AAPP = ZERO
 627          AAQQ = ONE
 628          CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ )
 629          IF ( AAPP .GT. BIG ) THEN
 630             INFO = - 9
 631             CALL XERBLA( 'DGEJSV', -INFO )
 632             RETURN
 633          END IF
 634          AAQQ = DSQRT(AAQQ)
 635          IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL  ) THEN
 636             SVA(p)  = AAPP * AAQQ
 637          ELSE
 638             NOSCAL  = .FALSE.
 639             SVA(p)  = AAPP * ( AAQQ * SCALEM )
 640             IF ( GOSCAL ) THEN
 641                GOSCAL = .FALSE.
 642                CALL DSCAL( p-1, SCALEM, SVA, 1 )
 643             END IF
 644          END IF
 645  1874 CONTINUE
 646 *
 647       IF ( NOSCAL ) SCALEM = ONE
 648 *
 649       AAPP = ZERO
 650       AAQQ = BIG
 651       DO 4781 p = 1, N
 652          AAPP = MAX( AAPP, SVA(p) )
 653          IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) )
 654  4781 CONTINUE
 655 *
 656 *     Quick return for zero M x N matrix
 657 * #:)
 658       IF ( AAPP .EQ. ZERO ) THEN
 659          IF ( LSVEC ) CALL DLASET( 'G', M, N1, ZERO, ONE, U, LDU )
 660          IF ( RSVEC ) CALL DLASET( 'G', N, N,  ZERO, ONE, V, LDV )
 661          WORK(1) = ONE
 662          WORK(2) = ONE
 663          IF ( ERREST ) WORK(3) = ONE
 664          IF ( LSVEC .AND. RSVEC ) THEN
 665             WORK(4) = ONE
 666             WORK(5) = ONE
 667          END IF
 668          IF ( L2TRAN ) THEN
 669             WORK(6) = ZERO
 670             WORK(7) = ZERO
 671          END IF
 672          IWORK(1) = 0
 673          IWORK(2) = 0
 674          IWORK(3) = 0
 675          RETURN
 676       END IF
 677 *
 678 *     Issue warning if denormalized column norms detected. Override the
 679 *     high relative accuracy request. Issue licence to kill columns
 680 *     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
 681 * #:(
 682       WARNING = 0
 683       IF ( AAQQ .LE. SFMIN ) THEN
 684          L2RANK = .TRUE.
 685          L2KILL = .TRUE.
 686          WARNING = 1
 687       END IF
 688 *
 689 *     Quick return for one-column matrix
 690 * #:)
 691       IF ( N .EQ. 1 ) THEN
 692 *
 693          IF ( LSVEC ) THEN
 694             CALL DLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
 695             CALL DLACPY( 'A', M, 1, A, LDA, U, LDU )
 696 *           computing all M left singular vectors of the M x 1 matrix
 697             IF ( N1 .NE. N  ) THEN
 698                CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
 699                CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
 700                CALL DCOPY( M, A(1,1), 1, U(1,1), 1 )
 701             END IF
 702          END IF
 703          IF ( RSVEC ) THEN
 704              V(1,1) = ONE
 705          END IF
 706          IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
 707             SVA(1)  = SVA(1) / SCALEM
 708             SCALEM  = ONE
 709          END IF
 710          WORK(1) = ONE / SCALEM
 711          WORK(2) = ONE
 712          IF ( SVA(1) .NE. ZERO ) THEN
 713             IWORK(1) = 1
 714             IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
 715                IWORK(2) = 1
 716             ELSE
 717                IWORK(2) = 0
 718             END IF
 719          ELSE
 720             IWORK(1) = 0
 721             IWORK(2) = 0
 722          END IF
 723          IWORK(3) = 0
 724          IF ( ERREST ) WORK(3) = ONE
 725          IF ( LSVEC .AND. RSVEC ) THEN
 726             WORK(4) = ONE
 727             WORK(5) = ONE
 728          END IF
 729          IF ( L2TRAN ) THEN
 730             WORK(6) = ZERO
 731             WORK(7) = ZERO
 732          END IF
 733          RETURN
 734 *
 735       END IF
 736 *
 737       TRANSP = .FALSE.
 738       L2TRAN = L2TRAN .AND. ( M .EQ. N )
 739 *
 740       AATMAX = -ONE
 741       AATMIN =  BIG
 742       IF ( ROWPIV .OR. L2TRAN ) THEN
 743 *
 744 *     Compute the row norms, needed to determine row pivoting sequence
 745 *     (in the case of heavily row weighted A, row pivoting is strongly
 746 *     advised) and to collect information needed to compare the
 747 *     structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.).
 748 *
 749          IF ( L2TRAN ) THEN
 750             DO 1950 p = 1, M
 751                XSC   = ZERO
 752                TEMP1 = ONE
 753                CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
 754 *              DLASSQ gets both the ell_2 and the ell_infinity norm
 755 *              in one pass through the vector
 756                WORK(M+N+p)  = XSC * SCALEM
 757                WORK(N+p)    = XSC * (SCALEM*DSQRT(TEMP1))
 758                AATMAX = MAX( AATMAX, WORK(N+p) )
 759                IF (WORK(N+p) .NE. ZERO) AATMIN = MIN(AATMIN,WORK(N+p))
 760  1950       CONTINUE
 761          ELSE
 762             DO 1904 p = 1, M
 763                WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) )
 764                AATMAX = MAX( AATMAX, WORK(M+N+p) )
 765                AATMIN = MIN( AATMIN, WORK(M+N+p) )
 766  1904       CONTINUE
 767          END IF
 768 *
 769       END IF
 770 *
 771 *     For square matrix A try to determine whether A^t  would be  better
 772 *     input for the preconditioned Jacobi SVD, with faster convergence.
 773 *     The decision is based on an O(N) function of the vector of column
 774 *     and row norms of A, based on the Shannon entropy. This should give
 775 *     the right choice in most cases when the difference actually matters.
 776 *     It may fail and pick the slower converging side.
 777 *
 778       ENTRA  = ZERO
 779       ENTRAT = ZERO
 780       IF ( L2TRAN ) THEN
 781 *
 782          XSC   = ZERO
 783          TEMP1 = ONE
 784          CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
 785          TEMP1 = ONE / TEMP1
 786 *
 787          ENTRA = ZERO
 788          DO 1113 p = 1, N
 789             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
 790             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
 791  1113    CONTINUE
 792          ENTRA = - ENTRA / DLOG(DBLE(N))
 793 *
 794 *        Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex.
 795 *        It is derived from the diagonal of  A^t * A.  Do the same with the
 796 *        diagonal of A * A^t, compute the entropy of the corresponding
 797 *        probability distribution. Note that A * A^t and A^t * A have the
 798 *        same trace.
 799 *
 800          ENTRAT = ZERO
 801          DO 1114 p = N+1, N+M
 802             BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1
 803             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
 804  1114    CONTINUE
 805          ENTRAT = - ENTRAT / DLOG(DBLE(M))
 806 *
 807 *        Analyze the entropies and decide A or A^t. Smaller entropy
 808 *        usually means better input for the algorithm.
 809 *
 810          TRANSP = ( ENTRAT .LT. ENTRA )
 811 *
 812 *        If A^t is better than A, transpose A.
 813 *
 814          IF ( TRANSP ) THEN
 815 *           In an optimal implementation, this trivial transpose
 816 *           should be replaced with faster transpose.
 817             DO 1115 p = 1, N - 1
 818                DO 1116 q = p + 1, N
 819                    TEMP1 = A(q,p)
 820                   A(q,p) = A(p,q)
 821                   A(p,q) = TEMP1
 822  1116          CONTINUE
 823  1115       CONTINUE
 824             DO 1117 p = 1, N
 825                WORK(M+N+p) = SVA(p)
 826                SVA(p)      = WORK(N+p)
 827  1117       CONTINUE
 828             TEMP1  = AAPP
 829             AAPP   = AATMAX
 830             AATMAX = TEMP1
 831             TEMP1  = AAQQ
 832             AAQQ   = AATMIN
 833             AATMIN = TEMP1
 834             KILL   = LSVEC
 835             LSVEC  = RSVEC
 836             RSVEC  = KILL
 837             IF ( LSVEC ) N1 = N
 838 *
 839             ROWPIV = .TRUE.
 840          END IF
 841 *
 842       END IF
 843 *     END IF L2TRAN
 844 *
 845 *     Scale the matrix so that its maximal singular value remains less
 846 *     than DSQRT(BIG) -- the matrix is scaled so that its maximal column
 847 *     has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep
 848 *     DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and
 849 *     BLAS routines that, in some implementations, are not capable of
 850 *     working in the full interval [SFMIN,BIG] and that they may provoke
 851 *     overflows in the intermediate results. If the singular values spread
 852 *     from SFMIN to BIG, then DGESVJ will compute them. So, in that case,
 853 *     one should use DGESVJ instead of DGEJSV.
 854 *
 855       BIG1   = DSQRT( BIG )
 856       TEMP1  = DSQRT( BIG / DBLE(N) )
 857 *
 858       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
 859       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
 860           AAQQ = ( AAQQ / AAPP ) * TEMP1
 861       ELSE
 862           AAQQ = ( AAQQ * TEMP1 ) / AAPP
 863       END IF
 864       TEMP1 = TEMP1 * SCALEM
 865       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
 866 *
 867 *     To undo scaling at the end of this procedure, multiply the
 868 *     computed singular values with USCAL2 / USCAL1.
 869 *
 870       USCAL1 = TEMP1
 871       USCAL2 = AAPP
 872 *
 873       IF ( L2KILL ) THEN
 874 *        L2KILL enforces computation of nonzero singular values in
 875 *        the restricted range of condition number of the initial A,
 876 *        sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN).
 877          XSC = DSQRT( SFMIN )
 878       ELSE
 879          XSC = SMALL
 880 *
 881 *        Now, if the condition number of A is too big,
 882 *        sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN,
 883 *        as a precaution measure, the full SVD is computed using DGESVJ
 884 *        with accumulated Jacobi rotations. This provides numerically
 885 *        more robust computation, at the cost of slightly increased run
 886 *        time. Depending on the concrete implementation of BLAS and LAPACK
 887 *        (i.e. how they behave in presence of extreme ill-conditioning) the
 888 *        implementor may decide to remove this switch.
 889          IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
 890             JRACC = .TRUE.
 891          END IF
 892 *
 893       END IF
 894       IF ( AAQQ .LT. XSC ) THEN
 895          DO 700 p = 1, N
 896             IF ( SVA(p) .LT. XSC ) THEN
 897                CALL DLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA )
 898                SVA(p) = ZERO
 899             END IF
 900  700     CONTINUE
 901       END IF
 902 *
 903 *     Preconditioning using QR factorization with pivoting
 904 *
 905       IF ( ROWPIV ) THEN
 906 *        Optional row permutation (Bjoerck row pivoting):
 907 *        A result by Cox and Higham shows that the Bjoerck's
 908 *        row pivoting combined with standard column pivoting
 909 *        has similar effect as Powell-Reid complete pivoting.
 910 *        The ell-infinity norms of A are made nonincreasing.
 911          DO 1952 p = 1, M - 1
 912             q = IDAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1
 913             IWORK(2*N+p) = q
 914             IF ( p .NE. q ) THEN
 915                TEMP1       = WORK(M+N+p)
 916                WORK(M+N+p) = WORK(M+N+q)
 917                WORK(M+N+q) = TEMP1
 918             END IF
 919  1952    CONTINUE
 920          CALL DLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 )
 921       END IF
 922 *
 923 *     End of the preparation phase (scaling, optional sorting and
 924 *     transposing, optional flushing of small columns).
 925 *
 926 *     Preconditioning
 927 *
 928 *     If the full SVD is needed, the right singular vectors are computed
 929 *     from a matrix equation, and for that we need theoretical analysis
 930 *     of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF.
 931 *     In all other cases the first RR QRF can be chosen by other criteria
 932 *     (eg speed by replacing global with restricted window pivoting, such
 933 *     as in SGEQPX from TOMS # 782). Good results will be obtained using
 934 *     SGEQPX with properly (!) chosen numerical parameters.
 935 *     Any improvement of DGEQP3 improves overal performance of DGEJSV.
 936 *
 937 *     A * P1 = Q1 * [ R1^t 0]^t:
 938       DO 1963 p = 1, N
 939 *        .. all columns are free columns
 940          IWORK(p) = 0
 941  1963 CONTINUE
 942       CALL DGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR )
 943 *
 944 *     The upper triangular matrix R1 from the first QRF is inspected for
 945 *     rank deficiency and possibilities for deflation, or possible
 946 *     ill-conditioning. Depending on the user specified flag L2RANK,
 947 *     the procedure explores possibilities to reduce the numerical
 948 *     rank by inspecting the computed upper triangular factor. If
 949 *     L2RANK or L2ABER are up, then DGEJSV will compute the SVD of
 950 *     A + dA, where ||dA|| <= f(M,N)*EPSLN.
 951 *
 952       NR = 1
 953       IF ( L2ABER ) THEN
 954 *        Standard absolute error bound suffices. All sigma_i with
 955 *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
 956 *        agressive enforcement of lower numerical rank by introducing a
 957 *        backward error of the order of N*EPSLN*||A||.
 958          TEMP1 = DSQRT(DBLE(N))*EPSLN
 959          DO 3001 p = 2, N
 960             IF ( DABS(A(p,p)) .GE. (TEMP1*DABS(A(1,1))) ) THEN
 961                NR = NR + 1
 962             ELSE
 963                GO TO 3002
 964             END IF
 965  3001    CONTINUE
 966  3002    CONTINUE
 967       ELSE IF ( L2RANK ) THEN
 968 *        .. similarly as above, only slightly more gentle (less agressive).
 969 *        Sudden drop on the diagonal of R1 is used as the criterion for
 970 *        close-to-rank-defficient.
 971          TEMP1 = DSQRT(SFMIN)
 972          DO 3401 p = 2, N
 973             IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR.
 974      $           ( DABS(A(p,p)) .LT. SMALL ) .OR.
 975      $           ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
 976             NR = NR + 1
 977  3401    CONTINUE
 978  3402    CONTINUE
 979 *
 980       ELSE
 981 *        The goal is high relative accuracy. However, if the matrix
 982 *        has high scaled condition number the relative accuracy is in
 983 *        general not feasible. Later on, a condition number estimator
 984 *        will be deployed to estimate the scaled condition number.
 985 *        Here we just remove the underflowed part of the triangular
 986 *        factor. This prevents the situation in which the code is
 987 *        working hard to get the accuracy not warranted by the data.
 988          TEMP1  = DSQRT(SFMIN)
 989          DO 3301 p = 2, N
 990             IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR.
 991      $          ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
 992             NR = NR + 1
 993  3301    CONTINUE
 994  3302    CONTINUE
 995 *
 996       END IF
 997 *
 998       ALMORT = .FALSE.
 999       IF ( NR .EQ. N ) THEN
1000          MAXPRJ = ONE
1001          DO 3051 p = 2, N
1002             TEMP1  = DABS(A(p,p)) / SVA(IWORK(p))
1003             MAXPRJ = MIN( MAXPRJ, TEMP1 )
1004  3051    CONTINUE
1005          IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
1006       END IF
1007 *
1008 *
1009       SCONDA = - ONE
1010       CONDR1 = - ONE
1011       CONDR2 = - ONE
1012 *
1013       IF ( ERREST ) THEN
1014          IF ( N .EQ. NR ) THEN
1015             IF ( RSVEC ) THEN
1016 *              .. V is available as workspace
1017                CALL DLACPY( 'U', N, N, A, LDA, V, LDV )
1018                DO 3053 p = 1, N
1019                   TEMP1 = SVA(IWORK(p))
1020                   CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 )
1021  3053          CONTINUE
1022                CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1,
1023      $              WORK(N+1), IWORK(2*N+M+1), IERR )
1024             ELSE IF ( LSVEC ) THEN
1025 *              .. U is available as workspace
1026                CALL DLACPY( 'U', N, N, A, LDA, U, LDU )
1027                DO 3054 p = 1, N
1028                   TEMP1 = SVA(IWORK(p))
1029                   CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 )
1030  3054          CONTINUE
1031                CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1,
1032      $              WORK(N+1), IWORK(2*N+M+1), IERR )
1033             ELSE
1034                CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N )
1035                DO 3052 p = 1, N
1036                   TEMP1 = SVA(IWORK(p))
1037                   CALL DSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 )
1038  3052          CONTINUE
1039 *           .. the columns of R are scaled to have unit Euclidean lengths.
1040                CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1,
1041      $              WORK(N+N*N+1), IWORK(2*N+M+1), IERR )
1042             END IF
1043             SCONDA = ONE / DSQRT(TEMP1)
1044 *           SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
1045 *           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
1046          ELSE
1047             SCONDA = - ONE
1048          END IF
1049       END IF
1050 *
1051       L2PERT = L2PERT .AND. ( DABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) )
1052 *     If there is no violent scaling, artificial perturbation is not needed.
1053 *
1054 *     Phase 3:
1055 *
1056       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
1057 *
1058 *         Singular Values only
1059 *
1060 *         .. transpose A(1:NR,1:N)
1061          DO 1946 p = 1, MIN( N-1, NR )
1062             CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
1063  1946    CONTINUE
1064 *
1065 *        The following two DO-loops introduce small relative perturbation
1066 *        into the strict upper triangle of the lower triangular matrix.
1067 *        Small entries below the main diagonal are also changed.
1068 *        This modification is useful if the computing environment does not
1069 *        provide/allow FLUSH TO ZERO underflow, for it prevents many
1070 *        annoying denormalized numbers in case of strongly scaled matrices.
1071 *        The perturbation is structured so that it does not introduce any
1072 *        new perturbation of the singular values, and it does not destroy
1073 *        the job done by the preconditioner.
1074 *        The licence for this perturbation is in the variable L2PERT, which
1075 *        should be .FALSE. if FLUSH TO ZERO underflow is active.
1076 *
1077          IF ( .NOT. ALMORT ) THEN
1078 *
1079             IF ( L2PERT ) THEN
1080 *              XSC = DSQRT(SMALL)
1081                XSC = EPSLN / DBLE(N)
1082                DO 4947 q = 1, NR
1083                   TEMP1 = XSC*DABS(A(q,q))
1084                   DO 4949 p = 1, N
1085                      IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
1086      $                    .OR. ( p .LT. q ) )
1087      $                     A(p,q) = DSIGN( TEMP1, A(p,q) )
1088  4949             CONTINUE
1089  4947          CONTINUE
1090             ELSE
1091                CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA )
1092             END IF
1093 *
1094 *            .. second preconditioning using the QR factorization
1095 *
1096             CALL DGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR )
1097 *
1098 *           .. and transpose upper to lower triangular
1099             DO 1948 p = 1, NR - 1
1100                CALL DCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
1101  1948       CONTINUE
1102 *
1103          END IF
1104 *
1105 *           Row-cyclic Jacobi SVD algorithm with column pivoting
1106 *
1107 *           .. again some perturbation (a "background noise") is added
1108 *           to drown denormals
1109             IF ( L2PERT ) THEN
1110 *              XSC = DSQRT(SMALL)
1111                XSC = EPSLN / DBLE(N)
1112                DO 1947 q = 1, NR
1113                   TEMP1 = XSC*DABS(A(q,q))
1114                   DO 1949 p = 1, NR
1115                      IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
1116      $                       .OR. ( p .LT. q ) )
1117      $                   A(p,q) = DSIGN( TEMP1, A(p,q) )
1118  1949             CONTINUE
1119  1947          CONTINUE
1120             ELSE
1121                CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
1122             END IF
1123 *
1124 *           .. and one-sided Jacobi rotations are started on a lower
1125 *           triangular matrix (plus perturbation which is ignored in
1126 *           the part which destroys triangular form (confusing?!))
1127 *
1128             CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA,
1129      $                      N, V, LDV, WORK, LWORK, INFO )
1130 *
1131             SCALEM  = WORK(1)
1132             NUMRANK = IDNINT(WORK(2))
1133 *
1134 *
1135       ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
1136 *
1137 *        -> Singular Values and Right Singular Vectors <-
1138 *
1139          IF ( ALMORT ) THEN
1140 *
1141 *           .. in this case NR equals N
1142             DO 1998 p = 1, NR
1143                CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1144  1998       CONTINUE
1145             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1146 *
1147             CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
1148      $                  WORK, LWORK, INFO )
1149             SCALEM  = WORK(1)
1150             NUMRANK = IDNINT(WORK(2))
1151
1152          ELSE
1153 *
1154 *        .. two more QR factorizations ( one QRF is not enough, two require
1155 *        accumulated product of Jacobi rotations, three are perfect )
1156 *
1157             CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
1158             CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
1159             CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
1160             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1161             CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1162      $                   LWORK-2*N, IERR )
1163             DO 8998 p = 1, NR
1164                CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
1165  8998       CONTINUE
1166             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1167 *
1168             CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
1169      $                  LDU, WORK(N+1), LWORK, INFO )
1170             SCALEM  = WORK(N+1)
1171             NUMRANK = IDNINT(WORK(N+2))
1172             IF ( NR .LT. N ) THEN
1173                CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),   LDV )
1174                CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),   LDV )
1175                CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
1176             END IF
1177 *
1178          CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
1179      $               V, LDV, WORK(N+1), LWORK-N, IERR )
1180 *
1181          END IF
1182 *
1183          DO 8991 p = 1, N
1184             CALL DCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
1185  8991    CONTINUE
1186          CALL DLACPY( 'All', N, N, A, LDA, V, LDV )
1187 *
1188          IF ( TRANSP ) THEN
1189             CALL DLACPY( 'All', N, N, V, LDV, U, LDU )
1190          END IF
1191 *
1192       ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
1193 *
1194 *        .. Singular Values and Left Singular Vectors                 ..
1195 *
1196 *        .. second preconditioning step to avoid need to accumulate
1197 *        Jacobi rotations in the Jacobi iterations.
1198          DO 1965 p = 1, NR
1199             CALL DCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
1200  1965    CONTINUE
1201          CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1202 *
1203          CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1),
1204      $              LWORK-2*N, IERR )
1205 *
1206          DO 1967 p = 1, NR - 1
1207             CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
1208  1967    CONTINUE
1209          CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1210 *
1211          CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
1212      $        LDA, WORK(N+1), LWORK-N, INFO )
1213          SCALEM  = WORK(N+1)
1214          NUMRANK = IDNINT(WORK(N+2))
1215 *
1216          IF ( NR .LT. M ) THEN
1217             CALL DLASET( 'A',  M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
1218             IF ( NR .LT. N1 ) THEN
1219                CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
1220                CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
1221             END IF
1222          END IF
1223 *
1224          CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1225      $               LDU, WORK(N+1), LWORK-N, IERR )
1226 *
1227          IF ( ROWPIV )
1228      $       CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1229 *
1230          DO 1974 p = 1, N1
1231             XSC = ONE / DNRM2( M, U(1,p), 1 )
1232             CALL DSCAL( M, XSC, U(1,p), 1 )
1233  1974    CONTINUE
1234 *
1235          IF ( TRANSP ) THEN
1236             CALL DLACPY( 'All', N, N, U, LDU, V, LDV )
1237          END IF
1238 *
1239       ELSE
1240 *
1241 *        .. Full SVD ..
1242 *
1243          IF ( .NOT. JRACC ) THEN
1244 *
1245          IF ( .NOT. ALMORT ) THEN
1246 *
1247 *           Second Preconditioning Step (QRF [with pivoting])
1248 *           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
1249 *           equivalent to an LQF CALL. Since in many libraries the QRF
1250 *           seems to be better optimized than the LQF, we do explicit
1251 *           transpose and use the QRF. This is subject to changes in an
1252 *           optimized implementation of DGEJSV.
1253 *
1254             DO 1968 p = 1, NR
1255                CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1256  1968       CONTINUE
1257 *
1258 *           .. the following two loops perturb small entries to avoid
1259 *           denormals in the second QR factorization, where they are
1260 *           as good as zeros. This is done to avoid painfully slow
1261 *           computation with denormals. The relative size of the perturbation
1262 *           is a parameter that can be changed by the implementer.
1263 *           This perturbation device will be obsolete on machines with
1264 *           properly implemented arithmetic.
1265 *           To switch it off, set L2PERT=.FALSE. To remove it from  the
1266 *           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
1267 *           The following two loops should be blocked and fused with the
1268 *           transposed copy above.
1269 *
1270             IF ( L2PERT ) THEN
1271                XSC = DSQRT(SMALL)
1272                DO 2969 q = 1, NR
1273                   TEMP1 = XSC*DABS( V(q,q) )
1274                   DO 2968 p = 1, N
1275                      IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
1276      $                   .OR. ( p .LT. q ) )
1277      $                   V(p,q) = DSIGN( TEMP1, V(p,q) )
1278                      IF ( p .LT. q ) V(p,q) = - V(p,q)
1279  2968             CONTINUE
1280  2969          CONTINUE
1281             ELSE
1282                CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1283             END IF
1284 *
1285 *           Estimate the row scaled condition number of R1
1286 *           (If R1 is rectangular, N > NR, then the condition number
1287 *           of the leading NR x NR submatrix is estimated.)
1288 *
1289             CALL DLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR )
1290             DO 3950 p = 1, NR
1291                TEMP1 = DNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1)
1292                CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1)
1293  3950       CONTINUE
1294             CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
1295      $                   WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
1296             CONDR1 = ONE / DSQRT(TEMP1)
1297 *           .. here need a second oppinion on the condition number
1298 *           .. then assume worst case scenario
1299 *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
1300 *           more conservative    <=> CONDR1 .LT. DSQRT(DBLE(N))
1301 *
1302             COND_OK = DSQRT(DBLE(NR))
1303 *[TP]       COND_OK is a tuning parameter.
1304
1305             IF ( CONDR1 .LT. COND_OK ) THEN
1306 *              .. the second QRF without pivoting. Note: in an optimized
1307 *              implementation, this QRF should be implemented as the QRF
1308 *              of a lower triangular matrix.
1309 *              R1^t = Q2 * R2
1310                CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1311      $              LWORK-2*N, IERR )
1312 *
1313                IF ( L2PERT ) THEN
1314                   XSC = DSQRT(SMALL)/EPSLN
1315                   DO 3959 p = 2, NR
1316                      DO 3958 q = 1, p - 1
1317                         TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
1318                         IF ( DABS(V(q,p)) .LE. TEMP1 )
1319      $                     V(q,p) = DSIGN( TEMP1, V(q,p) )
1320  3958                CONTINUE
1321  3959             CONTINUE
1322                END IF
1323 *
1324                IF ( NR .NE. N )
1325      $         CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
1326 *              .. save ...
1327 *
1328 *           .. this transposed copy should be better than naive
1329                DO 1969 p = 1, NR - 1
1330                   CALL DCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
1331  1969          CONTINUE
1332 *
1333                CONDR2 = CONDR1
1334 *
1335             ELSE
1336 *
1337 *              .. ill-conditioned case: second QRF with pivoting
1338 *              Note that windowed pivoting would be equaly good
1339 *              numerically, and more run-time efficient. So, in
1340 *              an optimal implementation, the next call to DGEQP3
1341 *              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
1342 *              with properly (carefully) chosen parameters.
1343 *
1344 *              R1^t * P2 = Q2 * R2
1345                DO 3003 p = 1, NR
1346                   IWORK(N+p) = 0
1347  3003          CONTINUE
1348                CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1),
1349      $                  WORK(2*N+1), LWORK-2*N, IERR )
1350 **               CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1351 **     $              LWORK-2*N, IERR )
1352                IF ( L2PERT ) THEN
1353                   XSC = DSQRT(SMALL)
1354                   DO 3969 p = 2, NR
1355                      DO 3968 q = 1, p - 1
1356                         TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
1357                         IF ( DABS(V(q,p)) .LE. TEMP1 )
1358      $                     V(q,p) = DSIGN( TEMP1, V(q,p) )
1359  3968                CONTINUE
1360  3969             CONTINUE
1361                END IF
1362 *
1363                CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
1364 *
1365                IF ( L2PERT ) THEN
1366                   XSC = DSQRT(SMALL)
1367                   DO 8970 p = 2, NR
1368                      DO 8971 q = 1, p - 1
1369                         TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q)))
1370                         V(p,q) = - DSIGN( TEMP1, V(q,p) )
1371  8971                CONTINUE
1372  8970             CONTINUE
1373                ELSE
1374                   CALL DLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV )
1375                END IF
1376 *              Now, compute R2 = L3 * Q3, the LQ factorization.
1377                CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1),
1378      $               WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
1379 *              .. and estimate the condition number
1380                CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR )
1381                DO 4950 p = 1, NR
1382                   TEMP1 = DNRM2( p, WORK(2*N+N*NR+NR+p), NR )
1383                   CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR )
1384  4950          CONTINUE
1385                CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
1386      $              WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR )
1387                CONDR2 = ONE / DSQRT(TEMP1)
1388 *
1389                IF ( CONDR2 .GE. COND_OK ) THEN
1390 *                 .. save the Householder vectors used for Q3
1391 *                 (this overwrittes the copy of R2, as it will not be
1392 *                 needed in this branch, but it does not overwritte the
1393 *                 Huseholder vectors of Q2.).
1394                   CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
1395 *                 .. and the rest of the information on Q3 is in
1396 *                 WORK(2*N+N*NR+1:2*N+N*NR+N)
1397                END IF
1398 *
1399             END IF
1400 *
1401             IF ( L2PERT ) THEN
1402                XSC = DSQRT(SMALL)
1403                DO 4968 q = 2, NR
1404                   TEMP1 = XSC * V(q,q)
1405                   DO 4969 p = 1, q - 1
1406 *                    V(p,q) = - DSIGN( TEMP1, V(q,p) )
1407                      V(p,q) = - DSIGN( TEMP1, V(p,q) )
1408  4969             CONTINUE
1409  4968          CONTINUE
1410             ELSE
1411                CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
1412             END IF
1413 *
1414 *        Second preconditioning finished; continue with Jacobi SVD
1415 *        The input matrix is lower trinagular.
1416 *
1417 *        Recover the right singular vectors as solution of a well
1418 *        conditioned triangular matrix equation.
1419 *
1420             IF ( CONDR1 .LT. COND_OK ) THEN
1421 *
1422                CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U,
1423      $              LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO )
1424                SCALEM  = WORK(2*N+N*NR+NR+1)
1425                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1426                DO 3970 p = 1, NR
1427                   CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
1428                   CALL DSCAL( NR, SVA(p),    V(1,p), 1 )
1429  3970          CONTINUE
1430
1431 *        .. pick the right matrix equation and solve it
1432 *
1433                IF ( NR .EQ. N ) THEN
1434 * :))             .. best case, R1 is inverted. The solution of this matrix
1435 *                 equation is Q2*V2 = the product of the Jacobi rotations
1436 *                 used in DGESVJ, premultiplied with the orthogonal matrix
1437 *                 from the second QR factorization.
1438                   CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
1439                ELSE
1440 *                 .. R1 is well conditioned, but non-square. Transpose(R2)
1441 *                 is inverted to get the product of the Jacobi rotations
1442 *                 used in DGESVJ. The Q-factor from the second QR
1443 *                 factorization is then built in explicitly.
1444                   CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1),
1445      $                 N,V,LDV)
1446                   IF ( NR .LT. N ) THEN
1447                     CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
1448                     CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
1449                     CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
1450                   END IF
1451                   CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1452      $                 V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
1453                END IF
1454 *
1455             ELSE IF ( CONDR2 .LT. COND_OK ) THEN
1456 *
1457 * :)           .. the input matrix A is very likely a relative of
1458 *              the Kahan matrix :)
1459 *              The matrix R2 is inverted. The solution of the matrix equation
1460 *              is Q3^T*V3 = the product of the Jacobi rotations (appplied to
1461 *              the lower triangular L3 from the LQ factorization of
1462 *              R2=L3*Q3), pre-multiplied with the transposed Q3.
1463                CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
1464      $              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
1465                SCALEM  = WORK(2*N+N*NR+NR+1)
1466                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1467                DO 3870 p = 1, NR
1468                   CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
1469                   CALL DSCAL( NR, SVA(p),    U(1,p), 1 )
1470  3870          CONTINUE
1471                CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
1472 *              .. apply the permutation from the second QR factorization
1473                DO 873 q = 1, NR
1474                   DO 872 p = 1, NR
1475                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1476  872              CONTINUE
1477                   DO 874 p = 1, NR
1478                      U(p,q) = WORK(2*N+N*NR+NR+p)
1479  874              CONTINUE
1480  873           CONTINUE
1481                IF ( NR .LT. N ) THEN
1482                   CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1483                   CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1484                   CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1485                END IF
1486                CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1487      $              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1488             ELSE
1489 *              Last line of defense.
1490 * #:(          This is a rather pathological case: no scaled condition
1491 *              improvement after two pivoted QR factorizations. Other
1492 *              possibility is that the rank revealing QR factorization
1493 *              or the condition estimator has failed, or the COND_OK
1494 *              is set very close to ONE (which is unnecessary). Normally,
1495 *              this branch should never be executed, but in rare cases of
1496 *              failure of the RRQR or condition estimator, the last line of
1497 *              defense ensures that DGEJSV completes the task.
1498 *              Compute the full SVD of L3 using DGESVJ with explicit
1499 *              accumulation of Jacobi rotations.
1500                CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
1501      $              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
1502                SCALEM  = WORK(2*N+N*NR+NR+1)
1503                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1504                IF ( NR .LT. N ) THEN
1505                   CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1506                   CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1507                   CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1508                END IF
1509                CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1510      $              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1511 *
1512                CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N,
1513      $              WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1),
1514      $              LWORK-2*N-N*NR-NR, IERR )
1515                DO 773 q = 1, NR
1516                   DO 772 p = 1, NR
1517                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1518  772              CONTINUE
1519                   DO 774 p = 1, NR
1520                      U(p,q) = WORK(2*N+N*NR+NR+p)
1521  774              CONTINUE
1522  773           CONTINUE
1523 *
1524             END IF
1525 *
1526 *           Permute the rows of V using the (column) permutation from the
1527 *           first QRF. Also, scale the columns to make them unit in
1528 *           Euclidean norm. This applies to all cases.
1529 *
1530             TEMP1 = DSQRT(DBLE(N)) * EPSLN
1531             DO 1972 q = 1, N
1532                DO 972 p = 1, N
1533                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1534   972          CONTINUE
1535                DO 973 p = 1, N
1536                   V(p,q) = WORK(2*N+N*NR+NR+p)
1537   973          CONTINUE
1538                XSC = ONE / DNRM2( N, V(1,q), 1 )
1539                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1540      $           CALL DSCAL( N, XSC, V(1,q), 1 )
1541  1972       CONTINUE
1542 *           At this moment, V contains the right singular vectors of A.
1543 *           Next, assemble the left singular vector matrix U (M x N).
1544             IF ( NR .LT. M ) THEN
1545                CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
1546                IF ( NR .LT. N1 ) THEN
1547                   CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
1548                   CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
1549                END IF
1550             END IF
1551 *
1552 *           The Q matrix from the first QRF is built into the left singular
1553 *           matrix U. This applies to all cases.
1554 *
1555             CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U,
1556      $           LDU, WORK(N+1), LWORK-N, IERR )
1557
1558 *           The columns of U are normalized. The cost is O(M*N) flops.
1559             TEMP1 = DSQRT(DBLE(M)) * EPSLN
1560             DO 1973 p = 1, NR
1561                XSC = ONE / DNRM2( M, U(1,p), 1 )
1562                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1563      $          CALL DSCAL( M, XSC, U(1,p), 1 )
1564  1973       CONTINUE
1565 *
1566 *           If the initial QRF is computed with row pivoting, the left
1567 *           singular vectors must be adjusted.
1568 *
1569             IF ( ROWPIV )
1570      $          CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1571 *
1572          ELSE
1573 *
1574 *        .. the initial matrix A has almost orthogonal columns and
1575 *        the second QRF is not needed
1576 *
1577             CALL DLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N )
1578             IF ( L2PERT ) THEN
1579                XSC = DSQRT(SMALL)
1580                DO 5970 p = 2, N
1581                   TEMP1 = XSC * WORK( N + (p-1)*N + p )
1582                   DO 5971 q = 1, p - 1
1583                      WORK(N+(q-1)*N+p)=-DSIGN(TEMP1,WORK(N+(p-1)*N+q))
1584  5971             CONTINUE
1585  5970          CONTINUE
1586             ELSE
1587                CALL DLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N )
1588             END IF
1589 *
1590             CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA,
1591      $           N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO )
1592 *
1593             SCALEM  = WORK(N+N*N+1)
1594             NUMRANK = IDNINT(WORK(N+N*N+2))
1595             DO 6970 p = 1, N
1596                CALL DCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 )
1597                CALL DSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 )
1598  6970       CONTINUE
1599 *
1600             CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N,
1601      $           ONE, A, LDA, WORK(N+1), N )
1602             DO 6972 p = 1, N
1603                CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV )
1604  6972       CONTINUE
1605             TEMP1 = DSQRT(DBLE(N))*EPSLN
1606             DO 6971 p = 1, N
1607                XSC = ONE / DNRM2( N, V(1,p), 1 )
1608                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1609      $            CALL DSCAL( N, XSC, V(1,p), 1 )
1610  6971       CONTINUE
1611 *
1612 *           Assemble the left singular vector matrix U (M x N).
1613 *
1614             IF ( N .LT. M ) THEN
1615                CALL DLASET( 'A',  M-N, N, ZERO, ZERO, U(N+1,1), LDU )
1616                IF ( N .LT. N1 ) THEN
1617                   CALL DLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
1618                   CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU )
1619                END IF
1620             END IF
1621             CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1622      $           LDU, WORK(N+1), LWORK-N, IERR )
1623             TEMP1 = DSQRT(DBLE(M))*EPSLN
1624             DO 6973 p = 1, N1
1625                XSC = ONE / DNRM2( M, U(1,p), 1 )
1626                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1627      $            CALL DSCAL( M, XSC, U(1,p), 1 )
1628  6973       CONTINUE
1629 *
1630             IF ( ROWPIV )
1631      $         CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1632 *
1633          END IF
1634 *
1635 *        end of the  >> almost orthogonal case <<  in the full SVD
1636 *
1637          ELSE
1638 *
1639 *        This branch deploys a preconditioned Jacobi SVD with explicitly
1640 *        accumulated rotations. It is included as optional, mainly for
1641 *        experimental purposes. It does perfom well, and can also be used.
1642 *        In this implementation, this branch will be automatically activated
1643 *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
1644 *        to be greater than the overflow threshold. This is because the
1645 *        a posteriori computation of the singular vectors assumes robust
1646 *        implementation of BLAS and some LAPACK procedures, capable of working
1647 *        in presence of extreme values. Since that is not always the case, ...
1648 *
1649          DO 7968 p = 1, NR
1650             CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1651  7968    CONTINUE
1652 *
1653          IF ( L2PERT ) THEN
1654             XSC = DSQRT(SMALL/EPSLN)
1655             DO 5969 q = 1, NR
1656                TEMP1 = XSC*DABS( V(q,q) )
1657                DO 5968 p = 1, N
1658                   IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
1659      $                .OR. ( p .LT. q ) )
1660      $                V(p,q) = DSIGN( TEMP1, V(p,q) )
1661                   IF ( p .LT. q ) V(p,q) = - V(p,q)
1662  5968          CONTINUE
1663  5969       CONTINUE
1664          ELSE
1665             CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1666          END IF
1667
1668          CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1669      $        LWORK-2*N, IERR )
1670          CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N )
1671 *
1672          DO 7969 p = 1, NR
1673             CALL DCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
1674  7969    CONTINUE
1675
1676          IF ( L2PERT ) THEN
1677             XSC = DSQRT(SMALL/EPSLN)
1678             DO 9970 q = 2, NR
1679                DO 9971 p = 1, q - 1
1680                   TEMP1 = XSC * MIN(DABS(U(p,p)),DABS(U(q,q)))
1681                   U(p,q) = - DSIGN( TEMP1, U(q,p) )
1682  9971          CONTINUE
1683  9970       CONTINUE
1684          ELSE
1685             CALL DLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1686          END IF
1687
1688          CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA,
1689      $        N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO )
1690          SCALEM  = WORK(2*N+N*NR+1)
1691          NUMRANK = IDNINT(WORK(2*N+N*NR+2))
1692
1693          IF ( NR .LT. N ) THEN
1694             CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1695             CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1696             CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1697          END IF
1698
1699          CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1700      $        V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1701 *
1702 *           Permute the rows of V using the (column) permutation from the
1703 *           first QRF. Also, scale the columns to make them unit in
1704 *           Euclidean norm. This applies to all cases.
1705 *
1706             TEMP1 = DSQRT(DBLE(N)) * EPSLN
1707             DO 7972 q = 1, N
1708                DO 8972 p = 1, N
1709                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1710  8972          CONTINUE
1711                DO 8973 p = 1, N
1712                   V(p,q) = WORK(2*N+N*NR+NR+p)
1713  8973          CONTINUE
1714                XSC = ONE / DNRM2( N, V(1,q), 1 )
1715                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1716      $           CALL DSCAL( N, XSC, V(1,q), 1 )
1717  7972       CONTINUE
1718 *
1719 *           At this moment, V contains the right singular vectors of A.
1720 *           Next, assemble the left singular vector matrix U (M x N).
1721 *
1722          IF ( NR .LT. M ) THEN
1723             CALL DLASET( 'A',  M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
1724             IF ( NR .LT. N1 ) THEN
1725                CALL DLASET( 'A',NR,  N1-NR, ZERO, ZERO,  U(1,NR+1),LDU )
1726                CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU )
1727             END IF
1728          END IF
1729 *
1730          CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1731      $        LDU, WORK(N+1), LWORK-N, IERR )
1732 *
1733             IF ( ROWPIV )
1734      $         CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1735 *
1736 *
1737          END IF
1738          IF ( TRANSP ) THEN
1739 *           .. swap U and V because the procedure worked on A^t
1740             DO 6974 p = 1, N
1741                CALL DSWAP( N, U(1,p), 1, V(1,p), 1 )
1742  6974       CONTINUE
1743          END IF
1744 *
1745       END IF
1746 *     end of the full SVD
1747 *
1748 *     Undo scaling, if necessary (and possible)
1749 *
1750       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
1751          CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
1752          USCAL1 = ONE
1753          USCAL2 = ONE
1754       END IF
1755 *
1756       IF ( NR .LT. N ) THEN
1757          DO 3004 p = NR+1, N
1758             SVA(p) = ZERO
1759  3004    CONTINUE
1760       END IF
1761 *
1762       WORK(1) = USCAL2 * SCALEM
1763       WORK(2) = USCAL1
1764       IF ( ERREST ) WORK(3) = SCONDA
1765       IF ( LSVEC .AND. RSVEC ) THEN
1766          WORK(4) = CONDR1
1767          WORK(5) = CONDR2
1768       END IF
1769       IF ( L2TRAN ) THEN
1770          WORK(6) = ENTRA
1771          WORK(7) = ENTRAT
1772       END IF
1773 *
1774       IWORK(1) = NR
1775       IWORK(2) = NUMRANK
1776       IWORK(3) = WARNING
1777 *
1778       RETURN
1779 *     ..
1780 *     .. END OF DGEJSV
1781 *     ..
1782       END
1783 *