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