0001:       SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
0002:      &                   M, N, A, LDA, SVA, U, LDU, V, LDV,
0003:      &                   WORK, LWORK, IWORK, INFO )
0004: *
0005: *  -- LAPACK routine (version 3.2)                                    --
0006: *
0007: *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
0008: *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
0009: *  -- November 2008                                                   --
0010: *
0011: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
0012: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
0013: *
0014: * This routine is also part of SIGMA (version 1.23, October 23. 2008.)
0015: * SIGMA is a library of algorithms for highly accurate algorithms for
0016: * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the
0017: * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
0018: *
0019: *     -#- Scalar Arguments -#-
0020: *
0021:       IMPLICIT    NONE
0022:       INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
0023: *
0024: *     -#- Array Arguments -#-
0025: *
0026:       DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ),
0027:      &            WORK( LWORK )
0028:       INTEGER     IWORK( * )
0029:       CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
0030: *     ..
0031: *
0032: *  Purpose
0033: *  ~~~~~~~
0034: *  DGEJSV computes the singular value decomposition (SVD) of a real M-by-N
0035: *  matrix [A], where M >= N. The SVD of [A] is written as
0036: *
0037: *               [A] = [U] * [SIGMA] * [V]^t,
0038: *
0039: *  where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
0040: *  diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and
0041: *  [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are
0042: *  the singular values of [A]. The columns of [U] and [V] are the left and
0043: *  the right singular vectors of [A], respectively. The matrices [U] and [V]
0044: *  are computed and stored in the arrays U and V, respectively. The diagonal
0045: *  of [SIGMA] is computed and stored in the array SVA.
0046: *
0047: *  Further details
0048: *  ~~~~~~~~~~~~~~~
0049: *  DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3,
0050: *  SGEQRF, and SGELQF as preprocessors and preconditioners. Optionally, an
0051: *  additional row pivoting can be used as a preprocessor, which in some
0052: *  cases results in much higher accuracy. An example is matrix A with the
0053: *  structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned
0054: *  diagonal matrices and C is well-conditioned matrix. In that case, complete
0055: *  pivoting in the first QR factorizations provides accuracy dependent on the
0056: *  condition number of C, and independent of D1, D2. Such higher accuracy is
0057: *  not completely understood theoretically, but it works well in practice.
0058: *  Further, if A can be written as A = B*D, with well-conditioned B and some
0059: *  diagonal D, then the high accuracy is guaranteed, both theoretically and
0060: *  in software, independent of D. For more details see [1], [2].
0061: *     The computational range for the singular values can be the full range
0062: *  ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS
0063: *  & LAPACK routines called by DGEJSV are implemented to work in that range.
0064: *  If that is not the case, then the restriction for safe computation with
0065: *  the singular values in the range of normalized IEEE numbers is that the
0066: *  spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not
0067: *  overflow. This code (DGEJSV) is best used in this restricted range,
0068: *  meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are
0069: *  returned as zeros. See JOBR for details on this.
0070: *     Further, this implementation is somewhat slower than the one described
0071: *  in [1,2] due to replacement of some non-LAPACK components, and because
0072: *  the choice of some tuning parameters in the iterative part (DGESVJ) is
0073: *  left to the implementer on a particular machine.
0074: *     The rank revealing QR factorization (in this code: SGEQP3) should be
0075: *  implemented as in [3]. We have a new version of SGEQP3 under development
0076: *  that is more robust than the current one in LAPACK, with a cleaner cut in
0077: *  rank defficient cases. It will be available in the SIGMA library [4].
0078: *  If M is much larger than N, it is obvious that the inital QRF with
0079: *  column pivoting can be preprocessed by the QRF without pivoting. That
0080: *  well known trick is not used in DGEJSV because in some cases heavy row
0081: *  weighting can be treated with complete pivoting. The overhead in cases
0082: *  M much larger than N is then only due to pivoting, but the benefits in
0083: *  terms of accuracy have prevailed. The implementer/user can incorporate
0084: *  this extra QRF step easily. The implementer can also improve data movement
0085: *  (matrix transpose, matrix copy, matrix transposed copy) - this
0086: *  implementation of DGEJSV uses only the simplest, naive data movement.
0087: *
0088: *  Contributors
0089: *  ~~~~~~~~~~~~
0090: *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
0091: *
0092: *  References
0093: *  ~~~~~~~~~~
0094: * [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
0095: *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
0096: *     LAPACK Working note 169.
0097: * [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
0098: *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
0099: *     LAPACK Working note 170.
0100: * [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
0101: *     factorization software - a case study.
0102: *     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
0103: *     LAPACK Working note 176.
0104: * [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
0105: *     QSVD, (H,K)-SVD computations.
0106: *     Department of Mathematics, University of Zagreb, 2008.
0107: *
0108: *  Bugs, examples and comments
0109: *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
0110: *  Please report all bugs and send interesting examples and/or comments to
0111: *  drmac@math.hr. Thank you.
0112: *
0113: *  Arguments
0114: *  ~~~~~~~~~
0115: *............................................................................
0116: *. JOBA   (input) CHARACTER*1
0117: *.        Specifies the level of accuracy:
0118: *.      = 'C': This option works well (high relative accuracy) if A = B * D,
0119: *.             with well-conditioned B and arbitrary diagonal matrix D.
0120: *.             The accuracy cannot be spoiled by COLUMN scaling. The
0121: *.             accuracy of the computed output depends on the condition of
0122: *.             B, and the procedure aims at the best theoretical accuracy.
0123: *.             The relative error max_{i=1:N}|d sigma_i| / sigma_i is
0124: *.             bounded by f(M,N)*epsilon* cond(B), independent of D.
0125: *.             The input matrix is preprocessed with the QRF with column
0126: *.             pivoting. This initial preprocessing and preconditioning by
0127: *.             a rank revealing QR factorization is common for all values of
0128: *.             JOBA. Additional actions are specified as follows:
0129: *.      = 'E': Computation as with 'C' with an additional estimate of the
0130: *.             condition number of B. It provides a realistic error bound.
0131: *.      = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings
0132: *.             D1, D2, and well-conditioned matrix C, this option gives
0133: *.             higher accuracy than the 'C' option. If the structure of the
0134: *.             input matrix is not known, and relative accuracy is
0135: *.             desirable, then this option is advisable. The input matrix A
0136: *.             is preprocessed with QR factorization with FULL (row and
0137: *.             column) pivoting.
0138: *.      = 'G'  Computation as with 'F' with an additional estimate of the
0139: *.             condition number of B, where A=D*B. If A has heavily weighted
0140: *.             rows, then using this condition number gives too pessimistic
0141: *.             error bound.
0142: *.      = 'A': Small singular values are the noise and the matrix is treated
0143: *.             as numerically rank defficient. The error in the computed
0144: *.             singular values is bounded by f(m,n)*epsilon*||A||.
0145: *.             The computed SVD A = U * S * V^t restores A up to
0146: *.             f(m,n)*epsilon*||A||.
0147: *.             This gives the procedure the licence to discard (set to zero)
0148: *.             all singular values below N*epsilon*||A||.
0149: *.      = 'R': Similar as in 'A'. Rank revealing property of the initial
0150: *.             QR factorization is used do reveal (using triangular factor)
0151: *.             a gap sigma_{r+1} < epsilon * sigma_r in which case the
0152: *.             numerical RANK is declared to be r. The SVD is computed with
0153: *.             absolute error bounds, but more accurately than with 'A'.
0154: *.
0155: *. JOBU   (input) CHARACTER*1
0156: *.        Specifies whether to compute the columns of U:
0157: *.      = 'U': N columns of U are returned in the array U.
0158: *.      = 'F': full set of M left sing. vectors is returned in the array U.
0159: *.      = 'W': U may be used as workspace of length M*N. See the description
0160: *.             of U.
0161: *.      = 'N': U is not computed.
0162: *.
0163: *. JOBV   (input) CHARACTER*1
0164: *.        Specifies whether to compute the matrix V:
0165: *.      = 'V': N columns of V are returned in the array V; Jacobi rotations
0166: *.             are not explicitly accumulated.
0167: *.      = 'J': N columns of V are returned in the array V, but they are
0168: *.             computed as the product of Jacobi rotations. This option is
0169: *.             allowed only if JOBU .NE. 'N', i.e. in computing the full SVD.
0170: *.      = 'W': V may be used as workspace of length N*N. See the description
0171: *.             of V.
0172: *.      = 'N': V is not computed.
0173: *.
0174: *. JOBR   (input) CHARACTER*1
0175: *.        Specifies the RANGE for the singular values. Issues the licence to
0176: *.        set to zero small positive singular values if they are outside
0177: *.        specified range. If A .NE. 0 is scaled so that the largest singular
0178: *.        value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
0179: *.        the licence to kill columns of A whose norm in c*A is less than
0180: *.        DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,
0181: *.        where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
0182: *.      = 'N': Do not kill small columns of c*A. This option assumes that
0183: *.             BLAS and QR factorizations and triangular solvers are
0184: *.             implemented to work in that range. If the condition of A
0185: *.             is greater than BIG, use DGESVJ.
0186: *.      = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)]
0187: *.             (roughly, as described above). This option is recommended.
0188: *.                                            ~~~~~~~~~~~~~~~~~~~~~~~~~~~
0189: *.        For computing the singular values in the FULL range [SFMIN,BIG]
0190: *.        use DGESVJ.
0191: *.
0192: *. JOBT   (input) CHARACTER*1
0193: *.        If the matrix is square then the procedure may determine to use
0194: *.        transposed A if A^t seems to be better with respect to convergence.
0195: *.        If the matrix is not square, JOBT is ignored. This is subject to
0196: *.        changes in the future.
0197: *.        The decision is based on two values of entropy over the adjoint
0198: *.        orbit of A^t * A. See the descriptions of WORK(6) and WORK(7).
0199: *.      = 'T': transpose if entropy test indicates possibly faster
0200: *.        convergence of Jacobi process if A^t is taken as input. If A is
0201: *.        replaced with A^t, then the row pivoting is included automatically.
0202: *.      = 'N': do not speculate.
0203: *.        This option can be used to compute only the singular values, or the
0204: *.        full SVD (U, SIGMA and V). For only one set of singular vectors
0205: *.        (U or V), the caller should provide both U and V, as one of the
0206: *.        matrices is used as workspace if the matrix A is transposed.
0207: *.        The implementer can easily remove this constraint and make the
0208: *.        code more complicated. See the descriptions of U and V.
0209: *.
0210: *. JOBP   (input) CHARACTER*1
0211: *.        Issues the licence to introduce structured perturbations to drown
0212: *.        denormalized numbers. This licence should be active if the
0213: *.        denormals are poorly implemented, causing slow computation,
0214: *.        especially in cases of fast convergence (!). For details see [1,2].
0215: *.        For the sake of simplicity, this perturbations are included only
0216: *.        when the full SVD or only the singular values are requested. The
0217: *.        implementer/user can easily add the perturbation for the cases of
0218: *.        computing one set of singular vectors.
0219: *.      = 'P': introduce perturbation
0220: *.      = 'N': do not perturb
0221: *............................................................................
0222: *
0223: *  M      (input) INTEGER
0224: *         The number of rows of the input matrix A.  M >= 0.
0225: *
0226: *  N      (input) INTEGER
0227: *         The number of columns of the input matrix A. M >= N >= 0.
0228: *
0229: *  A       (input/workspace) REAL array, dimension (LDA,N)
0230: *          On entry, the M-by-N matrix A.
0231: *
0232: *  LDA     (input) INTEGER
0233: *          The leading dimension of the array A.  LDA >= max(1,M).
0234: *
0235: *  SVA     (workspace/output) REAL array, dimension (N)
0236: *          On exit,
0237: *          - For WORK(1)/WORK(2) = ONE: The singular values of A. During the
0238: *            computation SVA contains Euclidean column norms of the
0239: *            iterated matrices in the array A.
0240: *          - For WORK(1) .NE. WORK(2): The singular values of A are
0241: *            (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if
0242: *            sigma_max(A) overflows or if small singular values have been
0243: *            saved from underflow by scaling the input matrix A.
0244: *          - If JOBR='R' then some of the singular values may be returned
0245: *            as exact zeros obtained by "set to zero" because they are
0246: *            below the numerical rank threshold or are denormalized numbers.
0247: *
0248: *  U       (workspace/output) REAL array, dimension ( LDU, N )
0249: *          If JOBU = 'U', then U contains on exit the M-by-N matrix of
0250: *                         the left singular vectors.
0251: *          If JOBU = 'F', then U contains on exit the M-by-M matrix of
0252: *                         the left singular vectors, including an ONB
0253: *                         of the orthogonal complement of the Range(A).
0254: *          If JOBU = 'W'  .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),
0255: *                         then U is used as workspace if the procedure
0256: *                         replaces A with A^t. In that case, [V] is computed
0257: *                         in U as left singular vectors of A^t and then
0258: *                         copied back to the V array. This 'W' option is just
0259: *                         a reminder to the caller that in this case U is
0260: *                         reserved as workspace of length N*N.
0261: *          If JOBU = 'N'  U is not referenced.
0262: *
0263: * LDU      (input) INTEGER
0264: *          The leading dimension of the array U,  LDU >= 1.
0265: *          IF  JOBU = 'U' or 'F' or 'W',  then LDU >= M.
0266: *
0267: *  V       (workspace/output) REAL array, dimension ( LDV, N )
0268: *          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
0269: *                         the right singular vectors;
0270: *          If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
0271: *                         then V is used as workspace if the pprocedure
0272: *                         replaces A with A^t. In that case, [U] is computed
0273: *                         in V as right singular vectors of A^t and then
0274: *                         copied back to the U array. This 'W' option is just
0275: *                         a reminder to the caller that in this case V is
0276: *                         reserved as workspace of length N*N.
0277: *          If JOBV = 'N'  V is not referenced.
0278: *
0279: *  LDV     (input) INTEGER
0280: *          The leading dimension of the array V,  LDV >= 1.
0281: *          If JOBV = 'V' or 'J' or 'W', then LDV >= N.
0282: *
0283: *  WORK    (workspace/output) REAL array, dimension at least LWORK.
0284: *          On exit,
0285: *          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
0286: *                    that SCALE*SVA(1:N) are the computed singular values
0287: *                    of A. (See the description of SVA().)
0288: *          WORK(2) = See the description of WORK(1).
0289: *          WORK(3) = SCONDA is an estimate for the condition number of
0290: *                    column equilibrated A. (If JOBA .EQ. 'E' or 'G')
0291: *                    SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
0292: *                    It is computed using DPOCON. It holds
0293: *                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
0294: *                    where R is the triangular factor from the QRF of A.
0295: *                    However, if R is truncated and the numerical rank is
0296: *                    determined to be strictly smaller than N, SCONDA is
0297: *                    returned as -1, thus indicating that the smallest
0298: *                    singular values might be lost.
0299: *
0300: *          If full SVD is needed, the following two condition numbers are
0301: *          useful for the analysis of the algorithm. They are provied for
0302: *          a developer/implementer who is familiar with the details of
0303: *          the method.
0304: *
0305: *          WORK(4) = an estimate of the scaled condition number of the
0306: *                    triangular factor in the first QR factorization.
0307: *          WORK(5) = an estimate of the scaled condition number of the
0308: *                    triangular factor in the second QR factorization.
0309: *          The following two parameters are computed if JOBT .EQ. 'T'.
0310: *          They are provided for a developer/implementer who is familiar
0311: *          with the details of the method.
0312: *
0313: *          WORK(6) = the entropy of A^t*A :: this is the Shannon entropy
0314: *                    of diag(A^t*A) / Trace(A^t*A) taken as point in the
0315: *                    probability simplex.
0316: *          WORK(7) = the entropy of A*A^t.
0317: *
0318: *  LWORK   (input) INTEGER
0319: *          Length of WORK to confirm proper allocation of work space.
0320: *          LWORK depends on the job:
0321: *
0322: *          If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and
0323: *            -> .. no scaled condition estimate required ( JOBE.EQ.'N'):
0324: *               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
0325: *               For optimal performance (blocked code) the optimal value
0326: *               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
0327: *               block size for xGEQP3/xGEQRF.
0328: *            -> .. an estimate of the scaled condition number of A is
0329: *               required (JOBA='E', 'G'). In this case, LWORK is the maximum
0330: *               of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7).
0331: *
0332: *          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),
0333: *            -> the minimal requirement is LWORK >= max(2*N+M,7).
0334: *            -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7),
0335: *               where NB is the optimal block size.
0336: *
0337: *          If SIGMA and the left singular vectors are needed
0338: *            -> the minimal requirement is LWORK >= max(2*N+M,7).
0339: *            -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7),
0340: *               where NB is the optimal block size.
0341: *
0342: *          If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and
0343: *            -> .. the singular vectors are computed without explicit
0344: *               accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N
0345: *            -> .. in the iterative part, the Jacobi rotations are
0346: *               explicitly accumulated (option, see the description of JOBV),
0347: *               then the minimal requirement is LWORK >= max(M+3*N+N*N,7).
0348: *               For better performance, if NB is the optimal block size,
0349: *               LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7).
0350: *
0351: *  IWORK   (workspace/output) INTEGER array, dimension M+3*N.
0352: *          On exit,
0353: *          IWORK(1) = the numerical rank determined after the initial
0354: *                     QR factorization with pivoting. See the descriptions
0355: *                     of JOBA and JOBR.
0356: *          IWORK(2) = the number of the computed nonzero singular values
0357: *          IWORK(3) = if nonzero, a warning message:
0358: *                     If IWORK(3).EQ.1 then some of the column norms of A
0359: *                     were denormalized floats. The requested high accuracy
0360: *                     is not warranted by the data.
0361: *
0362: *  INFO    (output) INTEGER
0363: *           < 0  : if INFO = -i, then the i-th argument had an illegal value.
0364: *           = 0 :  successfull exit;
0365: *           > 0 :  DGEJSV  did not converge in the maximal allowed number
0366: *                  of sweeps. The computed values may be inaccurate.
0367: *
0368: *............................................................................
0369: *
0370: *     Local Parameters:
0371: *
0372:       DOUBLE PRECISION   ZERO,  ONE
0373:       PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
0374: *
0375: *     Local Scalars:
0376: *
0377:       DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
0378:      &        CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,  MAXPRJ, SCALEM,
0379:      &        SCONDA, SFMIN,  SMALL,  TEMP1,  USCAL1, USCAL2, XSC
0380:       INTEGER IERR,   N1,     NR,     NUMRANK,        p, q,   WARNING
0381:       LOGICAL ALMORT, DEFR,   ERREST, GOSCAL, JRACC,  KILL,   LSVEC,
0382:      &        L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN,
0383:      &        NOSCAL, ROWPIV, RSVEC,  TRANSP
0384: *
0385: *     Intrinsic Functions:
0386: *
0387:       INTRINSIC DABS,  DLOG, DMAX1, DMIN1, DBLE,
0388:      &          MAX0, MIN0, IDNINT,  DSIGN,  DSQRT
0389: *
0390: *     External Functions:
0391: *
0392:       DOUBLE PRECISION  DLAMCH, DNRM2
0393:       INTEGER   IDAMAX
0394:       LOGICAL   LSAME
0395:       EXTERNAL  IDAMAX, LSAME, DLAMCH, DNRM2
0396: *
0397: *     External Subroutines ( BLAS, LAPACK ):
0398: *
0399:       EXTERNAL  DCOPY,  DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL,
0400:      &          DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ,
0401:      &          DORMQR, DPOCON, DSCAL,  DSWAP,  DTRSM,  XERBLA
0402: *
0403:       EXTERNAL  DGESVJ
0404: *
0405: *............................................................................
0406: *
0407: *     Test the input arguments
0408: *
0409:       LSVEC  = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
0410:       JRACC  = LSAME( JOBV, 'J' )
0411:       RSVEC  = LSAME( JOBV, 'V' ) .OR. JRACC
0412:       ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' )
0413:       L2RANK = LSAME( JOBA, 'R' )
0414:       L2ABER = LSAME( JOBA, 'A' )
0415:       ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' )
0416:       L2TRAN = LSAME( JOBT, 'T' )
0417:       L2KILL = LSAME( JOBR, 'R' )
0418:       DEFR   = LSAME( JOBR, 'N' )
0419:       L2PERT = LSAME( JOBP, 'P' )
0420: *
0421:       IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR.
0422:      &     ERREST .OR. LSAME( JOBA, 'C' ) )) THEN
0423:          INFO = - 1
0424:       ELSE IF ( .NOT.( LSVEC  .OR. LSAME( JOBU, 'N' ) .OR.
0425:      &                             LSAME( JOBU, 'W' )) ) THEN
0426:          INFO = - 2
0427:       ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR.
0428:      &   LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN
0429:          INFO = - 3
0430:       ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) )    THEN
0431:          INFO = - 4
0432:       ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN
0433:          INFO = - 5
0434:       ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN
0435:          INFO = - 6
0436:       ELSE IF ( M .LT. 0 ) THEN
0437:          INFO = - 7
0438:       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN
0439:          INFO = - 8
0440:       ELSE IF ( LDA .LT. M ) THEN
0441:          INFO = - 10
0442:       ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN
0443:          INFO = - 13
0444:       ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN
0445:          INFO = - 14
0446:       ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND.
0447:      &                           (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR.
0448:      & (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND.
0449:      &                         (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR.
0450:      & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR.
0451:      & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR.
0452:      & (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N))
0453:      & .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N)))
0454:      &   THEN
0455:          INFO = - 17
0456:       ELSE
0457: *        #:)
0458:          INFO = 0
0459:       END IF
0460: *
0461:       IF ( INFO .NE. 0 ) THEN
0462: *       #:(
0463:          CALL XERBLA( 'DGEJSV', - INFO )
0464:       END IF
0465: *
0466: *     Quick return for void matrix (Y3K safe)
0467: * #:)
0468:       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN
0469: *
0470: *     Determine whether the matrix U should be M x N or M x M
0471: *
0472:       IF ( LSVEC ) THEN
0473:          N1 = N
0474:          IF ( LSAME( JOBU, 'F' ) ) N1 = M
0475:       END IF
0476: *
0477: *     Set numerical parameters
0478: *
0479: *!    NOTE: Make sure DLAMCH() does not fail on the target architecture.
0480: *
0481: 
0482:       EPSLN = DLAMCH('Epsilon')
0483:       SFMIN = DLAMCH('SafeMinimum')
0484:       SMALL = SFMIN / EPSLN
0485:       BIG   = DLAMCH('O')
0486: *     BIG   = ONE / SFMIN
0487: *
0488: *     Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N
0489: *
0490: *(!)  If necessary, scale SVA() to protect the largest norm from
0491: *     overflow. It is possible that this scaling pushes the smallest
0492: *     column norm left from the underflow threshold (extreme case).
0493: *
0494:       SCALEM  = ONE / DSQRT(DBLE(M)*DBLE(N))
0495:       NOSCAL  = .TRUE.
0496:       GOSCAL  = .TRUE.
0497:       DO 1874 p = 1, N
0498:          AAPP = ZERO
0499:          AAQQ = ZERO
0500:          CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ )
0501:          IF ( AAPP .GT. BIG ) THEN
0502:             INFO = - 9
0503:             CALL XERBLA( 'DGEJSV', -INFO )
0504:             RETURN
0505:          END IF
0506:          AAQQ = DSQRT(AAQQ)
0507:          IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL  ) THEN
0508:             SVA(p)  = AAPP * AAQQ
0509:          ELSE
0510:             NOSCAL  = .FALSE.
0511:             SVA(p)  = AAPP * ( AAQQ * SCALEM )
0512:             IF ( GOSCAL ) THEN
0513:                GOSCAL = .FALSE.
0514:                CALL DSCAL( p-1, SCALEM, SVA, 1 )
0515:             END IF
0516:          END IF
0517:  1874 CONTINUE
0518: *
0519:       IF ( NOSCAL ) SCALEM = ONE
0520: *
0521:       AAPP = ZERO
0522:       AAQQ = BIG
0523:       DO 4781 p = 1, N
0524:          AAPP = DMAX1( AAPP, SVA(p) )
0525:          IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) )
0526:  4781 CONTINUE
0527: *
0528: *     Quick return for zero M x N matrix
0529: * #:)
0530:       IF ( AAPP .EQ. ZERO ) THEN
0531:          IF ( LSVEC ) CALL DLASET( 'G', M, N1, ZERO, ONE, U, LDU )
0532:          IF ( RSVEC ) CALL DLASET( 'G', N, N,  ZERO, ONE, V, LDV )
0533:          WORK(1) = ONE
0534:          WORK(2) = ONE
0535:          IF ( ERREST ) WORK(3) = ONE
0536:          IF ( LSVEC .AND. RSVEC ) THEN
0537:             WORK(4) = ONE
0538:             WORK(5) = ONE
0539:          END IF
0540:          IF ( L2TRAN ) THEN
0541:             WORK(6) = ZERO
0542:             WORK(7) = ZERO
0543:          END IF
0544:          IWORK(1) = 0
0545:          IWORK(2) = 0
0546:          RETURN
0547:       END IF
0548: *
0549: *     Issue warning if denormalized column norms detected. Override the
0550: *     high relative accuracy request. Issue licence to kill columns
0551: *     (set them to zero) whose norm is less than sigma_max / BIG (roughly).
0552: * #:(
0553:       WARNING = 0
0554:       IF ( AAQQ .LE. SFMIN ) THEN
0555:          L2RANK = .TRUE.
0556:          L2KILL = .TRUE.
0557:          WARNING = 1
0558:       END IF
0559: *
0560: *     Quick return for one-column matrix
0561: * #:)
0562:       IF ( N .EQ. 1 ) THEN
0563: *
0564:          IF ( LSVEC ) THEN
0565:             CALL DLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR )
0566:             CALL DLACPY( 'A', M, 1, A, LDA, U, LDU )
0567: *           computing all M left singular vectors of the M x 1 matrix
0568:             IF ( N1 .NE. N  ) THEN
0569:                CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR )
0570:                CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR )
0571:                CALL DCOPY( M, A(1,1), 1, U(1,1), 1 )
0572:             END IF
0573:          END IF
0574:          IF ( RSVEC ) THEN
0575:              V(1,1) = ONE
0576:          END IF
0577:          IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN
0578:             SVA(1)  = SVA(1) / SCALEM
0579:             SCALEM  = ONE
0580:          END IF
0581:          WORK(1) = ONE / SCALEM
0582:          WORK(2) = ONE
0583:          IF ( SVA(1) .NE. ZERO ) THEN
0584:             IWORK(1) = 1
0585:             IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN
0586:                IWORK(2) = 1
0587:             ELSE
0588:                IWORK(2) = 0
0589:             END IF
0590:          ELSE
0591:             IWORK(1) = 0
0592:             IWORK(2) = 0
0593:          END IF
0594:          IF ( ERREST ) WORK(3) = ONE
0595:          IF ( LSVEC .AND. RSVEC ) THEN
0596:             WORK(4) = ONE
0597:             WORK(5) = ONE
0598:          END IF
0599:          IF ( L2TRAN ) THEN
0600:             WORK(6) = ZERO
0601:             WORK(7) = ZERO
0602:          END IF
0603:          RETURN
0604: *
0605:       END IF
0606: *
0607:       TRANSP = .FALSE.
0608:       L2TRAN = L2TRAN .AND. ( M .EQ. N )
0609: *
0610:       AATMAX = -ONE
0611:       AATMIN =  BIG
0612:       IF ( ROWPIV .OR. L2TRAN ) THEN
0613: *
0614: *     Compute the row norms, needed to determine row pivoting sequence
0615: *     (in the case of heavily row weighted A, row pivoting is strongly
0616: *     advised) and to collect information needed to compare the
0617: *     structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.).
0618: *
0619:          IF ( L2TRAN ) THEN
0620:             DO 1950 p = 1, M
0621:                XSC   = ZERO
0622:                TEMP1 = ZERO
0623:                CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 )
0624: *              DLASSQ gets both the ell_2 and the ell_infinity norm
0625: *              in one pass through the vector
0626:                WORK(M+N+p)  = XSC * SCALEM
0627:                WORK(N+p)    = XSC * (SCALEM*DSQRT(TEMP1))
0628:                AATMAX = DMAX1( AATMAX, WORK(N+p) )
0629:                IF (WORK(N+p) .NE. ZERO) AATMIN = DMIN1(AATMIN,WORK(N+p))
0630:  1950       CONTINUE
0631:          ELSE
0632:             DO 1904 p = 1, M
0633:                WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) )
0634:                AATMAX = DMAX1( AATMAX, WORK(M+N+p) )
0635:                AATMIN = DMIN1( AATMIN, WORK(M+N+p) )
0636:  1904       CONTINUE
0637:          END IF
0638: *
0639:       END IF
0640: *
0641: *     For square matrix A try to determine whether A^t  would be  better
0642: *     input for the preconditioned Jacobi SVD, with faster convergence.
0643: *     The decision is based on an O(N) function of the vector of column
0644: *     and row norms of A, based on the Shannon entropy. This should give
0645: *     the right choice in most cases when the difference actually matters.
0646: *     It may fail and pick the slower converging side.
0647: *
0648:       ENTRA  = ZERO
0649:       ENTRAT = ZERO
0650:       IF ( L2TRAN ) THEN
0651: *
0652:          XSC   = ZERO
0653:          TEMP1 = ZERO
0654:          CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
0655:          TEMP1 = ONE / TEMP1
0656: *
0657:          ENTRA = ZERO
0658:          DO 1113 p = 1, N
0659:             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
0660:             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
0661:  1113    CONTINUE
0662:          ENTRA = - ENTRA / DLOG(DBLE(N))
0663: *
0664: *        Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex.
0665: *        It is derived from the diagonal of  A^t * A.  Do the same with the
0666: *        diagonal of A * A^t, compute the entropy of the corresponding
0667: *        probability distribution. Note that A * A^t and A^t * A have the
0668: *        same trace.
0669: *
0670:          ENTRAT = ZERO
0671:          DO 1114 p = N+1, N+M
0672:             BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1
0673:             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
0674:  1114    CONTINUE
0675:          ENTRAT = - ENTRAT / DLOG(DBLE(M))
0676: *
0677: *        Analyze the entropies and decide A or A^t. Smaller entropy
0678: *        usually means better input for the algorithm.
0679: *
0680:          TRANSP = ( ENTRAT .LT. ENTRA )
0681: *
0682: *        If A^t is better than A, transpose A.
0683: *
0684:          IF ( TRANSP ) THEN
0685: *           In an optimal implementation, this trivial transpose
0686: *           should be replaced with faster transpose.
0687:             DO 1115 p = 1, N - 1
0688:                DO 1116 q = p + 1, N
0689:                    TEMP1 = A(q,p)
0690:                   A(q,p) = A(p,q)
0691:                   A(p,q) = TEMP1
0692:  1116          CONTINUE
0693:  1115       CONTINUE
0694:             DO 1117 p = 1, N
0695:                WORK(M+N+p) = SVA(p)
0696:                SVA(p)      = WORK(N+p)
0697:  1117       CONTINUE
0698:             TEMP1  = AAPP
0699:             AAPP   = AATMAX
0700:             AATMAX = TEMP1
0701:             TEMP1  = AAQQ
0702:             AAQQ   = AATMIN
0703:             AATMIN = TEMP1
0704:             KILL   = LSVEC
0705:             LSVEC  = RSVEC
0706:             RSVEC  = KILL
0707: *
0708:             ROWPIV = .TRUE.
0709:          END IF
0710: *
0711:       END IF
0712: *     END IF L2TRAN
0713: *
0714: *     Scale the matrix so that its maximal singular value remains less
0715: *     than DSQRT(BIG) -- the matrix is scaled so that its maximal column
0716: *     has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep
0717: *     DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and
0718: *     BLAS routines that, in some implementations, are not capable of
0719: *     working in the full interval [SFMIN,BIG] and that they may provoke
0720: *     overflows in the intermediate results. If the singular values spread
0721: *     from SFMIN to BIG, then DGESVJ will compute them. So, in that case,
0722: *     one should use DGESVJ instead of DGEJSV.
0723: *
0724:       BIG1   = DSQRT( BIG )
0725:       TEMP1  = DSQRT( BIG / DBLE(N) )
0726: *
0727:       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
0728:       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
0729:           AAQQ = ( AAQQ / AAPP ) * TEMP1
0730:       ELSE
0731:           AAQQ = ( AAQQ * TEMP1 ) / AAPP
0732:       END IF
0733:       TEMP1 = TEMP1 * SCALEM
0734:       CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR )
0735: *
0736: *     To undo scaling at the end of this procedure, multiply the
0737: *     computed singular values with USCAL2 / USCAL1.
0738: *
0739:       USCAL1 = TEMP1
0740:       USCAL2 = AAPP
0741: *
0742:       IF ( L2KILL ) THEN
0743: *        L2KILL enforces computation of nonzero singular values in
0744: *        the restricted range of condition number of the initial A,
0745: *        sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN).
0746:          XSC = DSQRT( SFMIN )
0747:       ELSE
0748:          XSC = SMALL
0749: *
0750: *        Now, if the condition number of A is too big,
0751: *        sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN,
0752: *        as a precaution measure, the full SVD is computed using DGESVJ
0753: *        with accumulated Jacobi rotations. This provides numerically
0754: *        more robust computation, at the cost of slightly increased run
0755: *        time. Depending on the concrete implementation of BLAS and LAPACK
0756: *        (i.e. how they behave in presence of extreme ill-conditioning) the
0757: *        implementor may decide to remove this switch.
0758:          IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN
0759:             JRACC = .TRUE.
0760:          END IF
0761: *
0762:       END IF
0763:       IF ( AAQQ .LT. XSC ) THEN
0764:          DO 700 p = 1, N
0765:             IF ( SVA(p) .LT. XSC ) THEN
0766:                CALL DLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA )
0767:                SVA(p) = ZERO
0768:             END IF
0769:  700     CONTINUE
0770:       END IF
0771: *
0772: *     Preconditioning using QR factorization with pivoting
0773: *
0774:       IF ( ROWPIV ) THEN
0775: *        Optional row permutation (Bjoerck row pivoting):
0776: *        A result by Cox and Higham shows that the Bjoerck's
0777: *        row pivoting combined with standard column pivoting
0778: *        has similar effect as Powell-Reid complete pivoting.
0779: *        The ell-infinity norms of A are made nonincreasing.
0780:          DO 1952 p = 1, M - 1
0781:             q = IDAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1
0782:             IWORK(2*N+p) = q
0783:             IF ( p .NE. q ) THEN
0784:                TEMP1       = WORK(M+N+p)
0785:                WORK(M+N+p) = WORK(M+N+q)
0786:                WORK(M+N+q) = TEMP1
0787:             END IF
0788:  1952    CONTINUE
0789:          CALL DLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 )
0790:       END IF
0791: *
0792: *     End of the preparation phase (scaling, optional sorting and
0793: *     transposing, optional flushing of small columns).
0794: *
0795: *     Preconditioning
0796: *
0797: *     If the full SVD is needed, the right singular vectors are computed
0798: *     from a matrix equation, and for that we need theoretical analysis
0799: *     of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF.
0800: *     In all other cases the first RR QRF can be chosen by other criteria
0801: *     (eg speed by replacing global with restricted window pivoting, such
0802: *     as in SGEQPX from TOMS # 782). Good results will be obtained using
0803: *     SGEQPX with properly (!) chosen numerical parameters.
0804: *     Any improvement of DGEQP3 improves overal performance of DGEJSV.
0805: *
0806: *     A * P1 = Q1 * [ R1^t 0]^t:
0807:       DO 1963 p = 1, N
0808: *        .. all columns are free columns
0809:          IWORK(p) = 0
0810:  1963 CONTINUE
0811:       CALL DGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR )
0812: *
0813: *     The upper triangular matrix R1 from the first QRF is inspected for
0814: *     rank deficiency and possibilities for deflation, or possible
0815: *     ill-conditioning. Depending on the user specified flag L2RANK,
0816: *     the procedure explores possibilities to reduce the numerical
0817: *     rank by inspecting the computed upper triangular factor. If
0818: *     L2RANK or L2ABER are up, then DGEJSV will compute the SVD of
0819: *     A + dA, where ||dA|| <= f(M,N)*EPSLN.
0820: *
0821:       NR = 1
0822:       IF ( L2ABER ) THEN
0823: *        Standard absolute error bound suffices. All sigma_i with
0824: *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
0825: *        agressive enforcement of lower numerical rank by introducing a
0826: *        backward error of the order of N*EPSLN*||A||.
0827:          TEMP1 = DSQRT(DBLE(N))*EPSLN
0828:          DO 3001 p = 2, N
0829:             IF ( DABS(A(p,p)) .GE. (TEMP1*DABS(A(1,1))) ) THEN
0830:                NR = NR + 1
0831:             ELSE
0832:                GO TO 3002
0833:             END IF
0834:  3001    CONTINUE
0835:  3002    CONTINUE
0836:       ELSE IF ( L2RANK ) THEN
0837: *        .. similarly as above, only slightly more gentle (less agressive).
0838: *        Sudden drop on the diagonal of R1 is used as the criterion for
0839: *        close-to-rank-defficient.
0840:          TEMP1 = DSQRT(SFMIN)
0841:          DO 3401 p = 2, N
0842:             IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR.
0843:      &           ( DABS(A(p,p)) .LT. SMALL ) .OR.
0844:      &           ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402
0845:             NR = NR + 1
0846:  3401    CONTINUE
0847:  3402    CONTINUE
0848: *
0849:       ELSE
0850: *        The goal is high relative accuracy. However, if the matrix
0851: *        has high scaled condition number the relative accuracy is in
0852: *        general not feasible. Later on, a condition number estimator
0853: *        will be deployed to estimate the scaled condition number.
0854: *        Here we just remove the underflowed part of the triangular
0855: *        factor. This prevents the situation in which the code is
0856: *        working hard to get the accuracy not warranted by the data.
0857:          TEMP1  = DSQRT(SFMIN)
0858:          DO 3301 p = 2, N
0859:             IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR.
0860:      &          ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302
0861:             NR = NR + 1
0862:  3301    CONTINUE
0863:  3302    CONTINUE
0864: *
0865:       END IF
0866: *
0867:       ALMORT = .FALSE.
0868:       IF ( NR .EQ. N ) THEN
0869:          MAXPRJ = ONE
0870:          DO 3051 p = 2, N
0871:             TEMP1  = DABS(A(p,p)) / SVA(IWORK(p))
0872:             MAXPRJ = DMIN1( MAXPRJ, TEMP1 )
0873:  3051    CONTINUE
0874:          IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
0875:       END IF
0876: *
0877: *
0878:       SCONDA = - ONE
0879:       CONDR1 = - ONE
0880:       CONDR2 = - ONE
0881: *
0882:       IF ( ERREST ) THEN
0883:          IF ( N .EQ. NR ) THEN
0884:             IF ( RSVEC ) THEN
0885: *              .. V is available as workspace
0886:                CALL DLACPY( 'U', N, N, A, LDA, V, LDV )
0887:                DO 3053 p = 1, N
0888:                   TEMP1 = SVA(IWORK(p))
0889:                   CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 )
0890:  3053          CONTINUE
0891:                CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1,
0892:      &              WORK(N+1), IWORK(2*N+M+1), IERR )
0893:             ELSE IF ( LSVEC ) THEN
0894: *              .. U is available as workspace
0895:                CALL DLACPY( 'U', N, N, A, LDA, U, LDU )
0896:                DO 3054 p = 1, N
0897:                   TEMP1 = SVA(IWORK(p))
0898:                   CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 )
0899:  3054          CONTINUE
0900:                CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1,
0901:      &              WORK(N+1), IWORK(2*N+M+1), IERR )
0902:             ELSE
0903:                CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N )
0904:                DO 3052 p = 1, N
0905:                   TEMP1 = SVA(IWORK(p))
0906:                   CALL DSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 )
0907:  3052          CONTINUE
0908: *           .. the columns of R are scaled to have unit Euclidean lengths.
0909:                CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1,
0910:      &              WORK(N+N*N+1), IWORK(2*N+M+1), IERR )
0911:             END IF
0912:             SCONDA = ONE / DSQRT(TEMP1)
0913: *           SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
0914: *           N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
0915:          ELSE
0916:             SCONDA = - ONE
0917:          END IF
0918:       END IF
0919: *
0920:       L2PERT = L2PERT .AND. ( DABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) )
0921: *     If there is no violent scaling, artificial perturbation is not needed.
0922: *
0923: *     Phase 3:
0924: *
0925: 
0926:       IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
0927: *
0928: *         Singular Values only
0929: *
0930: *         .. transpose A(1:NR,1:N)
0931:          DO 1946 p = 1, MIN0( N-1, NR )
0932:             CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 )
0933:  1946    CONTINUE
0934: *
0935: *        The following two DO-loops introduce small relative perturbation
0936: *        into the strict upper triangle of the lower triangular matrix.
0937: *        Small entries below the main diagonal are also changed.
0938: *        This modification is useful if the computing environment does not
0939: *        provide/allow FLUSH TO ZERO underflow, for it prevents many
0940: *        annoying denormalized numbers in case of strongly scaled matrices.
0941: *        The perturbation is structured so that it does not introduce any
0942: *        new perturbation of the singular values, and it does not destroy
0943: *        the job done by the preconditioner.
0944: *        The licence for this perturbation is in the variable L2PERT, which
0945: *        should be .FALSE. if FLUSH TO ZERO underflow is active.
0946: *
0947:          IF ( .NOT. ALMORT ) THEN
0948: *
0949:             IF ( L2PERT ) THEN
0950: *              XSC = DSQRT(SMALL)
0951:                XSC = EPSLN / DBLE(N)
0952:                DO 4947 q = 1, NR
0953:                   TEMP1 = XSC*DABS(A(q,q))
0954:                   DO 4949 p = 1, N
0955:                      IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
0956:      &                    .OR. ( p .LT. q ) )
0957:      &                     A(p,q) = DSIGN( TEMP1, A(p,q) )
0958:  4949             CONTINUE
0959:  4947          CONTINUE
0960:             ELSE
0961:                CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA )
0962:             END IF
0963: *
0964: *            .. second preconditioning using the QR factorization
0965: *
0966:             CALL DGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR )
0967: *
0968: *           .. and transpose upper to lower triangular
0969:             DO 1948 p = 1, NR - 1
0970:                CALL DCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 )
0971:  1948       CONTINUE
0972: *
0973:          END IF
0974: *
0975: *           Row-cyclic Jacobi SVD algorithm with column pivoting
0976: *
0977: *           .. again some perturbation (a "background noise") is added
0978: *           to drown denormals
0979:             IF ( L2PERT ) THEN
0980: *              XSC = DSQRT(SMALL)
0981:                XSC = EPSLN / DBLE(N)
0982:                DO 1947 q = 1, NR
0983:                   TEMP1 = XSC*DABS(A(q,q))
0984:                   DO 1949 p = 1, NR
0985:                      IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) )
0986:      &                       .OR. ( p .LT. q ) )
0987:      &                   A(p,q) = DSIGN( TEMP1, A(p,q) )
0988:  1949             CONTINUE
0989:  1947          CONTINUE
0990:             ELSE
0991:                CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA )
0992:             END IF
0993: *
0994: *           .. and one-sided Jacobi rotations are started on a lower
0995: *           triangular matrix (plus perturbation which is ignored in
0996: *           the part which destroys triangular form (confusing?!))
0997: *
0998:             CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA,
0999:      &                      N, V, LDV, WORK, LWORK, INFO )
1000: *
1001:             SCALEM  = WORK(1)
1002:             NUMRANK = IDNINT(WORK(2))
1003: *
1004: *
1005:       ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
1006: *
1007: *        -> Singular Values and Right Singular Vectors <-
1008: *
1009:          IF ( ALMORT ) THEN
1010: *
1011: *           .. in this case NR equals N
1012:             DO 1998 p = 1, NR
1013:                CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1014:  1998       CONTINUE
1015:             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1016: *
1017:             CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA,
1018:      &                  WORK, LWORK, INFO )
1019:             SCALEM  = WORK(1)
1020:             NUMRANK = IDNINT(WORK(2))
1021: 
1022:          ELSE
1023: *
1024: *        .. two more QR factorizations ( one QRF is not enough, two require
1025: *        accumulated product of Jacobi rotations, three are perfect )
1026: *
1027:             CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA )
1028:             CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR)
1029:             CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV )
1030:             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1031:             CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1032:      &                   LWORK-2*N, IERR )
1033:             DO 8998 p = 1, NR
1034:                CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
1035:  8998       CONTINUE
1036:             CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1037: *
1038:             CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
1039:      &                  LDU, WORK(N+1), LWORK, INFO )
1040:             SCALEM  = WORK(N+1)
1041:             NUMRANK = IDNINT(WORK(N+2))
1042:             IF ( NR .LT. N ) THEN
1043:                CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1),   LDV )
1044:                CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1),   LDV )
1045:                CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV )
1046:             END IF
1047: *
1048:          CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK,
1049:      &               V, LDV, WORK(N+1), LWORK-N, IERR )
1050: *
1051:          END IF
1052: *
1053:          DO 8991 p = 1, N
1054:             CALL DCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA )
1055:  8991    CONTINUE
1056:          CALL DLACPY( 'All', N, N, A, LDA, V, LDV )
1057: *
1058:          IF ( TRANSP ) THEN
1059:             CALL DLACPY( 'All', N, N, V, LDV, U, LDU )
1060:          END IF
1061: *
1062:       ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN
1063: *
1064: *        -#- Singular Values and Left Singular Vectors                 -#-
1065: *
1066: *        .. second preconditioning step to avoid need to accumulate
1067: *        Jacobi rotations in the Jacobi iterations.
1068:          DO 1965 p = 1, NR
1069:             CALL DCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 )
1070:  1965    CONTINUE
1071:          CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1072: *
1073:          CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1),
1074:      &              LWORK-2*N, IERR )
1075: *
1076:          DO 1967 p = 1, NR - 1
1077:             CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 )
1078:  1967    CONTINUE
1079:          CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1080: *
1081:          CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A,
1082:      &        LDA, WORK(N+1), LWORK-N, INFO )
1083:          SCALEM  = WORK(N+1)
1084:          NUMRANK = IDNINT(WORK(N+2))
1085: *
1086:          IF ( NR .LT. M ) THEN
1087:             CALL DLASET( 'A',  M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU )
1088:             IF ( NR .LT. N1 ) THEN
1089:                CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU )
1090:                CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU )
1091:             END IF
1092:          END IF
1093: *
1094:          CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1095:      &               LDU, WORK(N+1), LWORK-N, IERR )
1096: *
1097:          IF ( ROWPIV )
1098:      &       CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1099: *
1100:          DO 1974 p = 1, N1
1101:             XSC = ONE / DNRM2( M, U(1,p), 1 )
1102:             CALL DSCAL( M, XSC, U(1,p), 1 )
1103:  1974    CONTINUE
1104: *
1105:          IF ( TRANSP ) THEN
1106:             CALL DLACPY( 'All', N, N, U, LDU, V, LDV )
1107:          END IF
1108: *
1109:       ELSE
1110: *
1111: *        -#- Full SVD -#-
1112: *
1113:          IF ( .NOT. JRACC ) THEN
1114: *
1115:          IF ( .NOT. ALMORT ) THEN
1116: *
1117: *           Second Preconditioning Step (QRF [with pivoting])
1118: *           Note that the composition of TRANSPOSE, QRF and TRANSPOSE is
1119: *           equivalent to an LQF CALL. Since in many libraries the QRF
1120: *           seems to be better optimized than the LQF, we do explicit
1121: *           transpose and use the QRF. This is subject to changes in an
1122: *           optimized implementation of DGEJSV.
1123: *
1124:             DO 1968 p = 1, NR
1125:                CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1126:  1968       CONTINUE
1127: *
1128: *           .. the following two loops perturb small entries to avoid
1129: *           denormals in the second QR factorization, where they are
1130: *           as good as zeros. This is done to avoid painfully slow
1131: *           computation with denormals. The relative size of the perturbation
1132: *           is a parameter that can be changed by the implementer.
1133: *           This perturbation device will be obsolete on machines with
1134: *           properly implemented arithmetic.
1135: *           To switch it off, set L2PERT=.FALSE. To remove it from  the
1136: *           code, remove the action under L2PERT=.TRUE., leave the ELSE part.
1137: *           The following two loops should be blocked and fused with the
1138: *           transposed copy above.
1139: *
1140:             IF ( L2PERT ) THEN
1141:                XSC = DSQRT(SMALL)
1142:                DO 2969 q = 1, NR
1143:                   TEMP1 = XSC*DABS( V(q,q) )
1144:                   DO 2968 p = 1, N
1145:                      IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
1146:      &                   .OR. ( p .LT. q ) )
1147:      &                   V(p,q) = DSIGN( TEMP1, V(p,q) )
1148:                      IF ( p. LT. q ) V(p,q) = - V(p,q)
1149:  2968             CONTINUE
1150:  2969          CONTINUE
1151:             ELSE
1152:                CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1153:             END IF
1154: *
1155: *           Estimate the row scaled condition number of R1
1156: *           (If R1 is rectangular, N > NR, then the condition number
1157: *           of the leading NR x NR submatrix is estimated.)
1158: *
1159:             CALL DLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR )
1160:             DO 3950 p = 1, NR
1161:                TEMP1 = DNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1)
1162:                CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1)
1163:  3950       CONTINUE
1164:             CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
1165:      &                   WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
1166:             CONDR1 = ONE / DSQRT(TEMP1)
1167: *           .. here need a second oppinion on the condition number
1168: *           .. then assume worst case scenario
1169: *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
1170: *           more conservative    <=> CONDR1 .LT. DSQRT(DBLE(N))
1171: *
1172:             COND_OK = DSQRT(DBLE(NR))
1173: *[TP]       COND_OK is a tuning parameter.
1174: 
1175:             IF ( CONDR1 .LT. COND_OK ) THEN
1176: *              .. the second QRF without pivoting. Note: in an optimized
1177: *              implementation, this QRF should be implemented as the QRF
1178: *              of a lower triangular matrix.
1179: *              R1^t = Q2 * R2
1180:                CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1181:      &              LWORK-2*N, IERR )
1182: *
1183:                IF ( L2PERT ) THEN
1184:                   XSC = DSQRT(SMALL)/EPSLN
1185:                   DO 3959 p = 2, NR
1186:                      DO 3958 q = 1, p - 1
1187:                         TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
1188:                         IF ( DABS(V(q,p)) .LE. TEMP1 )
1189:      &                     V(q,p) = DSIGN( TEMP1, V(q,p) )
1190:  3958                CONTINUE
1191:  3959             CONTINUE
1192:                END IF
1193: *
1194:                IF ( NR .NE. N )
1195: *              .. save ...
1196:      &         CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
1197: *
1198: *           .. this transposed copy should be better than naive
1199:                DO 1969 p = 1, NR - 1
1200:                   CALL DCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 )
1201:  1969          CONTINUE
1202: *
1203:                CONDR2 = CONDR1
1204: *
1205:             ELSE
1206: *
1207: *              .. ill-conditioned case: second QRF with pivoting
1208: *              Note that windowed pivoting would be equaly good
1209: *              numerically, and more run-time efficient. So, in
1210: *              an optimal implementation, the next call to DGEQP3
1211: *              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
1212: *              with properly (carefully) chosen parameters.
1213: *
1214: *              R1^t * P2 = Q2 * R2
1215:                DO 3003 p = 1, NR
1216:                   IWORK(N+p) = 0
1217:  3003          CONTINUE
1218:                CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1),
1219:      &                  WORK(2*N+1), LWORK-2*N, IERR )
1220: **               CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1221: **     &              LWORK-2*N, IERR )
1222:                IF ( L2PERT ) THEN
1223:                   XSC = DSQRT(SMALL)
1224:                   DO 3969 p = 2, NR
1225:                      DO 3968 q = 1, p - 1
1226:                         TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
1227:                         IF ( DABS(V(q,p)) .LE. TEMP1 )
1228:      &                     V(q,p) = DSIGN( TEMP1, V(q,p) )
1229:  3968                CONTINUE
1230:  3969             CONTINUE
1231:                END IF
1232: *
1233:                CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N )
1234: *
1235:                IF ( L2PERT ) THEN
1236:                   XSC = DSQRT(SMALL)
1237:                   DO 8970 p = 2, NR
1238:                      DO 8971 q = 1, p - 1
1239:                         TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q)))
1240:                         V(p,q) = - DSIGN( TEMP1, V(q,p) )
1241:  8971                CONTINUE
1242:  8970             CONTINUE
1243:                ELSE
1244:                   CALL DLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV )
1245:                END IF
1246: *              Now, compute R2 = L3 * Q3, the LQ factorization.
1247:                CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1),
1248:      &               WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR )
1249: *              .. and estimate the condition number
1250:                CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR )
1251:                DO 4950 p = 1, NR
1252:                   TEMP1 = DNRM2( p, WORK(2*N+N*NR+NR+p), NR )
1253:                   CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR )
1254:  4950          CONTINUE
1255:                CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1,
1256:      &              WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR )
1257:                CONDR2 = ONE / DSQRT(TEMP1)
1258: *
1259:                IF ( CONDR2 .GE. COND_OK ) THEN
1260: *                 .. save the Householder vectors used for Q3
1261: *                 (this overwrittes the copy of R2, as it will not be
1262: *                 needed in this branch, but it does not overwritte the
1263: *                 Huseholder vectors of Q2.).
1264:                   CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
1265: *                 .. and the rest of the information on Q3 is in
1266: *                 WORK(2*N+N*NR+1:2*N+N*NR+N)
1267:                END IF
1268: *
1269:             END IF
1270: *
1271:             IF ( L2PERT ) THEN
1272:                XSC = DSQRT(SMALL)
1273:                DO 4968 q = 2, NR
1274:                   TEMP1 = XSC * V(q,q)
1275:                   DO 4969 p = 1, q - 1
1276: *                    V(p,q) = - DSIGN( TEMP1, V(q,p) )
1277:                      V(p,q) = - DSIGN( TEMP1, V(p,q) )
1278:  4969             CONTINUE
1279:  4968          CONTINUE
1280:             ELSE
1281:                CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV )
1282:             END IF
1283: *
1284: *        Second preconditioning finished; continue with Jacobi SVD
1285: *        The input matrix is lower trinagular.
1286: *
1287: *        Recover the right singular vectors as solution of a well
1288: *        conditioned triangular matrix equation.
1289: *
1290:             IF ( CONDR1 .LT. COND_OK ) THEN
1291: *
1292:                CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U,
1293:      &              LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO )
1294:                SCALEM  = WORK(2*N+N*NR+NR+1)
1295:                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1296:                DO 3970 p = 1, NR
1297:                   CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
1298:                   CALL DSCAL( NR, SVA(p),    V(1,p), 1 )
1299:  3970          CONTINUE
1300: 
1301: *        .. pick the right matrix equation and solve it
1302: *
1303:                IF ( NR. EQ. N ) THEN
1304: * :))             .. best case, R1 is inverted. The solution of this matrix
1305: *                 equation is Q2*V2 = the product of the Jacobi rotations
1306: *                 used in DGESVJ, premultiplied with the orthogonal matrix
1307: *                 from the second QR factorization.
1308:                   CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV )
1309:                ELSE
1310: *                 .. R1 is well conditioned, but non-square. Transpose(R2)
1311: *                 is inverted to get the product of the Jacobi rotations
1312: *                 used in DGESVJ. The Q-factor from the second QR
1313: *                 factorization is then built in explicitly.
1314:                   CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1),
1315:      &                 N,V,LDV)
1316:                   IF ( NR .LT. N ) THEN
1317:                     CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV)
1318:                     CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV)
1319:                     CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV)
1320:                   END IF
1321:                   CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1322:      &                 V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
1323:                END IF
1324: *
1325:             ELSE IF ( CONDR2 .LT. COND_OK ) THEN
1326: *
1327: * :)           .. the input matrix A is very likely a relative of
1328: *              the Kahan matrix :)
1329: *              The matrix R2 is inverted. The solution of the matrix equation
1330: *              is Q3^T*V3 = the product of the Jacobi rotations (appplied to
1331: *              the lower triangular L3 from the LQ factorization of
1332: *              R2=L3*Q3), pre-multiplied with the transposed Q3.
1333:                CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U,
1334:      &              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
1335:                SCALEM  = WORK(2*N+N*NR+NR+1)
1336:                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1337:                DO 3870 p = 1, NR
1338:                   CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 )
1339:                   CALL DSCAL( NR, SVA(p),    U(1,p), 1 )
1340:  3870          CONTINUE
1341:                CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU)
1342: *              .. apply the permutation from the second QR factorization
1343:                DO 873 q = 1, NR
1344:                   DO 872 p = 1, NR
1345:                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1346:  872              CONTINUE
1347:                   DO 874 p = 1, NR
1348:                      U(p,q) = WORK(2*N+N*NR+NR+p)
1349:  874              CONTINUE
1350:  873           CONTINUE
1351:                IF ( NR .LT. N ) THEN
1352:                   CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1353:                   CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1354:                   CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1355:                END IF
1356:                CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1357:      &              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1358:             ELSE
1359: *              Last line of defense.
1360: * #:(          This is a rather pathological case: no scaled condition
1361: *              improvement after two pivoted QR factorizations. Other
1362: *              possibility is that the rank revealing QR factorization
1363: *              or the condition estimator has failed, or the COND_OK
1364: *              is set very close to ONE (which is unnecessary). Normally,
1365: *              this branch should never be executed, but in rare cases of
1366: *              failure of the RRQR or condition estimator, the last line of
1367: *              defense ensures that DGEJSV completes the task.
1368: *              Compute the full SVD of L3 using DGESVJ with explicit
1369: *              accumulation of Jacobi rotations.
1370:                CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U,
1371:      &              LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO )
1372:                SCALEM  = WORK(2*N+N*NR+NR+1)
1373:                NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2))
1374:                IF ( NR .LT. N ) THEN
1375:                   CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1376:                   CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1377:                   CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1378:                END IF
1379:                CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1380:      &              V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1381: *
1382:                CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N,
1383:      &              WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1),
1384:      &              LWORK-2*N-N*NR-NR, IERR )
1385:                DO 773 q = 1, NR
1386:                   DO 772 p = 1, NR
1387:                      WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q)
1388:  772              CONTINUE
1389:                   DO 774 p = 1, NR
1390:                      U(p,q) = WORK(2*N+N*NR+NR+p)
1391:  774              CONTINUE
1392:  773           CONTINUE
1393: *
1394:             END IF
1395: *
1396: *           Permute the rows of V using the (column) permutation from the
1397: *           first QRF. Also, scale the columns to make them unit in
1398: *           Euclidean norm. This applies to all cases.
1399: *
1400:             TEMP1 = DSQRT(DBLE(N)) * EPSLN
1401:             DO 1972 q = 1, N
1402:                DO 972 p = 1, N
1403:                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1404:   972          CONTINUE
1405:                DO 973 p = 1, N
1406:                   V(p,q) = WORK(2*N+N*NR+NR+p)
1407:   973          CONTINUE
1408:                XSC = ONE / DNRM2( N, V(1,q), 1 )
1409:                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1410:      &           CALL DSCAL( N, XSC, V(1,q), 1 )
1411:  1972       CONTINUE
1412: *           At this moment, V contains the right singular vectors of A.
1413: *           Next, assemble the left singular vector matrix U (M x N).
1414:             IF ( NR .LT. M ) THEN
1415:                CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU )
1416:                IF ( NR .LT. N1 ) THEN
1417:                   CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU)
1418:                   CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU)
1419:                END IF
1420:             END IF
1421: *
1422: *           The Q matrix from the first QRF is built into the left singular
1423: *           matrix U. This applies to all cases.
1424: *
1425:             CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U,
1426:      &           LDU, WORK(N+1), LWORK-N, IERR )
1427: 
1428: *           The columns of U are normalized. The cost is O(M*N) flops.
1429:             TEMP1 = DSQRT(DBLE(M)) * EPSLN
1430:             DO 1973 p = 1, NR
1431:                XSC = ONE / DNRM2( M, U(1,p), 1 )
1432:                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1433:      &          CALL DSCAL( M, XSC, U(1,p), 1 )
1434:  1973       CONTINUE
1435: *
1436: *           If the initial QRF is computed with row pivoting, the left
1437: *           singular vectors must be adjusted.
1438: *
1439:             IF ( ROWPIV )
1440:      &          CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1441: *
1442:          ELSE
1443: *
1444: *        .. the initial matrix A has almost orthogonal columns and
1445: *        the second QRF is not needed
1446: *
1447:             CALL DLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N )
1448:             IF ( L2PERT ) THEN
1449:                XSC = DSQRT(SMALL)
1450:                DO 5970 p = 2, N
1451:                   TEMP1 = XSC * WORK( N + (p-1)*N + p )
1452:                   DO 5971 q = 1, p - 1
1453:                      WORK(N+(q-1)*N+p)=-DSIGN(TEMP1,WORK(N+(p-1)*N+q))
1454:  5971             CONTINUE
1455:  5970          CONTINUE
1456:             ELSE
1457:                CALL DLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N )
1458:             END IF
1459: *
1460:             CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA,
1461:      &           N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO )
1462: *
1463:             SCALEM  = WORK(N+N*N+1)
1464:             NUMRANK = IDNINT(WORK(N+N*N+2))
1465:             DO 6970 p = 1, N
1466:                CALL DCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 )
1467:                CALL DSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 )
1468:  6970       CONTINUE
1469: *
1470:             CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N,
1471:      &           ONE, A, LDA, WORK(N+1), N )
1472:             DO 6972 p = 1, N
1473:                CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV )
1474:  6972       CONTINUE
1475:             TEMP1 = DSQRT(DBLE(N))*EPSLN
1476:             DO 6971 p = 1, N
1477:                XSC = ONE / DNRM2( N, V(1,p), 1 )
1478:                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1479:      &            CALL DSCAL( N, XSC, V(1,p), 1 )
1480:  6971       CONTINUE
1481: *
1482: *           Assemble the left singular vector matrix U (M x N).
1483: *
1484:             IF ( N .LT. M ) THEN
1485:                CALL DLASET( 'A',  M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
1486:                IF ( N .LT. N1 ) THEN
1487:                   CALL DLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
1488:                   CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(NR+1,N+1),LDU )
1489:                END IF
1490:             END IF
1491:             CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1492:      &           LDU, WORK(N+1), LWORK-N, IERR )
1493:             TEMP1 = DSQRT(DBLE(M))*EPSLN
1494:             DO 6973 p = 1, N1
1495:                XSC = ONE / DNRM2( M, U(1,p), 1 )
1496:                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1497:      &            CALL DSCAL( M, XSC, U(1,p), 1 )
1498:  6973       CONTINUE
1499: *
1500:             IF ( ROWPIV )
1501:      &         CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1502: *
1503:          END IF
1504: *
1505: *        end of the  >> almost orthogonal case <<  in the full SVD
1506: *
1507:          ELSE
1508: *
1509: *        This branch deploys a preconditioned Jacobi SVD with explicitly
1510: *        accumulated rotations. It is included as optional, mainly for
1511: *        experimental purposes. It does perfom well, and can also be used.
1512: *        In this implementation, this branch will be automatically activated
1513: *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
1514: *        to be greater than the overflow threshold. This is because the
1515: *        a posteriori computation of the singular vectors assumes robust
1516: *        implementation of BLAS and some LAPACK procedures, capable of working
1517: *        in presence of extreme values. Since that is not always the case, ...
1518: *
1519:          DO 7968 p = 1, NR
1520:             CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 )
1521:  7968    CONTINUE
1522: *
1523:          IF ( L2PERT ) THEN
1524:             XSC = DSQRT(SMALL/EPSLN)
1525:             DO 5969 q = 1, NR
1526:                TEMP1 = XSC*DABS( V(q,q) )
1527:                DO 5968 p = 1, N
1528:                   IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 )
1529:      &                .OR. ( p .LT. q ) )
1530:      &                V(p,q) = DSIGN( TEMP1, V(p,q) )
1531:                   IF ( p. LT. q ) V(p,q) = - V(p,q)
1532:  5968          CONTINUE
1533:  5969       CONTINUE
1534:          ELSE
1535:             CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
1536:          END IF
1537: 
1538:          CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1),
1539:      &        LWORK-2*N, IERR )
1540:          CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N )
1541: *
1542:          DO 7969 p = 1, NR
1543:             CALL DCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 )
1544:  7969    CONTINUE
1545: 
1546:          IF ( L2PERT ) THEN
1547:             XSC = DSQRT(SMALL/EPSLN)
1548:             DO 9970 q = 2, NR
1549:                DO 9971 p = 1, q - 1
1550:                   TEMP1 = XSC * DMIN1(DABS(U(p,p)),DABS(U(q,q)))
1551:                   U(p,q) = - DSIGN( TEMP1, U(q,p) )
1552:  9971          CONTINUE
1553:  9970       CONTINUE
1554:          ELSE
1555:             CALL DLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU )
1556:          END IF
1557: 
1558:          CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA,
1559:      &        N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO )
1560:          SCALEM  = WORK(2*N+N*NR+1)
1561:          NUMRANK = IDNINT(WORK(2*N+N*NR+2))
1562: 
1563:          IF ( NR .LT. N ) THEN
1564:             CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
1565:             CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
1566:             CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
1567:          END IF
1568: 
1569:          CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1),
1570:      &        V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR )
1571: *
1572: *           Permute the rows of V using the (column) permutation from the
1573: *           first QRF. Also, scale the columns to make them unit in
1574: *           Euclidean norm. This applies to all cases.
1575: *
1576:             TEMP1 = DSQRT(DBLE(N)) * EPSLN
1577:             DO 7972 q = 1, N
1578:                DO 8972 p = 1, N
1579:                   WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
1580:  8972          CONTINUE
1581:                DO 8973 p = 1, N
1582:                   V(p,q) = WORK(2*N+N*NR+NR+p)
1583:  8973          CONTINUE
1584:                XSC = ONE / DNRM2( N, V(1,q), 1 )
1585:                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
1586:      &           CALL DSCAL( N, XSC, V(1,q), 1 )
1587:  7972       CONTINUE
1588: *
1589: *           At this moment, V contains the right singular vectors of A.
1590: *           Next, assemble the left singular vector matrix U (M x N).
1591: *
1592:          IF ( N .LT. M ) THEN
1593:             CALL DLASET( 'A',  M-N, N, ZERO, ZERO, U(NR+1,1), LDU )
1594:             IF ( N .LT. N1 ) THEN
1595:                CALL DLASET( 'A',N,  N1-N, ZERO, ZERO,  U(1,N+1),LDU )
1596:                CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(NR+1,N+1),LDU )
1597:             END IF
1598:          END IF
1599: *
1600:          CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U,
1601:      &        LDU, WORK(N+1), LWORK-N, IERR )
1602: *
1603:             IF ( ROWPIV )
1604:      &         CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 )
1605: *
1606: *
1607:          END IF
1608:          IF ( TRANSP ) THEN
1609: *           .. swap U and V because the procedure worked on A^t
1610:             DO 6974 p = 1, N
1611:                CALL DSWAP( N, U(1,p), 1, V(1,p), 1 )
1612:  6974       CONTINUE
1613:          END IF
1614: *
1615:       END IF
1616: *     end of the full SVD
1617: *
1618: *     Undo scaling, if necessary (and possible)
1619: *
1620:       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
1621:          CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
1622:          USCAL1 = ONE
1623:          USCAL2 = ONE
1624:       END IF
1625: *
1626:       IF ( NR .LT. N ) THEN
1627:          DO 3004 p = NR+1, N
1628:             SVA(p) = ZERO
1629:  3004    CONTINUE
1630:       END IF
1631: *
1632:       WORK(1) = USCAL2 * SCALEM
1633:       WORK(2) = USCAL1
1634:       IF ( ERREST ) WORK(3) = SCONDA
1635:       IF ( LSVEC .AND. RSVEC ) THEN
1636:          WORK(4) = CONDR1
1637:          WORK(5) = CONDR2
1638:       END IF
1639:       IF ( L2TRAN ) THEN
1640:          WORK(6) = ENTRA
1641:          WORK(7) = ENTRAT
1642:       END IF
1643: *
1644:       IWORK(1) = NR
1645:       IWORK(2) = NUMRANK
1646:       IWORK(3) = WARNING
1647: *
1648:       RETURN
1649: *     ..
1650: *     .. END OF DGEJSV
1651: *     ..
1652:       END
1653: *
1654: