0001:       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
0002:      +                   LDV, WORK, LWORK, INFO )
0003: *
0004: *  -- LAPACK routine (version 3.2.1)                                    --
0005: *
0006: *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
0007: *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
0008: *  -- April 2009                                                      --
0009: *
0010: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
0011: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
0012: *
0013: * This routine is also part of SIGMA (version 1.23, October 23. 2008.)
0014: * SIGMA is a library of algorithms for highly accurate algorithms for
0015: * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the
0016: * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
0017: *
0018:       IMPLICIT           NONE
0019: *     ..
0020: *     .. Scalar Arguments ..
0021:       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
0022:       CHARACTER*1        JOBA, JOBU, JOBV
0023: *     ..
0024: *     .. Array Arguments ..
0025:       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
0026:      +                   WORK( LWORK )
0027: *     ..
0028: *
0029: *  Purpose
0030: *  =======
0031: *
0032: *  DGESVJ computes the singular value decomposition (SVD) of a real
0033: *  M-by-N matrix A, where M >= N. The SVD of A is written as
0034: *                                     [++]   [xx]   [x0]   [xx]
0035: *               A = U * SIGMA * V^t,  [++] = [xx] * [ox] * [xx]
0036: *                                     [++]   [xx]
0037: *  where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
0038: *  matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
0039: *  of SIGMA are the singular values of A. The columns of U and V are the
0040: *  left and the right singular vectors of A, respectively.
0041: *
0042: *  Further Details
0043: *  ~~~~~~~~~~~~~~~
0044: *  The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
0045: *  rotations. The rotations are implemented as fast scaled rotations of
0046: *  Anda and Park [1]. In the case of underflow of the Jacobi angle, a
0047: *  modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
0048: *  column interchanges of de Rijk [2]. The relative accuracy of the computed
0049: *  singular values and the accuracy of the computed singular vectors (in
0050: *  angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
0051: *  The condition number that determines the accuracy in the full rank case
0052: *  is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
0053: *  spectral condition number. The best performance of this Jacobi SVD
0054: *  procedure is achieved if used in an  accelerated version of Drmac and
0055: *  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
0056: *  Some tunning parameters (marked with [TP]) are available for the
0057: *  implementer.
0058: *  The computational range for the nonzero singular values is the  machine
0059: *  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
0060: *  denormalized singular values can be computed with the corresponding
0061: *  gradual loss of accurate digits.
0062: *
0063: *  Contributors
0064: *  ~~~~~~~~~~~~
0065: *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
0066: *
0067: *  References
0068: *  ~~~~~~~~~~
0069: * [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
0070: *     SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
0071: * [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
0072: *     singular value decomposition on a vector computer.
0073: *     SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
0074: * [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
0075: * [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
0076: *     value computation in floating point arithmetic.
0077: *     SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
0078: * [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
0079: *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
0080: *     LAPACK Working note 169.
0081: * [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
0082: *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
0083: *     LAPACK Working note 170.
0084: * [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
0085: *     QSVD, (H,K)-SVD computations.
0086: *     Department of Mathematics, University of Zagreb, 2008.
0087: *
0088: *  Bugs, Examples and Comments
0089: *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
0090: *  Please report all bugs and send interesting test examples and comments to
0091: *  drmac@math.hr. Thank you.
0092: *
0093: *  Arguments
0094: *  =========
0095: *
0096: *  JOBA    (input) CHARACTER* 1
0097: *          Specifies the structure of A.
0098: *          = 'L': The input matrix A is lower triangular;
0099: *          = 'U': The input matrix A is upper triangular;
0100: *          = 'G': The input matrix A is general M-by-N matrix, M >= N.
0101: *
0102: *  JOBU    (input) CHARACTER*1
0103: *          Specifies whether to compute the left singular vectors
0104: *          (columns of U):
0105: *          = 'U': The left singular vectors corresponding to the nonzero
0106: *                 singular values are computed and returned in the leading
0107: *                 columns of A. See more details in the description of A.
0108: *                 The default numerical orthogonality threshold is set to
0109: *                 approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E').
0110: *          = 'C': Analogous to JOBU='U', except that user can control the
0111: *                 level of numerical orthogonality of the computed left
0112: *                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
0113: *                 CTOL is given on input in the array WORK.
0114: *                 No CTOL smaller than ONE is allowed. CTOL greater
0115: *                 than 1 / EPS is meaningless. The option 'C'
0116: *                 can be used if M*EPS is satisfactory orthogonality
0117: *                 of the computed left singular vectors, so CTOL=M could
0118: *                 save few sweeps of Jacobi rotations.
0119: *                 See the descriptions of A and WORK(1).
0120: *          = 'N': The matrix U is not computed. However, see the
0121: *                 description of A.
0122: *
0123: *  JOBV    (input) CHARACTER*1
0124: *          Specifies whether to compute the right singular vectors, that
0125: *          is, the matrix V:
0126: *          = 'V' : the matrix V is computed and returned in the array V
0127: *          = 'A' : the Jacobi rotations are applied to the MV-by-N
0128: *                  array V. In other words, the right singular vector
0129: *                  matrix V is not computed explicitly, instead it is
0130: *                  applied to an MV-by-N matrix initially stored in the
0131: *                  first MV rows of V.
0132: *          = 'N' : the matrix V is not computed and the array V is not
0133: *                  referenced
0134: *
0135: *  M       (input) INTEGER
0136: *          The number of rows of the input matrix A.  M >= 0.
0137: *
0138: *  N       (input) INTEGER
0139: *          The number of columns of the input matrix A.
0140: *          M >= N >= 0.
0141: *
0142: *  A       (input/output) REAL array, dimension (LDA,N)
0143: *          On entry, the M-by-N matrix A.
0144: *          On exit :
0145: *          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
0146: *                 If INFO .EQ. 0 :
0147: *                 RANKA orthonormal columns of U are returned in the
0148: *                 leading RANKA columns of the array A. Here RANKA <= N
0149: *                 is the number of computed singular values of A that are
0150: *                 above the underflow threshold DLAMCH('S'). The singular
0151: *                 vectors corresponding to underflowed or zero singular
0152: *                 values are not computed. The value of RANKA is returned
0153: *                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
0154: *                 descriptions of SVA and WORK. The computed columns of U
0155: *                 are mutually numerically orthogonal up to approximately
0156: *                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
0157: *                 see the description of JOBU.
0158: *                 If INFO .GT. 0 :
0159: *                 the procedure DGESVJ did not converge in the given number
0160: *                 of iterations (sweeps). In that case, the computed
0161: *                 columns of U may not be orthogonal up to TOL. The output
0162: *                 U (stored in A), SIGMA (given by the computed singular
0163: *                 values in SVA(1:N)) and V is still a decomposition of the
0164: *                 input matrix A in the sense that the residual
0165: *                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
0166: *
0167: *          If JOBU .EQ. 'N' :
0168: *                 If INFO .EQ. 0 :
0169: *                 Note that the left singular vectors are 'for free' in the
0170: *                 one-sided Jacobi SVD algorithm. However, if only the
0171: *                 singular values are needed, the level of numerical
0172: *                 orthogonality of U is not an issue and iterations are
0173: *                 stopped when the columns of the iterated matrix are
0174: *                 numerically orthogonal up to approximately M*EPS. Thus,
0175: *                 on exit, A contains the columns of U scaled with the
0176: *                 corresponding singular values.
0177: *                 If INFO .GT. 0 :
0178: *                 the procedure DGESVJ did not converge in the given number
0179: *                 of iterations (sweeps).
0180: *
0181: *  LDA     (input) INTEGER
0182: *          The leading dimension of the array A.  LDA >= max(1,M).
0183: *
0184: *  SVA     (workspace/output) REAL array, dimension (N)
0185: *          On exit :
0186: *          If INFO .EQ. 0 :
0187: *          depending on the value SCALE = WORK(1), we have:
0188: *                 If SCALE .EQ. ONE :
0189: *                 SVA(1:N) contains the computed singular values of A.
0190: *                 During the computation SVA contains the Euclidean column
0191: *                 norms of the iterated matrices in the array A.
0192: *                 If SCALE .NE. ONE :
0193: *                 The singular values of A are SCALE*SVA(1:N), and this
0194: *                 factored representation is due to the fact that some of the
0195: *                 singular values of A might underflow or overflow.
0196: *          If INFO .GT. 0 :
0197: *          the procedure DGESVJ did not converge in the given number of
0198: *          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
0199: *
0200: *  MV      (input) INTEGER
0201: *          If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ
0202: *          is applied to the first MV rows of V. See the description of JOBV.
0203: *
0204: *  V       (input/output) REAL array, dimension (LDV,N)
0205: *          If JOBV = 'V', then V contains on exit the N-by-N matrix of
0206: *                         the right singular vectors;
0207: *          If JOBV = 'A', then V contains the product of the computed right
0208: *                         singular vector matrix and the initial matrix in
0209: *                         the array V.
0210: *          If JOBV = 'N', then V is not referenced.
0211: *
0212: *  LDV     (input) INTEGER
0213: *          The leading dimension of the array V, LDV .GE. 1.
0214: *          If JOBV .EQ. 'V', then LDV .GE. max(1,N).
0215: *          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
0216: *
0217: *  WORK    (input/workspace/output) REAL array, dimension max(4,M+N).
0218: *          On entry :
0219: *          If JOBU .EQ. 'C' :
0220: *          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
0221: *                    The process stops if all columns of A are mutually
0222: *                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
0223: *                    It is required that CTOL >= ONE, i.e. it is not
0224: *                    allowed to force the routine to obtain orthogonality
0225: *                    below EPSILON.
0226: *          On exit :
0227: *          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
0228: *                    are the computed singular values of A.
0229: *                    (See description of SVA().)
0230: *          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
0231: *                    singular values.
0232: *          WORK(3) = NINT(WORK(3)) is the number of the computed singular
0233: *                    values that are larger than the underflow threshold.
0234: *          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
0235: *                    rotations needed for numerical convergence.
0236: *          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
0237: *                    This is useful information in cases when DGESVJ did
0238: *                    not converge, as it can be used to estimate whether
0239: *                    the output is stil useful and for post festum analysis.
0240: *          WORK(6) = the largest absolute value over all sines of the
0241: *                    Jacobi rotation angles in the last sweep. It can be
0242: *                    useful for a post festum analysis.
0243: *
0244: *  LWORK   length of WORK, WORK >= MAX(6,M+N)
0245: *
0246: *  INFO    (output) INTEGER
0247: *          = 0 : successful exit.
0248: *          < 0 : if INFO = -i, then the i-th argument had an illegal value
0249: *          > 0 : DGESVJ did not converge in the maximal allowed number (30)
0250: *                of sweeps. The output may still be useful. See the
0251: *                description of WORK.
0252: *
0253: *  =====================================================================
0254: *
0255: *     .. Local Parameters ..
0256:       DOUBLE PRECISION   ZERO, HALF, ONE, TWO
0257:       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
0258:      +                   TWO = 2.0D0 )
0259:       INTEGER            NSWEEP
0260:       PARAMETER          ( NSWEEP = 30 )
0261: *     ..
0262: *     .. Local Scalars ..
0263:       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
0264:      +                   BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
0265:      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
0266:      +                   SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
0267:      +                   THSIGN, TOL
0268:       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
0269:      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
0270:      +                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
0271:      +                   SWBAND
0272:       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
0273:      +                   RSVEC, UCTOL, UPPER
0274: *     ..
0275: *     .. Local Arrays ..
0276:       DOUBLE PRECISION   FASTR( 5 )
0277: *     ..
0278: *     .. Intrinsic Functions ..
0279:       INTRINSIC          DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT
0280: *     ..
0281: *     .. External Functions ..
0282: *     ..
0283: *     from BLAS
0284:       DOUBLE PRECISION   DDOT, DNRM2
0285:       EXTERNAL           DDOT, DNRM2
0286:       INTEGER            IDAMAX
0287:       EXTERNAL           IDAMAX
0288: *     from LAPACK
0289:       DOUBLE PRECISION   DLAMCH
0290:       EXTERNAL           DLAMCH
0291:       LOGICAL            LSAME
0292:       EXTERNAL           LSAME
0293: *     ..
0294: *     .. External Subroutines ..
0295: *     ..
0296: *     from BLAS
0297:       EXTERNAL           DAXPY, DCOPY, DROTM, DSCAL, DSWAP
0298: *     from LAPACK
0299:       EXTERNAL           DLASCL, DLASET, DLASSQ, XERBLA
0300: *
0301:       EXTERNAL           DGSVJ0, DGSVJ1
0302: *     ..
0303: *     .. Executable Statements ..
0304: *
0305: *     Test the input arguments
0306: *
0307:       LSVEC = LSAME( JOBU, 'U' )
0308:       UCTOL = LSAME( JOBU, 'C' )
0309:       RSVEC = LSAME( JOBV, 'V' )
0310:       APPLV = LSAME( JOBV, 'A' )
0311:       UPPER = LSAME( JOBA, 'U' )
0312:       LOWER = LSAME( JOBA, 'L' )
0313: *
0314:       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
0315:          INFO = -1
0316:       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
0317:          INFO = -2
0318:       ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
0319:          INFO = -3
0320:       ELSE IF( M.LT.0 ) THEN
0321:          INFO = -4
0322:       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
0323:          INFO = -5
0324:       ELSE IF( LDA.LT.M ) THEN
0325:          INFO = -7
0326:       ELSE IF( MV.LT.0 ) THEN
0327:          INFO = -9
0328:       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
0329:      +         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
0330:          INFO = -11
0331:       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
0332:          INFO = -12
0333:       ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN
0334:          INFO = -13
0335:       ELSE
0336:          INFO = 0
0337:       END IF
0338: *
0339: *     #:(
0340:       IF( INFO.NE.0 ) THEN
0341:          CALL XERBLA( 'DGESVJ', -INFO )
0342:          RETURN
0343:       END IF
0344: *
0345: * #:) Quick return for void matrix
0346: *
0347:       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
0348: *
0349: *     Set numerical parameters
0350: *     The stopping criterion for Jacobi rotations is
0351: *
0352: *     max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
0353: *
0354: *     where EPS is the round-off and CTOL is defined as follows:
0355: *
0356:       IF( UCTOL ) THEN
0357: *        ... user controlled
0358:          CTOL = WORK( 1 )
0359:       ELSE
0360: *        ... default
0361:          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
0362:             CTOL = DSQRT( DBLE( M ) )
0363:          ELSE
0364:             CTOL = DBLE( M )
0365:          END IF
0366:       END IF
0367: *     ... and the machine dependent parameters are
0368: *[!]  (Make sure that DLAMCH() works properly on the target machine.)
0369: *
0370:       EPSILON = DLAMCH( 'Epsilon' )
0371:       ROOTEPS = DSQRT( EPSILON )
0372:       SFMIN = DLAMCH( 'SafeMinimum' )
0373:       ROOTSFMIN = DSQRT( SFMIN )
0374:       SMALL = SFMIN / EPSILON
0375:       BIG = DLAMCH( 'Overflow' )
0376: *     BIG         = ONE    / SFMIN
0377:       ROOTBIG = ONE / ROOTSFMIN
0378:       LARGE = BIG / DSQRT( DBLE( M*N ) )
0379:       BIGTHETA = ONE / ROOTEPS
0380: *
0381:       TOL = CTOL*EPSILON
0382:       ROOTTOL = DSQRT( TOL )
0383: *
0384:       IF( DBLE( M )*EPSILON.GE.ONE ) THEN
0385:          INFO = -5
0386:          CALL XERBLA( 'DGESVJ', -INFO )
0387:          RETURN
0388:       END IF
0389: *
0390: *     Initialize the right singular vector matrix.
0391: *
0392:       IF( RSVEC ) THEN
0393:          MVL = N
0394:          CALL DLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
0395:       ELSE IF( APPLV ) THEN
0396:          MVL = MV
0397:       END IF
0398:       RSVEC = RSVEC .OR. APPLV
0399: *
0400: *     Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
0401: *(!)  If necessary, scale A to protect the largest singular value
0402: *     from overflow. It is possible that saving the largest singular
0403: *     value destroys the information about the small ones.
0404: *     This initial scaling is almost minimal in the sense that the
0405: *     goal is to make sure that no column norm overflows, and that
0406: *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
0407: *     in A are detected, the procedure returns with INFO=-6.
0408: *
0409:       SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) )
0410:       NOSCALE = .TRUE.
0411:       GOSCALE = .TRUE.
0412: *
0413:       IF( LOWER ) THEN
0414: *        the input matrix is M-by-N lower triangular (trapezoidal)
0415:          DO 1874 p = 1, N
0416:             AAPP = ZERO
0417:             AAQQ = ZERO
0418:             CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
0419:             IF( AAPP.GT.BIG ) THEN
0420:                INFO = -6
0421:                CALL XERBLA( 'DGESVJ', -INFO )
0422:                RETURN
0423:             END IF
0424:             AAQQ = DSQRT( AAQQ )
0425:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
0426:                SVA( p ) = AAPP*AAQQ
0427:             ELSE
0428:                NOSCALE = .FALSE.
0429:                SVA( p ) = AAPP*( AAQQ*SCALE )
0430:                IF( GOSCALE ) THEN
0431:                   GOSCALE = .FALSE.
0432:                   DO 1873 q = 1, p - 1
0433:                      SVA( q ) = SVA( q )*SCALE
0434:  1873             CONTINUE
0435:                END IF
0436:             END IF
0437:  1874    CONTINUE
0438:       ELSE IF( UPPER ) THEN
0439: *        the input matrix is M-by-N upper triangular (trapezoidal)
0440:          DO 2874 p = 1, N
0441:             AAPP = ZERO
0442:             AAQQ = ZERO
0443:             CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
0444:             IF( AAPP.GT.BIG ) THEN
0445:                INFO = -6
0446:                CALL XERBLA( 'DGESVJ', -INFO )
0447:                RETURN
0448:             END IF
0449:             AAQQ = DSQRT( AAQQ )
0450:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
0451:                SVA( p ) = AAPP*AAQQ
0452:             ELSE
0453:                NOSCALE = .FALSE.
0454:                SVA( p ) = AAPP*( AAQQ*SCALE )
0455:                IF( GOSCALE ) THEN
0456:                   GOSCALE = .FALSE.
0457:                   DO 2873 q = 1, p - 1
0458:                      SVA( q ) = SVA( q )*SCALE
0459:  2873             CONTINUE
0460:                END IF
0461:             END IF
0462:  2874    CONTINUE
0463:       ELSE
0464: *        the input matrix is M-by-N general dense
0465:          DO 3874 p = 1, N
0466:             AAPP = ZERO
0467:             AAQQ = ZERO
0468:             CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
0469:             IF( AAPP.GT.BIG ) THEN
0470:                INFO = -6
0471:                CALL XERBLA( 'DGESVJ', -INFO )
0472:                RETURN
0473:             END IF
0474:             AAQQ = DSQRT( AAQQ )
0475:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
0476:                SVA( p ) = AAPP*AAQQ
0477:             ELSE
0478:                NOSCALE = .FALSE.
0479:                SVA( p ) = AAPP*( AAQQ*SCALE )
0480:                IF( GOSCALE ) THEN
0481:                   GOSCALE = .FALSE.
0482:                   DO 3873 q = 1, p - 1
0483:                      SVA( q ) = SVA( q )*SCALE
0484:  3873             CONTINUE
0485:                END IF
0486:             END IF
0487:  3874    CONTINUE
0488:       END IF
0489: *
0490:       IF( NOSCALE )SCALE = ONE
0491: *
0492: *     Move the smaller part of the spectrum from the underflow threshold
0493: *(!)  Start by determining the position of the nonzero entries of the
0494: *     array SVA() relative to ( SFMIN, BIG ).
0495: *
0496:       AAPP = ZERO
0497:       AAQQ = BIG
0498:       DO 4781 p = 1, N
0499:          IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) )
0500:          AAPP = DMAX1( AAPP, SVA( p ) )
0501:  4781 CONTINUE
0502: *
0503: * #:) Quick return for zero matrix
0504: *
0505:       IF( AAPP.EQ.ZERO ) THEN
0506:          IF( LSVEC )CALL DLASET( 'G', M, N, ZERO, ONE, A, LDA )
0507:          WORK( 1 ) = ONE
0508:          WORK( 2 ) = ZERO
0509:          WORK( 3 ) = ZERO
0510:          WORK( 4 ) = ZERO
0511:          WORK( 5 ) = ZERO
0512:          WORK( 6 ) = ZERO
0513:          RETURN
0514:       END IF
0515: *
0516: * #:) Quick return for one-column matrix
0517: *
0518:       IF( N.EQ.1 ) THEN
0519:          IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
0520:      +                           A( 1, 1 ), LDA, IERR )
0521:          WORK( 1 ) = ONE / SCALE
0522:          IF( SVA( 1 ).GE.SFMIN ) THEN
0523:             WORK( 2 ) = ONE
0524:          ELSE
0525:             WORK( 2 ) = ZERO
0526:          END IF
0527:          WORK( 3 ) = ZERO
0528:          WORK( 4 ) = ZERO
0529:          WORK( 5 ) = ZERO
0530:          WORK( 6 ) = ZERO
0531:          RETURN
0532:       END IF
0533: *
0534: *     Protect small singular values from underflow, and try to
0535: *     avoid underflows/overflows in computing Jacobi rotations.
0536: *
0537:       SN = DSQRT( SFMIN / EPSILON )
0538:       TEMP1 = DSQRT( BIG / DBLE( N ) )
0539:       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
0540:      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
0541:          TEMP1 = DMIN1( BIG, TEMP1 / AAPP )
0542: *         AAQQ  = AAQQ*TEMP1
0543: *         AAPP  = AAPP*TEMP1
0544:       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
0545:          TEMP1 = DMIN1( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) )
0546: *         AAQQ  = AAQQ*TEMP1
0547: *         AAPP  = AAPP*TEMP1
0548:       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
0549:          TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP )
0550: *         AAQQ  = AAQQ*TEMP1
0551: *         AAPP  = AAPP*TEMP1
0552:       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
0553:          TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) )
0554: *         AAQQ  = AAQQ*TEMP1
0555: *         AAPP  = AAPP*TEMP1
0556:       ELSE
0557:          TEMP1 = ONE
0558:       END IF
0559: *
0560: *     Scale, if necessary
0561: *
0562:       IF( TEMP1.NE.ONE ) THEN
0563:          CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
0564:       END IF
0565:       SCALE = TEMP1*SCALE
0566:       IF( SCALE.NE.ONE ) THEN
0567:          CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
0568:          SCALE = ONE / SCALE
0569:       END IF
0570: *
0571: *     Row-cyclic Jacobi SVD algorithm with column pivoting
0572: *
0573:       EMPTSW = ( N*( N-1 ) ) / 2
0574:       NOTROT = 0
0575:       FASTR( 1 ) = ZERO
0576: *
0577: *     A is represented in factored form A = A * diag(WORK), where diag(WORK)
0578: *     is initialized to identity. WORK is updated during fast scaled
0579: *     rotations.
0580: *
0581:       DO 1868 q = 1, N
0582:          WORK( q ) = ONE
0583:  1868 CONTINUE
0584: *
0585: *
0586:       SWBAND = 3
0587: *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
0588: *     if DGESVJ is used as a computational routine in the preconditioned
0589: *     Jacobi SVD algorithm DGESVJ. For sweeps i=1:SWBAND the procedure
0590: *     works on pivots inside a band-like region around the diagonal.
0591: *     The boundaries are determined dynamically, based on the number of
0592: *     pivots above a threshold.
0593: *
0594:       KBL = MIN0( 8, N )
0595: *[TP] KBL is a tuning parameter that defines the tile size in the
0596: *     tiling of the p-q loops of pivot pairs. In general, an optimal
0597: *     value of KBL depends on the matrix dimensions and on the
0598: *     parameters of the computer's memory.
0599: *
0600:       NBL = N / KBL
0601:       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
0602: *
0603:       BLSKIP = KBL**2
0604: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
0605: *
0606:       ROWSKIP = MIN0( 5, KBL )
0607: *[TP] ROWSKIP is a tuning parameter.
0608: *
0609:       LKAHEAD = 1
0610: *[TP] LKAHEAD is a tuning parameter.
0611: *
0612: *     Quasi block transformations, using the lower (upper) triangular
0613: *     structure of the input matrix. The quasi-block-cycling usually
0614: *     invokes cubic convergence. Big part of this cycle is done inside
0615: *     canonical subspaces of dimensions less than M.
0616: *
0617:       IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
0618: *[TP] The number of partition levels and the actual partition are
0619: *     tuning parameters.
0620:          N4 = N / 4
0621:          N2 = N / 2
0622:          N34 = 3*N4
0623:          IF( APPLV ) THEN
0624:             q = 0
0625:          ELSE
0626:             q = 1
0627:          END IF
0628: *
0629:          IF( LOWER ) THEN
0630: *
0631: *     This works very well on lower triangular matrices, in particular
0632: *     in the framework of the preconditioned Jacobi SVD (xGEJSV).
0633: *     The idea is simple:
0634: *     [+ 0 0 0]   Note that Jacobi transformations of [0 0]
0635: *     [+ + 0 0]                                       [0 0]
0636: *     [+ + x 0]   actually work on [x 0]              [x 0]
0637: *     [+ + x x]                    [x x].             [x x]
0638: *
0639:             CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
0640:      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,
0641:      +                   V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
0642:      +                   2, WORK( N+1 ), LWORK-N, IERR )
0643: *
0644:             CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
0645:      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
0646:      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
0647:      +                   WORK( N+1 ), LWORK-N, IERR )
0648: *
0649:             CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
0650:      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
0651:      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
0652:      +                   WORK( N+1 ), LWORK-N, IERR )
0653: *
0654:             CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
0655:      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,
0656:      +                   V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
0657:      +                   WORK( N+1 ), LWORK-N, IERR )
0658: *
0659:             CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
0660:      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
0661:      +                   IERR )
0662: *
0663:             CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
0664:      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
0665:      +                   LWORK-N, IERR )
0666: *
0667: *
0668:          ELSE IF( UPPER ) THEN
0669: *
0670: *
0671:             CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
0672:      +                   EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
0673:      +                   IERR )
0674: *
0675:             CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
0676:      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
0677:      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
0678:      +                   IERR )
0679: *
0680:             CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
0681:      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
0682:      +                   LWORK-N, IERR )
0683: *
0684:             CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
0685:      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
0686:      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
0687:      +                   WORK( N+1 ), LWORK-N, IERR )
0688: 
0689:          END IF
0690: *
0691:       END IF
0692: *
0693: *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
0694: *
0695:       DO 1993 i = 1, NSWEEP
0696: *
0697: *     .. go go go ...
0698: *
0699:          MXAAPQ = ZERO
0700:          MXSINJ = ZERO
0701:          ISWROT = 0
0702: *
0703:          NOTROT = 0
0704:          PSKIPPED = 0
0705: *
0706: *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
0707: *     1 <= p < q <= N. This is the first step toward a blocked implementation
0708: *     of the rotations. New implementation, based on block transformations,
0709: *     is under development.
0710: *
0711:          DO 2000 ibr = 1, NBL
0712: *
0713:             igl = ( ibr-1 )*KBL + 1
0714: *
0715:             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
0716: *
0717:                igl = igl + ir1*KBL
0718: *
0719:                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
0720: *
0721: *     .. de Rijk's pivoting
0722: *
0723:                   q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
0724:                   IF( p.NE.q ) THEN
0725:                      CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
0726:                      IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
0727:      +                                      V( 1, q ), 1 )
0728:                      TEMP1 = SVA( p )
0729:                      SVA( p ) = SVA( q )
0730:                      SVA( q ) = TEMP1
0731:                      TEMP1 = WORK( p )
0732:                      WORK( p ) = WORK( q )
0733:                      WORK( q ) = TEMP1
0734:                   END IF
0735: *
0736:                   IF( ir1.EQ.0 ) THEN
0737: *
0738: *        Column norms are periodically updated by explicit
0739: *        norm computation.
0740: *        Caveat:
0741: *        Unfortunately, some BLAS implementations compute DNRM2(M,A(1,p),1)
0742: *        as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
0743: *        overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and to
0744: *        underflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
0745: *        Hence, DNRM2 cannot be trusted, not even in the case when
0746: *        the true norm is far from the under(over)flow boundaries.
0747: *        If properly implemented DNRM2 is available, the IF-THEN-ELSE
0748: *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".
0749: *
0750:                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
0751:      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
0752:                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
0753:                      ELSE
0754:                         TEMP1 = ZERO
0755:                         AAPP = ZERO
0756:                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
0757:                         SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
0758:                      END IF
0759:                      AAPP = SVA( p )
0760:                   ELSE
0761:                      AAPP = SVA( p )
0762:                   END IF
0763: *
0764:                   IF( AAPP.GT.ZERO ) THEN
0765: *
0766:                      PSKIPPED = 0
0767: *
0768:                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
0769: *
0770:                         AAQQ = SVA( q )
0771: *
0772:                         IF( AAQQ.GT.ZERO ) THEN
0773: *
0774:                            AAPP0 = AAPP
0775:                            IF( AAQQ.GE.ONE ) THEN
0776:                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
0777:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
0778:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
0779:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
0780:      +                                  AAQQ ) / AAPP
0781:                               ELSE
0782:                                  CALL DCOPY( M, A( 1, p ), 1,
0783:      +                                       WORK( N+1 ), 1 )
0784:                                  CALL DLASCL( 'G', 0, 0, AAPP,
0785:      +                                        WORK( p ), M, 1,
0786:      +                                        WORK( N+1 ), LDA, IERR )
0787:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
0788:      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
0789:                               END IF
0790:                            ELSE
0791:                               ROTOK = AAPP.LE.( AAQQ / SMALL )
0792:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
0793:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
0794:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
0795:      +                                  AAQQ ) / AAPP
0796:                               ELSE
0797:                                  CALL DCOPY( M, A( 1, q ), 1,
0798:      +                                       WORK( N+1 ), 1 )
0799:                                  CALL DLASCL( 'G', 0, 0, AAQQ,
0800:      +                                        WORK( q ), M, 1,
0801:      +                                        WORK( N+1 ), LDA, IERR )
0802:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
0803:      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP
0804:                               END IF
0805:                            END IF
0806: *
0807:                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
0808: *
0809: *        TO rotate or NOT to rotate, THAT is the question ...
0810: *
0811:                            IF( DABS( AAPQ ).GT.TOL ) THEN
0812: *
0813: *           .. rotate
0814: *[RTD]      ROTATED = ROTATED + ONE
0815: *
0816:                               IF( ir1.EQ.0 ) THEN
0817:                                  NOTROT = 0
0818:                                  PSKIPPED = 0
0819:                                  ISWROT = ISWROT + 1
0820:                               END IF
0821: *
0822:                               IF( ROTOK ) THEN
0823: *
0824:                                  AQOAP = AAQQ / AAPP
0825:                                  APOAQ = AAPP / AAQQ
0826:                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /
0827:      +                                   AAPQ
0828: *
0829:                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
0830: *
0831:                                     T = HALF / THETA
0832:                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
0833:                                     FASTR( 4 ) = -T*WORK( q ) /
0834:      +                                           WORK( p )
0835:                                     CALL DROTM( M, A( 1, p ), 1,
0836:      +                                          A( 1, q ), 1, FASTR )
0837:                                     IF( RSVEC )CALL DROTM( MVL,
0838:      +                                              V( 1, p ), 1,
0839:      +                                              V( 1, q ), 1,
0840:      +                                              FASTR )
0841:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
0842:      +                                         ONE+T*APOAQ*AAPQ ) )
0843:                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*
0844:      +                                     AAPQ )
0845:                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
0846: *
0847:                                  ELSE
0848: *
0849: *                 .. choose correct signum for THETA and rotate
0850: *
0851:                                     THSIGN = -DSIGN( ONE, AAPQ )
0852:                                     T = ONE / ( THETA+THSIGN*
0853:      +                                  DSQRT( ONE+THETA*THETA ) )
0854:                                     CS = DSQRT( ONE / ( ONE+T*T ) )
0855:                                     SN = T*CS
0856: *
0857:                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
0858:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
0859:      +                                         ONE+T*APOAQ*AAPQ ) )
0860:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
0861:      +                                     ONE-T*AQOAP*AAPQ ) )
0862: *
0863:                                     APOAQ = WORK( p ) / WORK( q )
0864:                                     AQOAP = WORK( q ) / WORK( p )
0865:                                     IF( WORK( p ).GE.ONE ) THEN
0866:                                        IF( WORK( q ).GE.ONE ) THEN
0867:                                           FASTR( 3 ) = T*APOAQ
0868:                                           FASTR( 4 ) = -T*AQOAP
0869:                                           WORK( p ) = WORK( p )*CS
0870:                                           WORK( q ) = WORK( q )*CS
0871:                                           CALL DROTM( M, A( 1, p ), 1,
0872:      +                                                A( 1, q ), 1,
0873:      +                                                FASTR )
0874:                                           IF( RSVEC )CALL DROTM( MVL,
0875:      +                                        V( 1, p ), 1, V( 1, q ),
0876:      +                                        1, FASTR )
0877:                                        ELSE
0878:                                           CALL DAXPY( M, -T*AQOAP,
0879:      +                                                A( 1, q ), 1,
0880:      +                                                A( 1, p ), 1 )
0881:                                           CALL DAXPY( M, CS*SN*APOAQ,
0882:      +                                                A( 1, p ), 1,
0883:      +                                                A( 1, q ), 1 )
0884:                                           WORK( p ) = WORK( p )*CS
0885:                                           WORK( q ) = WORK( q ) / CS
0886:                                           IF( RSVEC ) THEN
0887:                                              CALL DAXPY( MVL, -T*AQOAP,
0888:      +                                                   V( 1, q ), 1,
0889:      +                                                   V( 1, p ), 1 )
0890:                                              CALL DAXPY( MVL,
0891:      +                                                   CS*SN*APOAQ,
0892:      +                                                   V( 1, p ), 1,
0893:      +                                                   V( 1, q ), 1 )
0894:                                           END IF
0895:                                        END IF
0896:                                     ELSE
0897:                                        IF( WORK( q ).GE.ONE ) THEN
0898:                                           CALL DAXPY( M, T*APOAQ,
0899:      +                                                A( 1, p ), 1,
0900:      +                                                A( 1, q ), 1 )
0901:                                           CALL DAXPY( M, -CS*SN*AQOAP,
0902:      +                                                A( 1, q ), 1,
0903:      +                                                A( 1, p ), 1 )
0904:                                           WORK( p ) = WORK( p ) / CS
0905:                                           WORK( q ) = WORK( q )*CS
0906:                                           IF( RSVEC ) THEN
0907:                                              CALL DAXPY( MVL, T*APOAQ,
0908:      +                                                   V( 1, p ), 1,
0909:      +                                                   V( 1, q ), 1 )
0910:                                              CALL DAXPY( MVL,
0911:      +                                                   -CS*SN*AQOAP,
0912:      +                                                   V( 1, q ), 1,
0913:      +                                                   V( 1, p ), 1 )
0914:                                           END IF
0915:                                        ELSE
0916:                                           IF( WORK( p ).GE.WORK( q ) )
0917:      +                                        THEN
0918:                                              CALL DAXPY( M, -T*AQOAP,
0919:      +                                                   A( 1, q ), 1,
0920:      +                                                   A( 1, p ), 1 )
0921:                                              CALL DAXPY( M, CS*SN*APOAQ,
0922:      +                                                   A( 1, p ), 1,
0923:      +                                                   A( 1, q ), 1 )
0924:                                              WORK( p ) = WORK( p )*CS
0925:                                              WORK( q ) = WORK( q ) / CS
0926:                                              IF( RSVEC ) THEN
0927:                                                 CALL DAXPY( MVL,
0928:      +                                               -T*AQOAP,
0929:      +                                               V( 1, q ), 1,
0930:      +                                               V( 1, p ), 1 )
0931:                                                 CALL DAXPY( MVL,
0932:      +                                               CS*SN*APOAQ,
0933:      +                                               V( 1, p ), 1,
0934:      +                                               V( 1, q ), 1 )
0935:                                              END IF
0936:                                           ELSE
0937:                                              CALL DAXPY( M, T*APOAQ,
0938:      +                                                   A( 1, p ), 1,
0939:      +                                                   A( 1, q ), 1 )
0940:                                              CALL DAXPY( M,
0941:      +                                                   -CS*SN*AQOAP,
0942:      +                                                   A( 1, q ), 1,
0943:      +                                                   A( 1, p ), 1 )
0944:                                              WORK( p ) = WORK( p ) / CS
0945:                                              WORK( q ) = WORK( q )*CS
0946:                                              IF( RSVEC ) THEN
0947:                                                 CALL DAXPY( MVL,
0948:      +                                               T*APOAQ, V( 1, p ),
0949:      +                                               1, V( 1, q ), 1 )
0950:                                                 CALL DAXPY( MVL,
0951:      +                                               -CS*SN*AQOAP,
0952:      +                                               V( 1, q ), 1,
0953:      +                                               V( 1, p ), 1 )
0954:                                              END IF
0955:                                           END IF
0956:                                        END IF
0957:                                     END IF
0958:                                  END IF
0959: *
0960:                               ELSE
0961: *              .. have to use modified Gram-Schmidt like transformation
0962:                                  CALL DCOPY( M, A( 1, p ), 1,
0963:      +                                       WORK( N+1 ), 1 )
0964:                                  CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
0965:      +                                        1, WORK( N+1 ), LDA,
0966:      +                                        IERR )
0967:                                  CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
0968:      +                                        1, A( 1, q ), LDA, IERR )
0969:                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )
0970:                                  CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,
0971:      +                                       A( 1, q ), 1 )
0972:                                  CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
0973:      +                                        1, A( 1, q ), LDA, IERR )
0974:                                  SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
0975:      +                                      ONE-AAPQ*AAPQ ) )
0976:                                  MXSINJ = DMAX1( MXSINJ, SFMIN )
0977:                               END IF
0978: *           END IF ROTOK THEN ... ELSE
0979: *
0980: *           In the case of cancellation in updating SVA(q), SVA(p)
0981: *           recompute SVA(q), SVA(p).
0982: *
0983:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
0984:      +                            THEN
0985:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
0986:      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
0987:                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
0988:      +                                         WORK( q )
0989:                                  ELSE
0990:                                     T = ZERO
0991:                                     AAQQ = ZERO
0992:                                     CALL DLASSQ( M, A( 1, q ), 1, T,
0993:      +                                           AAQQ )
0994:                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
0995:                                  END IF
0996:                               END IF
0997:                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
0998:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
0999:      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1000:                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
1001:      +                                     WORK( p )
1002:                                  ELSE
1003:                                     T = ZERO
1004:                                     AAPP = ZERO
1005:                                     CALL DLASSQ( M, A( 1, p ), 1, T,
1006:      +                                           AAPP )
1007:                                     AAPP = T*DSQRT( AAPP )*WORK( p )
1008:                                  END IF
1009:                                  SVA( p ) = AAPP
1010:                               END IF
1011: *
1012:                            ELSE
1013: *        A(:,p) and A(:,q) already numerically orthogonal
1014:                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1015: *[RTD]      SKIPPED  = SKIPPED  + 1
1016:                               PSKIPPED = PSKIPPED + 1
1017:                            END IF
1018:                         ELSE
1019: *        A(:,q) is zero column
1020:                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1021:                            PSKIPPED = PSKIPPED + 1
1022:                         END IF
1023: *
1024:                         IF( ( i.LE.SWBAND ) .AND.
1025:      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1026:                            IF( ir1.EQ.0 )AAPP = -AAPP
1027:                            NOTROT = 0
1028:                            GO TO 2103
1029:                         END IF
1030: *
1031:  2002                CONTINUE
1032: *     END q-LOOP
1033: *
1034:  2103                CONTINUE
1035: *     bailed out of q-loop
1036: *
1037:                      SVA( p ) = AAPP
1038: *
1039:                   ELSE
1040:                      SVA( p ) = AAPP
1041:                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
1042:      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
1043:                   END IF
1044: *
1045:  2001          CONTINUE
1046: *     end of the p-loop
1047: *     end of doing the block ( ibr, ibr )
1048:  1002       CONTINUE
1049: *     end of ir1-loop
1050: *
1051: * ... go to the off diagonal blocks
1052: *
1053:             igl = ( ibr-1 )*KBL + 1
1054: *
1055:             DO 2010 jbc = ibr + 1, NBL
1056: *
1057:                jgl = ( jbc-1 )*KBL + 1
1058: *
1059: *        doing the block at ( ibr, jbc )
1060: *
1061:                IJBLSK = 0
1062:                DO 2100 p = igl, MIN0( igl+KBL-1, N )
1063: *
1064:                   AAPP = SVA( p )
1065:                   IF( AAPP.GT.ZERO ) THEN
1066: *
1067:                      PSKIPPED = 0
1068: *
1069:                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
1070: *
1071:                         AAQQ = SVA( q )
1072:                         IF( AAQQ.GT.ZERO ) THEN
1073:                            AAPP0 = AAPP
1074: *
1075: *     .. M x 2 Jacobi SVD ..
1076: *
1077: *        Safe Gram matrix computation
1078: *
1079:                            IF( AAQQ.GE.ONE ) THEN
1080:                               IF( AAPP.GE.AAQQ ) THEN
1081:                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
1082:                               ELSE
1083:                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
1084:                               END IF
1085:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
1086:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
1087:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
1088:      +                                  AAQQ ) / AAPP
1089:                               ELSE
1090:                                  CALL DCOPY( M, A( 1, p ), 1,
1091:      +                                       WORK( N+1 ), 1 )
1092:                                  CALL DLASCL( 'G', 0, 0, AAPP,
1093:      +                                        WORK( p ), M, 1,
1094:      +                                        WORK( N+1 ), LDA, IERR )
1095:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
1096:      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
1097:                               END IF
1098:                            ELSE
1099:                               IF( AAPP.GE.AAQQ ) THEN
1100:                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
1101:                               ELSE
1102:                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
1103:                               END IF
1104:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
1105:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
1106:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
1107:      +                                  AAQQ ) / AAPP
1108:                               ELSE
1109:                                  CALL DCOPY( M, A( 1, q ), 1,
1110:      +                                       WORK( N+1 ), 1 )
1111:                                  CALL DLASCL( 'G', 0, 0, AAQQ,
1112:      +                                        WORK( q ), M, 1,
1113:      +                                        WORK( N+1 ), LDA, IERR )
1114:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
1115:      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP
1116:                               END IF
1117:                            END IF
1118: *
1119:                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
1120: *
1121: *        TO rotate or NOT to rotate, THAT is the question ...
1122: *
1123:                            IF( DABS( AAPQ ).GT.TOL ) THEN
1124:                               NOTROT = 0
1125: *[RTD]      ROTATED  = ROTATED + 1
1126:                               PSKIPPED = 0
1127:                               ISWROT = ISWROT + 1
1128: *
1129:                               IF( ROTOK ) THEN
1130: *
1131:                                  AQOAP = AAQQ / AAPP
1132:                                  APOAQ = AAPP / AAQQ
1133:                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /
1134:      +                                   AAPQ
1135:                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
1136: *
1137:                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
1138:                                     T = HALF / THETA
1139:                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
1140:                                     FASTR( 4 ) = -T*WORK( q ) /
1141:      +                                           WORK( p )
1142:                                     CALL DROTM( M, A( 1, p ), 1,
1143:      +                                          A( 1, q ), 1, FASTR )
1144:                                     IF( RSVEC )CALL DROTM( MVL,
1145:      +                                              V( 1, p ), 1,
1146:      +                                              V( 1, q ), 1,
1147:      +                                              FASTR )
1148:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1149:      +                                         ONE+T*APOAQ*AAPQ ) )
1150:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
1151:      +                                     ONE-T*AQOAP*AAPQ ) )
1152:                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
1153:                                  ELSE
1154: *
1155: *                 .. choose correct signum for THETA and rotate
1156: *
1157:                                     THSIGN = -DSIGN( ONE, AAPQ )
1158:                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
1159:                                     T = ONE / ( THETA+THSIGN*
1160:      +                                  DSQRT( ONE+THETA*THETA ) )
1161:                                     CS = DSQRT( ONE / ( ONE+T*T ) )
1162:                                     SN = T*CS
1163:                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
1164:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1165:      +                                         ONE+T*APOAQ*AAPQ ) )
1166:                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*
1167:      +                                     AAPQ )
1168: *
1169:                                     APOAQ = WORK( p ) / WORK( q )
1170:                                     AQOAP = WORK( q ) / WORK( p )
1171:                                     IF( WORK( p ).GE.ONE ) THEN
1172: *
1173:                                        IF( WORK( q ).GE.ONE ) THEN
1174:                                           FASTR( 3 ) = T*APOAQ
1175:                                           FASTR( 4 ) = -T*AQOAP
1176:                                           WORK( p ) = WORK( p )*CS
1177:                                           WORK( q ) = WORK( q )*CS
1178:                                           CALL DROTM( M, A( 1, p ), 1,
1179:      +                                                A( 1, q ), 1,
1180:      +                                                FASTR )
1181:                                           IF( RSVEC )CALL DROTM( MVL,
1182:      +                                        V( 1, p ), 1, V( 1, q ),
1183:      +                                        1, FASTR )
1184:                                        ELSE
1185:                                           CALL DAXPY( M, -T*AQOAP,
1186:      +                                                A( 1, q ), 1,
1187:      +                                                A( 1, p ), 1 )
1188:                                           CALL DAXPY( M, CS*SN*APOAQ,
1189:      +                                                A( 1, p ), 1,
1190:      +                                                A( 1, q ), 1 )
1191:                                           IF( RSVEC ) THEN
1192:                                              CALL DAXPY( MVL, -T*AQOAP,
1193:      +                                                   V( 1, q ), 1,
1194:      +                                                   V( 1, p ), 1 )
1195:                                              CALL DAXPY( MVL,
1196:      +                                                   CS*SN*APOAQ,
1197:      +                                                   V( 1, p ), 1,
1198:      +                                                   V( 1, q ), 1 )
1199:                                           END IF
1200:                                           WORK( p ) = WORK( p )*CS
1201:                                           WORK( q ) = WORK( q ) / CS
1202:                                        END IF
1203:                                     ELSE
1204:                                        IF( WORK( q ).GE.ONE ) THEN
1205:                                           CALL DAXPY( M, T*APOAQ,
1206:      +                                                A( 1, p ), 1,
1207:      +                                                A( 1, q ), 1 )
1208:                                           CALL DAXPY( M, -CS*SN*AQOAP,
1209:      +                                                A( 1, q ), 1,
1210:      +                                                A( 1, p ), 1 )
1211:                                           IF( RSVEC ) THEN
1212:                                              CALL DAXPY( MVL, T*APOAQ,
1213:      +                                                   V( 1, p ), 1,
1214:      +                                                   V( 1, q ), 1 )
1215:                                              CALL DAXPY( MVL,
1216:      +                                                   -CS*SN*AQOAP,
1217:      +                                                   V( 1, q ), 1,
1218:      +                                                   V( 1, p ), 1 )
1219:                                           END IF
1220:                                           WORK( p ) = WORK( p ) / CS
1221:                                           WORK( q ) = WORK( q )*CS
1222:                                        ELSE
1223:                                           IF( WORK( p ).GE.WORK( q ) )
1224:      +                                        THEN
1225:                                              CALL DAXPY( M, -T*AQOAP,
1226:      +                                                   A( 1, q ), 1,
1227:      +                                                   A( 1, p ), 1 )
1228:                                              CALL DAXPY( M, CS*SN*APOAQ,
1229:      +                                                   A( 1, p ), 1,
1230:      +                                                   A( 1, q ), 1 )
1231:                                              WORK( p ) = WORK( p )*CS
1232:                                              WORK( q ) = WORK( q ) / CS
1233:                                              IF( RSVEC ) THEN
1234:                                                 CALL DAXPY( MVL,
1235:      +                                               -T*AQOAP,
1236:      +                                               V( 1, q ), 1,
1237:      +                                               V( 1, p ), 1 )
1238:                                                 CALL DAXPY( MVL,
1239:      +                                               CS*SN*APOAQ,
1240:      +                                               V( 1, p ), 1,
1241:      +                                               V( 1, q ), 1 )
1242:                                              END IF
1243:                                           ELSE
1244:                                              CALL DAXPY( M, T*APOAQ,
1245:      +                                                   A( 1, p ), 1,
1246:      +                                                   A( 1, q ), 1 )
1247:                                              CALL DAXPY( M,
1248:      +                                                   -CS*SN*AQOAP,
1249:      +                                                   A( 1, q ), 1,
1250:      +                                                   A( 1, p ), 1 )
1251:                                              WORK( p ) = WORK( p ) / CS
1252:                                              WORK( q ) = WORK( q )*CS
1253:                                              IF( RSVEC ) THEN
1254:                                                 CALL DAXPY( MVL,
1255:      +                                               T*APOAQ, V( 1, p ),
1256:      +                                               1, V( 1, q ), 1 )
1257:                                                 CALL DAXPY( MVL,
1258:      +                                               -CS*SN*AQOAP,
1259:      +                                               V( 1, q ), 1,
1260:      +                                               V( 1, p ), 1 )
1261:                                              END IF
1262:                                           END IF
1263:                                        END IF
1264:                                     END IF
1265:                                  END IF
1266: *
1267:                               ELSE
1268:                                  IF( AAPP.GT.AAQQ ) THEN
1269:                                     CALL DCOPY( M, A( 1, p ), 1,
1270:      +                                          WORK( N+1 ), 1 )
1271:                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
1272:      +                                           M, 1, WORK( N+1 ), LDA,
1273:      +                                           IERR )
1274:                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
1275:      +                                           M, 1, A( 1, q ), LDA,
1276:      +                                           IERR )
1277:                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )
1278:                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),
1279:      +                                          1, A( 1, q ), 1 )
1280:                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
1281:      +                                           M, 1, A( 1, q ), LDA,
1282:      +                                           IERR )
1283:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1284:      +                                         ONE-AAPQ*AAPQ ) )
1285:                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
1286:                                  ELSE
1287:                                     CALL DCOPY( M, A( 1, q ), 1,
1288:      +                                          WORK( N+1 ), 1 )
1289:                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
1290:      +                                           M, 1, WORK( N+1 ), LDA,
1291:      +                                           IERR )
1292:                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
1293:      +                                           M, 1, A( 1, p ), LDA,
1294:      +                                           IERR )
1295:                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )
1296:                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),
1297:      +                                          1, A( 1, p ), 1 )
1298:                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,
1299:      +                                           M, 1, A( 1, p ), LDA,
1300:      +                                           IERR )
1301:                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
1302:      +                                         ONE-AAPQ*AAPQ ) )
1303:                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
1304:                                  END IF
1305:                               END IF
1306: *           END IF ROTOK THEN ... ELSE
1307: *
1308: *           In the case of cancellation in updating SVA(q)
1309: *           .. recompute SVA(q)
1310:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1311:      +                            THEN
1312:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
1313:      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
1314:                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
1315:      +                                         WORK( q )
1316:                                  ELSE
1317:                                     T = ZERO
1318:                                     AAQQ = ZERO
1319:                                     CALL DLASSQ( M, A( 1, q ), 1, T,
1320:      +                                           AAQQ )
1321:                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
1322:                                  END IF
1323:                               END IF
1324:                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
1325:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
1326:      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1327:                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
1328:      +                                     WORK( p )
1329:                                  ELSE
1330:                                     T = ZERO
1331:                                     AAPP = ZERO
1332:                                     CALL DLASSQ( M, A( 1, p ), 1, T,
1333:      +                                           AAPP )
1334:                                     AAPP = T*DSQRT( AAPP )*WORK( p )
1335:                                  END IF
1336:                                  SVA( p ) = AAPP
1337:                               END IF
1338: *              end of OK rotation
1339:                            ELSE
1340:                               NOTROT = NOTROT + 1
1341: *[RTD]      SKIPPED  = SKIPPED  + 1
1342:                               PSKIPPED = PSKIPPED + 1
1343:                               IJBLSK = IJBLSK + 1
1344:                            END IF
1345:                         ELSE
1346:                            NOTROT = NOTROT + 1
1347:                            PSKIPPED = PSKIPPED + 1
1348:                            IJBLSK = IJBLSK + 1
1349:                         END IF
1350: *
1351:                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
1352:      +                      THEN
1353:                            SVA( p ) = AAPP
1354:                            NOTROT = 0
1355:                            GO TO 2011
1356:                         END IF
1357:                         IF( ( i.LE.SWBAND ) .AND.
1358:      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1359:                            AAPP = -AAPP
1360:                            NOTROT = 0
1361:                            GO TO 2203
1362:                         END IF
1363: *
1364:  2200                CONTINUE
1365: *        end of the q-loop
1366:  2203                CONTINUE
1367: *
1368:                      SVA( p ) = AAPP
1369: *
1370:                   ELSE
1371: *
1372:                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1373:      +                   MIN0( jgl+KBL-1, N ) - jgl + 1
1374:                      IF( AAPP.LT.ZERO )NOTROT = 0
1375: *
1376:                   END IF
1377: *
1378:  2100          CONTINUE
1379: *     end of the p-loop
1380:  2010       CONTINUE
1381: *     end of the jbc-loop
1382:  2011       CONTINUE
1383: *2011 bailed out of the jbc-loop
1384:             DO 2012 p = igl, MIN0( igl+KBL-1, N )
1385:                SVA( p ) = DABS( SVA( p ) )
1386:  2012       CONTINUE
1387: ***
1388:  2000    CONTINUE
1389: *2000 :: end of the ibr-loop
1390: *
1391: *     .. update SVA(N)
1392:          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
1393:      +       THEN
1394:             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
1395:          ELSE
1396:             T = ZERO
1397:             AAPP = ZERO
1398:             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
1399:             SVA( N ) = T*DSQRT( AAPP )*WORK( N )
1400:          END IF
1401: *
1402: *     Additional steering devices
1403: *
1404:          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
1405:      +       ( ISWROT.LE.N ) ) )SWBAND = i
1406: *
1407:          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
1408:      +       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1409:             GO TO 1994
1410:          END IF
1411: *
1412:          IF( NOTROT.GE.EMPTSW )GO TO 1994
1413: *
1414:  1993 CONTINUE
1415: *     end i=1:NSWEEP loop
1416: *
1417: * #:( Reaching this point means that the procedure has not converged.
1418:       INFO = NSWEEP - 1
1419:       GO TO 1995
1420: *
1421:  1994 CONTINUE
1422: * #:) Reaching this point means numerical convergence after the i-th
1423: *     sweep.
1424: *
1425:       INFO = 0
1426: * #:) INFO = 0 confirms successful iterations.
1427:  1995 CONTINUE
1428: *
1429: *     Sort the singular values and find how many are above
1430: *     the underflow threshold.
1431: *
1432:       N2 = 0
1433:       N4 = 0
1434:       DO 5991 p = 1, N - 1
1435:          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
1436:          IF( p.NE.q ) THEN
1437:             TEMP1 = SVA( p )
1438:             SVA( p ) = SVA( q )
1439:             SVA( q ) = TEMP1
1440:             TEMP1 = WORK( p )
1441:             WORK( p ) = WORK( q )
1442:             WORK( q ) = TEMP1
1443:             CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
1444:             IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
1445:          END IF
1446:          IF( SVA( p ).NE.ZERO ) THEN
1447:             N4 = N4 + 1
1448:             IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
1449:          END IF
1450:  5991 CONTINUE
1451:       IF( SVA( N ).NE.ZERO ) THEN
1452:          N4 = N4 + 1
1453:          IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
1454:       END IF
1455: *
1456: *     Normalize the left singular vectors.
1457: *
1458:       IF( LSVEC .OR. UCTOL ) THEN
1459:          DO 1998 p = 1, N2
1460:             CALL DSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
1461:  1998    CONTINUE
1462:       END IF
1463: *
1464: *     Scale the product of Jacobi rotations (assemble the fast rotations).
1465: *
1466:       IF( RSVEC ) THEN
1467:          IF( APPLV ) THEN
1468:             DO 2398 p = 1, N
1469:                CALL DSCAL( MVL, WORK( p ), V( 1, p ), 1 )
1470:  2398       CONTINUE
1471:          ELSE
1472:             DO 2399 p = 1, N
1473:                TEMP1 = ONE / DNRM2( MVL, V( 1, p ), 1 )
1474:                CALL DSCAL( MVL, TEMP1, V( 1, p ), 1 )
1475:  2399       CONTINUE
1476:          END IF
1477:       END IF
1478: *
1479: *     Undo scaling, if necessary (and possible).
1480:       IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
1481:      +    SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
1482:      +    ( SFMIN / SCALE ) ) ) ) THEN
1483:          DO 2400 p = 1, N
1484:             SVA( p ) = SCALE*SVA( p )
1485:  2400    CONTINUE
1486:          SCALE = ONE
1487:       END IF
1488: *
1489:       WORK( 1 ) = SCALE
1490: *     The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
1491: *     then some of the singular values may overflow or underflow and
1492: *     the spectrum is given in this factored representation.
1493: *
1494:       WORK( 2 ) = DBLE( N4 )
1495: *     N4 is the number of computed nonzero singular values of A.
1496: *
1497:       WORK( 3 ) = DBLE( N2 )
1498: *     N2 is the number of singular values of A greater than SFMIN.
1499: *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
1500: *     that may carry some information.
1501: *
1502:       WORK( 4 ) = DBLE( i )
1503: *     i is the index of the last sweep before declaring convergence.
1504: *
1505:       WORK( 5 ) = MXAAPQ
1506: *     MXAAPQ is the largest absolute value of scaled pivots in the
1507: *     last sweep
1508: *
1509:       WORK( 6 ) = MXSINJ
1510: *     MXSINJ is the largest absolute value of the sines of Jacobi angles
1511: *     in the last sweep
1512: *
1513:       RETURN
1514: *     ..
1515: *     .. END OF DGESVJ
1516: *     ..
1517:       END
1518: