0001:       SUBROUTINE SGESVJ( 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:       REAL               A( LDA, * ), SVA( N ), V( LDV, * ),
0026:      +                   WORK( LWORK )
0027: *     ..
0028: *
0029: *  Purpose
0030: *  =======
0031: *
0032: *  SGESVJ 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=SQRT(M), EPS=SLAMCH('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 SLAMCH('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=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
0157: *                 see the description of JOBU.
0158: *                 If INFO .GT. 0,
0159: *                 the procedure SGESVJ 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: *          If JOBU .EQ. 'N':
0167: *                 If INFO .EQ. 0 :
0168: *                 Note that the left singular vectors are 'for free' in the
0169: *                 one-sided Jacobi SVD algorithm. However, if only the
0170: *                 singular values are needed, the level of numerical
0171: *                 orthogonality of U is not an issue and iterations are
0172: *                 stopped when the columns of the iterated matrix are
0173: *                 numerically orthogonal up to approximately M*EPS. Thus,
0174: *                 on exit, A contains the columns of U scaled with the
0175: *                 corresponding singular values.
0176: *                 If INFO .GT. 0 :
0177: *                 the procedure SGESVJ did not converge in the given number
0178: *                 of iterations (sweeps).
0179: *
0180: *  LDA     (input) INTEGER
0181: *          The leading dimension of the array A.  LDA >= max(1,M).
0182: *
0183: *  SVA     (workspace/output) REAL array, dimension (N)
0184: *          On exit,
0185: *          If INFO .EQ. 0 :
0186: *          depending on the value SCALE = WORK(1), we have:
0187: *                 If SCALE .EQ. ONE:
0188: *                 SVA(1:N) contains the computed singular values of A.
0189: *                 During the computation SVA contains the Euclidean column
0190: *                 norms of the iterated matrices in the array A.
0191: *                 If SCALE .NE. ONE:
0192: *                 The singular values of A are SCALE*SVA(1:N), and this
0193: *                 factored representation is due to the fact that some of the
0194: *                 singular values of A might underflow or overflow.
0195: *
0196: *          If INFO .GT. 0 :
0197: *          the procedure SGESVJ 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 SGESVJ
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=SLAMCH('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 vcalues 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 SGESVJ 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 : SGESVJ 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: *     .. Local Parameters ..
0255:       REAL               ZERO, HALF, ONE, TWO
0256:       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
0257:      +                   TWO = 2.0E0 )
0258:       INTEGER            NSWEEP
0259:       PARAMETER          ( NSWEEP = 30 )
0260: *     ..
0261: *     .. Local Scalars ..
0262:       REAL               AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
0263:      +                   BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,
0264:      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
0265:      +                   SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,
0266:      +                   THSIGN, TOL
0267:       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
0268:      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
0269:      +                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
0270:      +                   SWBAND
0271:       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
0272:      +                   RSVEC, UCTOL, UPPER
0273: *     ..
0274: *     .. Local Arrays ..
0275:       REAL               FASTR( 5 )
0276: *     ..
0277: *     .. Intrinsic Functions ..
0278:       INTRINSIC          ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT
0279: *     ..
0280: *     .. External Functions ..
0281: *     from BLAS
0282:       REAL               SDOT, SNRM2
0283:       EXTERNAL           SDOT, SNRM2
0284:       INTEGER            ISAMAX
0285:       EXTERNAL           ISAMAX
0286: *     from LAPACK
0287:       REAL               SLAMCH
0288:       EXTERNAL           SLAMCH
0289:       LOGICAL            LSAME
0290:       EXTERNAL           LSAME
0291: *     ..
0292: *     .. External Subroutines ..
0293: *     from BLAS
0294:       EXTERNAL           SAXPY, SCOPY, SROTM, SSCAL, SSWAP
0295: *     from LAPACK
0296:       EXTERNAL           SLASCL, SLASET, SLASSQ, XERBLA
0297: *
0298:       EXTERNAL           SGSVJ0, SGSVJ1
0299: *     ..
0300: *     .. Executable Statements ..
0301: *
0302: *     Test the input arguments
0303: *
0304:       LSVEC = LSAME( JOBU, 'U' )
0305:       UCTOL = LSAME( JOBU, 'C' )
0306:       RSVEC = LSAME( JOBV, 'V' )
0307:       APPLV = LSAME( JOBV, 'A' )
0308:       UPPER = LSAME( JOBA, 'U' )
0309:       LOWER = LSAME( JOBA, 'L' )
0310: *
0311:       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
0312:          INFO = -1
0313:       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
0314:          INFO = -2
0315:       ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
0316:          INFO = -3
0317:       ELSE IF( M.LT.0 ) THEN
0318:          INFO = -4
0319:       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
0320:          INFO = -5
0321:       ELSE IF( LDA.LT.M ) THEN
0322:          INFO = -7
0323:       ELSE IF( MV.LT.0 ) THEN
0324:          INFO = -9
0325:       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
0326:      +         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
0327:          INFO = -11
0328:       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
0329:          INFO = -12
0330:       ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN
0331:          INFO = -13
0332:       ELSE
0333:          INFO = 0
0334:       END IF
0335: *
0336: *     #:(
0337:       IF( INFO.NE.0 ) THEN
0338:          CALL XERBLA( 'SGESVJ', -INFO )
0339:          RETURN
0340:       END IF
0341: *
0342: * #:) Quick return for void matrix
0343: *
0344:       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
0345: *
0346: *     Set numerical parameters
0347: *     The stopping criterion for Jacobi rotations is
0348: *
0349: *     max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
0350: *
0351: *     where EPS is the round-off and CTOL is defined as follows:
0352: *
0353:       IF( UCTOL ) THEN
0354: *        ... user controlled
0355:          CTOL = WORK( 1 )
0356:       ELSE
0357: *        ... default
0358:          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
0359:             CTOL = SQRT( FLOAT( M ) )
0360:          ELSE
0361:             CTOL = FLOAT( M )
0362:          END IF
0363:       END IF
0364: *     ... and the machine dependent parameters are
0365: *[!]  (Make sure that SLAMCH() works properly on the target machine.)
0366: *
0367:       EPSILON = SLAMCH( 'Epsilon' )
0368:       ROOTEPS = SQRT( EPSILON )
0369:       SFMIN = SLAMCH( 'SafeMinimum' )
0370:       ROOTSFMIN = SQRT( SFMIN )
0371:       SMALL = SFMIN / EPSILON
0372:       BIG = SLAMCH( 'Overflow' )
0373:       ROOTBIG = ONE / ROOTSFMIN
0374:       LARGE = BIG / SQRT( FLOAT( M*N ) )
0375:       BIGTHETA = ONE / ROOTEPS
0376: *
0377:       TOL = CTOL*EPSILON
0378:       ROOTTOL = SQRT( TOL )
0379: *
0380:       IF( FLOAT( M )*EPSILON.GE.ONE ) THEN
0381:          INFO = -5
0382:          CALL XERBLA( 'SGESVJ', -INFO )
0383:          RETURN
0384:       END IF
0385: *
0386: *     Initialize the right singular vector matrix.
0387: *
0388:       IF( RSVEC ) THEN
0389:          MVL = N
0390:          CALL SLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
0391:       ELSE IF( APPLV ) THEN
0392:          MVL = MV
0393:       END IF
0394:       RSVEC = RSVEC .OR. APPLV
0395: *
0396: *     Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
0397: *(!)  If necessary, scale A to protect the largest singular value
0398: *     from overflow. It is possible that saving the largest singular
0399: *     value destroys the information about the small ones.
0400: *     This initial scaling is almost minimal in the sense that the
0401: *     goal is to make sure that no column norm overflows, and that
0402: *     SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
0403: *     in A are detected, the procedure returns with INFO=-6.
0404: *
0405:       SCALE = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
0406:       NOSCALE = .TRUE.
0407:       GOSCALE = .TRUE.
0408: *
0409:       IF( LOWER ) THEN
0410: *        the input matrix is M-by-N lower triangular (trapezoidal)
0411:          DO 1874 p = 1, N
0412:             AAPP = ZERO
0413:             AAQQ = ZERO
0414:             CALL SLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
0415:             IF( AAPP.GT.BIG ) THEN
0416:                INFO = -6
0417:                CALL XERBLA( 'SGESVJ', -INFO )
0418:                RETURN
0419:             END IF
0420:             AAQQ = SQRT( AAQQ )
0421:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
0422:                SVA( p ) = AAPP*AAQQ
0423:             ELSE
0424:                NOSCALE = .FALSE.
0425:                SVA( p ) = AAPP*( AAQQ*SCALE )
0426:                IF( GOSCALE ) THEN
0427:                   GOSCALE = .FALSE.
0428:                   DO 1873 q = 1, p - 1
0429:                      SVA( q ) = SVA( q )*SCALE
0430:  1873             CONTINUE
0431:                END IF
0432:             END IF
0433:  1874    CONTINUE
0434:       ELSE IF( UPPER ) THEN
0435: *        the input matrix is M-by-N upper triangular (trapezoidal)
0436:          DO 2874 p = 1, N
0437:             AAPP = ZERO
0438:             AAQQ = ZERO
0439:             CALL SLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
0440:             IF( AAPP.GT.BIG ) THEN
0441:                INFO = -6
0442:                CALL XERBLA( 'SGESVJ', -INFO )
0443:                RETURN
0444:             END IF
0445:             AAQQ = SQRT( AAQQ )
0446:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
0447:                SVA( p ) = AAPP*AAQQ
0448:             ELSE
0449:                NOSCALE = .FALSE.
0450:                SVA( p ) = AAPP*( AAQQ*SCALE )
0451:                IF( GOSCALE ) THEN
0452:                   GOSCALE = .FALSE.
0453:                   DO 2873 q = 1, p - 1
0454:                      SVA( q ) = SVA( q )*SCALE
0455:  2873             CONTINUE
0456:                END IF
0457:             END IF
0458:  2874    CONTINUE
0459:       ELSE
0460: *        the input matrix is M-by-N general dense
0461:          DO 3874 p = 1, N
0462:             AAPP = ZERO
0463:             AAQQ = ZERO
0464:             CALL SLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
0465:             IF( AAPP.GT.BIG ) THEN
0466:                INFO = -6
0467:                CALL XERBLA( 'SGESVJ', -INFO )
0468:                RETURN
0469:             END IF
0470:             AAQQ = SQRT( AAQQ )
0471:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
0472:                SVA( p ) = AAPP*AAQQ
0473:             ELSE
0474:                NOSCALE = .FALSE.
0475:                SVA( p ) = AAPP*( AAQQ*SCALE )
0476:                IF( GOSCALE ) THEN
0477:                   GOSCALE = .FALSE.
0478:                   DO 3873 q = 1, p - 1
0479:                      SVA( q ) = SVA( q )*SCALE
0480:  3873             CONTINUE
0481:                END IF
0482:             END IF
0483:  3874    CONTINUE
0484:       END IF
0485: *
0486:       IF( NOSCALE )SCALE = ONE
0487: *
0488: *     Move the smaller part of the spectrum from the underflow threshold
0489: *(!)  Start by determining the position of the nonzero entries of the
0490: *     array SVA() relative to ( SFMIN, BIG ).
0491: *
0492:       AAPP = ZERO
0493:       AAQQ = BIG
0494:       DO 4781 p = 1, N
0495:          IF( SVA( p ).NE.ZERO )AAQQ = AMIN1( AAQQ, SVA( p ) )
0496:          AAPP = AMAX1( AAPP, SVA( p ) )
0497:  4781 CONTINUE
0498: *
0499: * #:) Quick return for zero matrix
0500: *
0501:       IF( AAPP.EQ.ZERO ) THEN
0502:          IF( LSVEC )CALL SLASET( 'G', M, N, ZERO, ONE, A, LDA )
0503:          WORK( 1 ) = ONE
0504:          WORK( 2 ) = ZERO
0505:          WORK( 3 ) = ZERO
0506:          WORK( 4 ) = ZERO
0507:          WORK( 5 ) = ZERO
0508:          WORK( 6 ) = ZERO
0509:          RETURN
0510:       END IF
0511: *
0512: * #:) Quick return for one-column matrix
0513: *
0514:       IF( N.EQ.1 ) THEN
0515:          IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,
0516:      +                           A( 1, 1 ), LDA, IERR )
0517:          WORK( 1 ) = ONE / SCALE
0518:          IF( SVA( 1 ).GE.SFMIN ) THEN
0519:             WORK( 2 ) = ONE
0520:          ELSE
0521:             WORK( 2 ) = ZERO
0522:          END IF
0523:          WORK( 3 ) = ZERO
0524:          WORK( 4 ) = ZERO
0525:          WORK( 5 ) = ZERO
0526:          WORK( 6 ) = ZERO
0527:          RETURN
0528:       END IF
0529: *
0530: *     Protect small singular values from underflow, and try to
0531: *     avoid underflows/overflows in computing Jacobi rotations.
0532: *
0533:       SN = SQRT( SFMIN / EPSILON )
0534:       TEMP1 = SQRT( BIG / FLOAT( N ) )
0535:       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
0536:      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
0537:          TEMP1 = AMIN1( BIG, TEMP1 / AAPP )
0538: *         AAQQ  = AAQQ*TEMP1
0539: *         AAPP  = AAPP*TEMP1
0540:       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
0541:          TEMP1 = AMIN1( SN / AAQQ, BIG / ( AAPP*SQRT( FLOAT( N ) ) ) )
0542: *         AAQQ  = AAQQ*TEMP1
0543: *         AAPP  = AAPP*TEMP1
0544:       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
0545:          TEMP1 = AMAX1( SN / AAQQ, TEMP1 / AAPP )
0546: *         AAQQ  = AAQQ*TEMP1
0547: *         AAPP  = AAPP*TEMP1
0548:       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
0549:          TEMP1 = AMIN1( SN / AAQQ, BIG / ( SQRT( FLOAT( N ) )*AAPP ) )
0550: *         AAQQ  = AAQQ*TEMP1
0551: *         AAPP  = AAPP*TEMP1
0552:       ELSE
0553:          TEMP1 = ONE
0554:       END IF
0555: *
0556: *     Scale, if necessary
0557: *
0558:       IF( TEMP1.NE.ONE ) THEN
0559:          CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
0560:       END IF
0561:       SCALE = TEMP1*SCALE
0562:       IF( SCALE.NE.ONE ) THEN
0563:          CALL SLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )
0564:          SCALE = ONE / SCALE
0565:       END IF
0566: *
0567: *     Row-cyclic Jacobi SVD algorithm with column pivoting
0568: *
0569:       EMPTSW = ( N*( N-1 ) ) / 2
0570:       NOTROT = 0
0571:       FASTR( 1 ) = ZERO
0572: *
0573: *     A is represented in factored form A = A * diag(WORK), where diag(WORK)
0574: *     is initialized to identity. WORK is updated during fast scaled
0575: *     rotations.
0576: *
0577:       DO 1868 q = 1, N
0578:          WORK( q ) = ONE
0579:  1868 CONTINUE
0580: *
0581: *
0582:       SWBAND = 3
0583: *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
0584: *     if SGESVJ is used as a computational routine in the preconditioned
0585: *     Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
0586: *     works on pivots inside a band-like region around the diagonal.
0587: *     The boundaries are determined dynamically, based on the number of
0588: *     pivots above a threshold.
0589: *
0590:       KBL = MIN0( 8, N )
0591: *[TP] KBL is a tuning parameter that defines the tile size in the
0592: *     tiling of the p-q loops of pivot pairs. In general, an optimal
0593: *     value of KBL depends on the matrix dimensions and on the
0594: *     parameters of the computer's memory.
0595: *
0596:       NBL = N / KBL
0597:       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
0598: *
0599:       BLSKIP = KBL**2
0600: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
0601: *
0602:       ROWSKIP = MIN0( 5, KBL )
0603: *[TP] ROWSKIP is a tuning parameter.
0604: *
0605:       LKAHEAD = 1
0606: *[TP] LKAHEAD is a tuning parameter.
0607: *
0608: *     Quasi block transformations, using the lower (upper) triangular
0609: *     structure of the input matrix. The quasi-block-cycling usually
0610: *     invokes cubic convergence. Big part of this cycle is done inside
0611: *     canonical subspaces of dimensions less than M.
0612: *
0613:       IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
0614: *[TP] The number of partition levels and the actual partition are
0615: *     tuning parameters.
0616:          N4 = N / 4
0617:          N2 = N / 2
0618:          N34 = 3*N4
0619:          IF( APPLV ) THEN
0620:             q = 0
0621:          ELSE
0622:             q = 1
0623:          END IF
0624: *
0625:          IF( LOWER ) THEN
0626: *
0627: *     This works very well on lower triangular matrices, in particular
0628: *     in the framework of the preconditioned Jacobi SVD (xGEJSV).
0629: *     The idea is simple:
0630: *     [+ 0 0 0]   Note that Jacobi transformations of [0 0]
0631: *     [+ + 0 0]                                       [0 0]
0632: *     [+ + x 0]   actually work on [x 0]              [x 0]
0633: *     [+ + x x]                    [x x].             [x x]
0634: *
0635:             CALL SGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
0636:      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,
0637:      +                   V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,
0638:      +                   2, WORK( N+1 ), LWORK-N, IERR )
0639: *
0640:             CALL SGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
0641:      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
0642:      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,
0643:      +                   WORK( N+1 ), LWORK-N, IERR )
0644: *
0645:             CALL SGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
0646:      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
0647:      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
0648:      +                   WORK( N+1 ), LWORK-N, IERR )
0649: *
0650:             CALL SGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
0651:      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,
0652:      +                   V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,
0653:      +                   WORK( N+1 ), LWORK-N, IERR )
0654: *
0655:             CALL SGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
0656:      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
0657:      +                   IERR )
0658: *
0659:             CALL SGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
0660:      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
0661:      +                   LWORK-N, IERR )
0662: *
0663: *
0664:          ELSE IF( UPPER ) THEN
0665: *
0666: *
0667:             CALL SGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
0668:      +                   EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
0669:      +                   IERR )
0670: *
0671:             CALL SGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
0672:      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
0673:      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
0674:      +                   IERR )
0675: *
0676:             CALL SGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
0677:      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),
0678:      +                   LWORK-N, IERR )
0679: *
0680:             CALL SGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
0681:      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
0682:      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,
0683:      +                   WORK( N+1 ), LWORK-N, IERR )
0684: 
0685:          END IF
0686: *
0687:       END IF
0688: *
0689: *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
0690: *
0691:       DO 1993 i = 1, NSWEEP
0692: *     .. go go go ...
0693: *
0694:          MXAAPQ = ZERO
0695:          MXSINJ = ZERO
0696:          ISWROT = 0
0697: *
0698:          NOTROT = 0
0699:          PSKIPPED = 0
0700: *
0701: *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
0702: *     1 <= p < q <= N. This is the first step toward a blocked implementation
0703: *     of the rotations. New implementation, based on block transformations,
0704: *     is under development.
0705: *
0706:          DO 2000 ibr = 1, NBL
0707: *
0708:             igl = ( ibr-1 )*KBL + 1
0709: *
0710:             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
0711: *
0712:                igl = igl + ir1*KBL
0713: *
0714:                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
0715: *
0716: *     .. de Rijk's pivoting
0717: *
0718:                   q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
0719:                   IF( p.NE.q ) THEN
0720:                      CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
0721:                      IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1,
0722:      +                                      V( 1, q ), 1 )
0723:                      TEMP1 = SVA( p )
0724:                      SVA( p ) = SVA( q )
0725:                      SVA( q ) = TEMP1
0726:                      TEMP1 = WORK( p )
0727:                      WORK( p ) = WORK( q )
0728:                      WORK( q ) = TEMP1
0729:                   END IF
0730: *
0731:                   IF( ir1.EQ.0 ) THEN
0732: *
0733: *        Column norms are periodically updated by explicit
0734: *        norm computation.
0735: *        Caveat:
0736: *        Unfortunately, some BLAS implementations compute SNRM2(M,A(1,p),1)
0737: *        as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
0738: *        overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
0739: *        underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
0740: *        Hence, SNRM2 cannot be trusted, not even in the case when
0741: *        the true norm is far from the under(over)flow boundaries.
0742: *        If properly implemented SNRM2 is available, the IF-THEN-ELSE
0743: *        below should read "AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)".
0744: *
0745:                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
0746:      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
0747:                         SVA( p ) = SNRM2( M, A( 1, p ), 1 )*WORK( p )
0748:                      ELSE
0749:                         TEMP1 = ZERO
0750:                         AAPP = ZERO
0751:                         CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
0752:                         SVA( p ) = TEMP1*SQRT( AAPP )*WORK( p )
0753:                      END IF
0754:                      AAPP = SVA( p )
0755:                   ELSE
0756:                      AAPP = SVA( p )
0757:                   END IF
0758: *
0759:                   IF( AAPP.GT.ZERO ) THEN
0760: *
0761:                      PSKIPPED = 0
0762: *
0763:                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
0764: *
0765:                         AAQQ = SVA( q )
0766: *
0767:                         IF( AAQQ.GT.ZERO ) THEN
0768: *
0769:                            AAPP0 = AAPP
0770:                            IF( AAQQ.GE.ONE ) THEN
0771:                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
0772:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
0773:                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
0774:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
0775:      +                                  AAQQ ) / AAPP
0776:                               ELSE
0777:                                  CALL SCOPY( M, A( 1, p ), 1,
0778:      +                                       WORK( N+1 ), 1 )
0779:                                  CALL SLASCL( 'G', 0, 0, AAPP,
0780:      +                                        WORK( p ), M, 1,
0781:      +                                        WORK( N+1 ), LDA, IERR )
0782:                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
0783:      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
0784:                               END IF
0785:                            ELSE
0786:                               ROTOK = AAPP.LE.( AAQQ / SMALL )
0787:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
0788:                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
0789:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
0790:      +                                  AAQQ ) / AAPP
0791:                               ELSE
0792:                                  CALL SCOPY( M, A( 1, q ), 1,
0793:      +                                       WORK( N+1 ), 1 )
0794:                                  CALL SLASCL( 'G', 0, 0, AAQQ,
0795:      +                                        WORK( q ), M, 1,
0796:      +                                        WORK( N+1 ), LDA, IERR )
0797:                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
0798:      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP
0799:                               END IF
0800:                            END IF
0801: *
0802:                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
0803: *
0804: *        TO rotate or NOT to rotate, THAT is the question ...
0805: *
0806:                            IF( ABS( AAPQ ).GT.TOL ) THEN
0807: *
0808: *           .. rotate
0809: *[RTD]      ROTATED = ROTATED + ONE
0810: *
0811:                               IF( ir1.EQ.0 ) THEN
0812:                                  NOTROT = 0
0813:                                  PSKIPPED = 0
0814:                                  ISWROT = ISWROT + 1
0815:                               END IF
0816: *
0817:                               IF( ROTOK ) THEN
0818: *
0819:                                  AQOAP = AAQQ / AAPP
0820:                                  APOAQ = AAPP / AAQQ
0821:                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
0822: *
0823:                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
0824: *
0825:                                     T = HALF / THETA
0826:                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
0827:                                     FASTR( 4 ) = -T*WORK( q ) /
0828:      +                                           WORK( p )
0829:                                     CALL SROTM( M, A( 1, p ), 1,
0830:      +                                          A( 1, q ), 1, FASTR )
0831:                                     IF( RSVEC )CALL SROTM( MVL,
0832:      +                                              V( 1, p ), 1,
0833:      +                                              V( 1, q ), 1,
0834:      +                                              FASTR )
0835:                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
0836:      +                                         ONE+T*APOAQ*AAPQ ) )
0837:                                     AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
0838:                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
0839: *
0840:                                  ELSE
0841: *
0842: *                 .. choose correct signum for THETA and rotate
0843: *
0844:                                     THSIGN = -SIGN( ONE, AAPQ )
0845:                                     T = ONE / ( THETA+THSIGN*
0846:      +                                  SQRT( ONE+THETA*THETA ) )
0847:                                     CS = SQRT( ONE / ( ONE+T*T ) )
0848:                                     SN = T*CS
0849: *
0850:                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
0851:                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
0852:      +                                         ONE+T*APOAQ*AAPQ ) )
0853:                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
0854:      +                                     ONE-T*AQOAP*AAPQ ) )
0855: *
0856:                                     APOAQ = WORK( p ) / WORK( q )
0857:                                     AQOAP = WORK( q ) / WORK( p )
0858:                                     IF( WORK( p ).GE.ONE ) THEN
0859:                                        IF( WORK( q ).GE.ONE ) THEN
0860:                                           FASTR( 3 ) = T*APOAQ
0861:                                           FASTR( 4 ) = -T*AQOAP
0862:                                           WORK( p ) = WORK( p )*CS
0863:                                           WORK( q ) = WORK( q )*CS
0864:                                           CALL SROTM( M, A( 1, p ), 1,
0865:      +                                                A( 1, q ), 1,
0866:      +                                                FASTR )
0867:                                           IF( RSVEC )CALL SROTM( MVL,
0868:      +                                        V( 1, p ), 1, V( 1, q ),
0869:      +                                        1, FASTR )
0870:                                        ELSE
0871:                                           CALL SAXPY( M, -T*AQOAP,
0872:      +                                                A( 1, q ), 1,
0873:      +                                                A( 1, p ), 1 )
0874:                                           CALL SAXPY( M, CS*SN*APOAQ,
0875:      +                                                A( 1, p ), 1,
0876:      +                                                A( 1, q ), 1 )
0877:                                           WORK( p ) = WORK( p )*CS
0878:                                           WORK( q ) = WORK( q ) / CS
0879:                                           IF( RSVEC ) THEN
0880:                                              CALL SAXPY( MVL, -T*AQOAP,
0881:      +                                                   V( 1, q ), 1,
0882:      +                                                   V( 1, p ), 1 )
0883:                                              CALL SAXPY( MVL,
0884:      +                                                   CS*SN*APOAQ,
0885:      +                                                   V( 1, p ), 1,
0886:      +                                                   V( 1, q ), 1 )
0887:                                           END IF
0888:                                        END IF
0889:                                     ELSE
0890:                                        IF( WORK( q ).GE.ONE ) THEN
0891:                                           CALL SAXPY( M, T*APOAQ,
0892:      +                                                A( 1, p ), 1,
0893:      +                                                A( 1, q ), 1 )
0894:                                           CALL SAXPY( M, -CS*SN*AQOAP,
0895:      +                                                A( 1, q ), 1,
0896:      +                                                A( 1, p ), 1 )
0897:                                           WORK( p ) = WORK( p ) / CS
0898:                                           WORK( q ) = WORK( q )*CS
0899:                                           IF( RSVEC ) THEN
0900:                                              CALL SAXPY( MVL, T*APOAQ,
0901:      +                                                   V( 1, p ), 1,
0902:      +                                                   V( 1, q ), 1 )
0903:                                              CALL SAXPY( MVL,
0904:      +                                                   -CS*SN*AQOAP,
0905:      +                                                   V( 1, q ), 1,
0906:      +                                                   V( 1, p ), 1 )
0907:                                           END IF
0908:                                        ELSE
0909:                                           IF( WORK( p ).GE.WORK( q ) )
0910:      +                                        THEN
0911:                                              CALL SAXPY( M, -T*AQOAP,
0912:      +                                                   A( 1, q ), 1,
0913:      +                                                   A( 1, p ), 1 )
0914:                                              CALL SAXPY( M, CS*SN*APOAQ,
0915:      +                                                   A( 1, p ), 1,
0916:      +                                                   A( 1, q ), 1 )
0917:                                              WORK( p ) = WORK( p )*CS
0918:                                              WORK( q ) = WORK( q ) / CS
0919:                                              IF( RSVEC ) THEN
0920:                                                 CALL SAXPY( MVL,
0921:      +                                               -T*AQOAP,
0922:      +                                               V( 1, q ), 1,
0923:      +                                               V( 1, p ), 1 )
0924:                                                 CALL SAXPY( MVL,
0925:      +                                               CS*SN*APOAQ,
0926:      +                                               V( 1, p ), 1,
0927:      +                                               V( 1, q ), 1 )
0928:                                              END IF
0929:                                           ELSE
0930:                                              CALL SAXPY( M, T*APOAQ,
0931:      +                                                   A( 1, p ), 1,
0932:      +                                                   A( 1, q ), 1 )
0933:                                              CALL SAXPY( M,
0934:      +                                                   -CS*SN*AQOAP,
0935:      +                                                   A( 1, q ), 1,
0936:      +                                                   A( 1, p ), 1 )
0937:                                              WORK( p ) = WORK( p ) / CS
0938:                                              WORK( q ) = WORK( q )*CS
0939:                                              IF( RSVEC ) THEN
0940:                                                 CALL SAXPY( MVL,
0941:      +                                               T*APOAQ, V( 1, p ),
0942:      +                                               1, V( 1, q ), 1 )
0943:                                                 CALL SAXPY( MVL,
0944:      +                                               -CS*SN*AQOAP,
0945:      +                                               V( 1, q ), 1,
0946:      +                                               V( 1, p ), 1 )
0947:                                              END IF
0948:                                           END IF
0949:                                        END IF
0950:                                     END IF
0951:                                  END IF
0952: *
0953:                               ELSE
0954: *              .. have to use modified Gram-Schmidt like transformation
0955:                                  CALL SCOPY( M, A( 1, p ), 1,
0956:      +                                       WORK( N+1 ), 1 )
0957:                                  CALL SLASCL( 'G', 0, 0, AAPP, ONE, M,
0958:      +                                        1, WORK( N+1 ), LDA,
0959:      +                                        IERR )
0960:                                  CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M,
0961:      +                                        1, A( 1, q ), LDA, IERR )
0962:                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )
0963:                                  CALL SAXPY( M, TEMP1, WORK( N+1 ), 1,
0964:      +                                       A( 1, q ), 1 )
0965:                                  CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M,
0966:      +                                        1, A( 1, q ), LDA, IERR )
0967:                                  SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
0968:      +                                      ONE-AAPQ*AAPQ ) )
0969:                                  MXSINJ = AMAX1( MXSINJ, SFMIN )
0970:                               END IF
0971: *           END IF ROTOK THEN ... ELSE
0972: *
0973: *           In the case of cancellation in updating SVA(q), SVA(p)
0974: *           recompute SVA(q), SVA(p).
0975: *
0976:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
0977:      +                            THEN
0978:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
0979:      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
0980:                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
0981:      +                                         WORK( q )
0982:                                  ELSE
0983:                                     T = ZERO
0984:                                     AAQQ = ZERO
0985:                                     CALL SLASSQ( M, A( 1, q ), 1, T,
0986:      +                                           AAQQ )
0987:                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
0988:                                  END IF
0989:                               END IF
0990:                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
0991:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
0992:      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
0993:                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
0994:      +                                     WORK( p )
0995:                                  ELSE
0996:                                     T = ZERO
0997:                                     AAPP = ZERO
0998:                                     CALL SLASSQ( M, A( 1, p ), 1, T,
0999:      +                                           AAPP )
1000:                                     AAPP = T*SQRT( AAPP )*WORK( p )
1001:                                  END IF
1002:                                  SVA( p ) = AAPP
1003:                               END IF
1004: *
1005:                            ELSE
1006: *        A(:,p) and A(:,q) already numerically orthogonal
1007:                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1008: *[RTD]      SKIPPED  = SKIPPED  + 1
1009:                               PSKIPPED = PSKIPPED + 1
1010:                            END IF
1011:                         ELSE
1012: *        A(:,q) is zero column
1013:                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1014:                            PSKIPPED = PSKIPPED + 1
1015:                         END IF
1016: *
1017:                         IF( ( i.LE.SWBAND ) .AND.
1018:      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1019:                            IF( ir1.EQ.0 )AAPP = -AAPP
1020:                            NOTROT = 0
1021:                            GO TO 2103
1022:                         END IF
1023: *
1024:  2002                CONTINUE
1025: *     END q-LOOP
1026: *
1027:  2103                CONTINUE
1028: *     bailed out of q-loop
1029: *
1030:                      SVA( p ) = AAPP
1031: *
1032:                   ELSE
1033:                      SVA( p ) = AAPP
1034:                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
1035:      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
1036:                   END IF
1037: *
1038:  2001          CONTINUE
1039: *     end of the p-loop
1040: *     end of doing the block ( ibr, ibr )
1041:  1002       CONTINUE
1042: *     end of ir1-loop
1043: *
1044: * ... go to the off diagonal blocks
1045: *
1046:             igl = ( ibr-1 )*KBL + 1
1047: *
1048:             DO 2010 jbc = ibr + 1, NBL
1049: *
1050:                jgl = ( jbc-1 )*KBL + 1
1051: *
1052: *        doing the block at ( ibr, jbc )
1053: *
1054:                IJBLSK = 0
1055:                DO 2100 p = igl, MIN0( igl+KBL-1, N )
1056: *
1057:                   AAPP = SVA( p )
1058:                   IF( AAPP.GT.ZERO ) THEN
1059: *
1060:                      PSKIPPED = 0
1061: *
1062:                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
1063: *
1064:                         AAQQ = SVA( q )
1065:                         IF( AAQQ.GT.ZERO ) THEN
1066:                            AAPP0 = AAPP
1067: *
1068: *     .. M x 2 Jacobi SVD ..
1069: *
1070: *        Safe Gram matrix computation
1071: *
1072:                            IF( AAQQ.GE.ONE ) THEN
1073:                               IF( AAPP.GE.AAQQ ) THEN
1074:                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
1075:                               ELSE
1076:                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
1077:                               END IF
1078:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
1079:                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
1080:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
1081:      +                                  AAQQ ) / AAPP
1082:                               ELSE
1083:                                  CALL SCOPY( M, A( 1, p ), 1,
1084:      +                                       WORK( N+1 ), 1 )
1085:                                  CALL SLASCL( 'G', 0, 0, AAPP,
1086:      +                                        WORK( p ), M, 1,
1087:      +                                        WORK( N+1 ), LDA, IERR )
1088:                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
1089:      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
1090:                               END IF
1091:                            ELSE
1092:                               IF( AAPP.GE.AAQQ ) THEN
1093:                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
1094:                               ELSE
1095:                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
1096:                               END IF
1097:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
1098:                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
1099:      +                                  q ), 1 )*WORK( p )*WORK( q ) /
1100:      +                                  AAQQ ) / AAPP
1101:                               ELSE
1102:                                  CALL SCOPY( M, A( 1, q ), 1,
1103:      +                                       WORK( N+1 ), 1 )
1104:                                  CALL SLASCL( 'G', 0, 0, AAQQ,
1105:      +                                        WORK( q ), M, 1,
1106:      +                                        WORK( N+1 ), LDA, IERR )
1107:                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
1108:      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP
1109:                               END IF
1110:                            END IF
1111: *
1112:                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
1113: *
1114: *        TO rotate or NOT to rotate, THAT is the question ...
1115: *
1116:                            IF( ABS( AAPQ ).GT.TOL ) THEN
1117:                               NOTROT = 0
1118: *[RTD]      ROTATED  = ROTATED + 1
1119:                               PSKIPPED = 0
1120:                               ISWROT = ISWROT + 1
1121: *
1122:                               IF( ROTOK ) THEN
1123: *
1124:                                  AQOAP = AAQQ / AAPP
1125:                                  APOAQ = AAPP / AAQQ
1126:                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
1127:                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
1128: *
1129:                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
1130:                                     T = HALF / THETA
1131:                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
1132:                                     FASTR( 4 ) = -T*WORK( q ) /
1133:      +                                           WORK( p )
1134:                                     CALL SROTM( M, A( 1, p ), 1,
1135:      +                                          A( 1, q ), 1, FASTR )
1136:                                     IF( RSVEC )CALL SROTM( MVL,
1137:      +                                              V( 1, p ), 1,
1138:      +                                              V( 1, q ), 1,
1139:      +                                              FASTR )
1140:                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
1141:      +                                         ONE+T*APOAQ*AAPQ ) )
1142:                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
1143:      +                                     ONE-T*AQOAP*AAPQ ) )
1144:                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
1145:                                  ELSE
1146: *
1147: *                 .. choose correct signum for THETA and rotate
1148: *
1149:                                     THSIGN = -SIGN( ONE, AAPQ )
1150:                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
1151:                                     T = ONE / ( THETA+THSIGN*
1152:      +                                  SQRT( ONE+THETA*THETA ) )
1153:                                     CS = SQRT( ONE / ( ONE+T*T ) )
1154:                                     SN = T*CS
1155:                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
1156:                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
1157:      +                                         ONE+T*APOAQ*AAPQ ) )
1158:                                     AAPP = AAPP*SQRT( ONE-T*AQOAP*AAPQ )
1159: *
1160:                                     APOAQ = WORK( p ) / WORK( q )
1161:                                     AQOAP = WORK( q ) / WORK( p )
1162:                                     IF( WORK( p ).GE.ONE ) THEN
1163: *
1164:                                        IF( WORK( q ).GE.ONE ) THEN
1165:                                           FASTR( 3 ) = T*APOAQ
1166:                                           FASTR( 4 ) = -T*AQOAP
1167:                                           WORK( p ) = WORK( p )*CS
1168:                                           WORK( q ) = WORK( q )*CS
1169:                                           CALL SROTM( M, A( 1, p ), 1,
1170:      +                                                A( 1, q ), 1,
1171:      +                                                FASTR )
1172:                                           IF( RSVEC )CALL SROTM( MVL,
1173:      +                                        V( 1, p ), 1, V( 1, q ),
1174:      +                                        1, FASTR )
1175:                                        ELSE
1176:                                           CALL SAXPY( M, -T*AQOAP,
1177:      +                                                A( 1, q ), 1,
1178:      +                                                A( 1, p ), 1 )
1179:                                           CALL SAXPY( M, CS*SN*APOAQ,
1180:      +                                                A( 1, p ), 1,
1181:      +                                                A( 1, q ), 1 )
1182:                                           IF( RSVEC ) THEN
1183:                                              CALL SAXPY( MVL, -T*AQOAP,
1184:      +                                                   V( 1, q ), 1,
1185:      +                                                   V( 1, p ), 1 )
1186:                                              CALL SAXPY( MVL,
1187:      +                                                   CS*SN*APOAQ,
1188:      +                                                   V( 1, p ), 1,
1189:      +                                                   V( 1, q ), 1 )
1190:                                           END IF
1191:                                           WORK( p ) = WORK( p )*CS
1192:                                           WORK( q ) = WORK( q ) / CS
1193:                                        END IF
1194:                                     ELSE
1195:                                        IF( WORK( q ).GE.ONE ) THEN
1196:                                           CALL SAXPY( M, T*APOAQ,
1197:      +                                                A( 1, p ), 1,
1198:      +                                                A( 1, q ), 1 )
1199:                                           CALL SAXPY( M, -CS*SN*AQOAP,
1200:      +                                                A( 1, q ), 1,
1201:      +                                                A( 1, p ), 1 )
1202:                                           IF( RSVEC ) THEN
1203:                                              CALL SAXPY( MVL, T*APOAQ,
1204:      +                                                   V( 1, p ), 1,
1205:      +                                                   V( 1, q ), 1 )
1206:                                              CALL SAXPY( MVL,
1207:      +                                                   -CS*SN*AQOAP,
1208:      +                                                   V( 1, q ), 1,
1209:      +                                                   V( 1, p ), 1 )
1210:                                           END IF
1211:                                           WORK( p ) = WORK( p ) / CS
1212:                                           WORK( q ) = WORK( q )*CS
1213:                                        ELSE
1214:                                           IF( WORK( p ).GE.WORK( q ) )
1215:      +                                        THEN
1216:                                              CALL SAXPY( M, -T*AQOAP,
1217:      +                                                   A( 1, q ), 1,
1218:      +                                                   A( 1, p ), 1 )
1219:                                              CALL SAXPY( M, CS*SN*APOAQ,
1220:      +                                                   A( 1, p ), 1,
1221:      +                                                   A( 1, q ), 1 )
1222:                                              WORK( p ) = WORK( p )*CS
1223:                                              WORK( q ) = WORK( q ) / CS
1224:                                              IF( RSVEC ) THEN
1225:                                                 CALL SAXPY( MVL,
1226:      +                                               -T*AQOAP,
1227:      +                                               V( 1, q ), 1,
1228:      +                                               V( 1, p ), 1 )
1229:                                                 CALL SAXPY( MVL,
1230:      +                                               CS*SN*APOAQ,
1231:      +                                               V( 1, p ), 1,
1232:      +                                               V( 1, q ), 1 )
1233:                                              END IF
1234:                                           ELSE
1235:                                              CALL SAXPY( M, T*APOAQ,
1236:      +                                                   A( 1, p ), 1,
1237:      +                                                   A( 1, q ), 1 )
1238:                                              CALL SAXPY( M,
1239:      +                                                   -CS*SN*AQOAP,
1240:      +                                                   A( 1, q ), 1,
1241:      +                                                   A( 1, p ), 1 )
1242:                                              WORK( p ) = WORK( p ) / CS
1243:                                              WORK( q ) = WORK( q )*CS
1244:                                              IF( RSVEC ) THEN
1245:                                                 CALL SAXPY( MVL,
1246:      +                                               T*APOAQ, V( 1, p ),
1247:      +                                               1, V( 1, q ), 1 )
1248:                                                 CALL SAXPY( MVL,
1249:      +                                               -CS*SN*AQOAP,
1250:      +                                               V( 1, q ), 1,
1251:      +                                               V( 1, p ), 1 )
1252:                                              END IF
1253:                                           END IF
1254:                                        END IF
1255:                                     END IF
1256:                                  END IF
1257: *
1258:                               ELSE
1259:                                  IF( AAPP.GT.AAQQ ) THEN
1260:                                     CALL SCOPY( M, A( 1, p ), 1,
1261:      +                                          WORK( N+1 ), 1 )
1262:                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
1263:      +                                           M, 1, WORK( N+1 ), LDA,
1264:      +                                           IERR )
1265:                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
1266:      +                                           M, 1, A( 1, q ), LDA,
1267:      +                                           IERR )
1268:                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )
1269:                                     CALL SAXPY( M, TEMP1, WORK( N+1 ),
1270:      +                                          1, A( 1, q ), 1 )
1271:                                     CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
1272:      +                                           M, 1, A( 1, q ), LDA,
1273:      +                                           IERR )
1274:                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
1275:      +                                         ONE-AAPQ*AAPQ ) )
1276:                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
1277:                                  ELSE
1278:                                     CALL SCOPY( M, A( 1, q ), 1,
1279:      +                                          WORK( N+1 ), 1 )
1280:                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
1281:      +                                           M, 1, WORK( N+1 ), LDA,
1282:      +                                           IERR )
1283:                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
1284:      +                                           M, 1, A( 1, p ), LDA,
1285:      +                                           IERR )
1286:                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )
1287:                                     CALL SAXPY( M, TEMP1, WORK( N+1 ),
1288:      +                                          1, A( 1, p ), 1 )
1289:                                     CALL SLASCL( 'G', 0, 0, ONE, AAPP,
1290:      +                                           M, 1, A( 1, p ), LDA,
1291:      +                                           IERR )
1292:                                     SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
1293:      +                                         ONE-AAPQ*AAPQ ) )
1294:                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
1295:                                  END IF
1296:                               END IF
1297: *           END IF ROTOK THEN ... ELSE
1298: *
1299: *           In the case of cancellation in updating SVA(q)
1300: *           .. recompute SVA(q)
1301:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1302:      +                            THEN
1303:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
1304:      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
1305:                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
1306:      +                                         WORK( q )
1307:                                  ELSE
1308:                                     T = ZERO
1309:                                     AAQQ = ZERO
1310:                                     CALL SLASSQ( M, A( 1, q ), 1, T,
1311:      +                                           AAQQ )
1312:                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
1313:                                  END IF
1314:                               END IF
1315:                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
1316:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
1317:      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1318:                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
1319:      +                                     WORK( p )
1320:                                  ELSE
1321:                                     T = ZERO
1322:                                     AAPP = ZERO
1323:                                     CALL SLASSQ( M, A( 1, p ), 1, T,
1324:      +                                           AAPP )
1325:                                     AAPP = T*SQRT( AAPP )*WORK( p )
1326:                                  END IF
1327:                                  SVA( p ) = AAPP
1328:                               END IF
1329: *              end of OK rotation
1330:                            ELSE
1331:                               NOTROT = NOTROT + 1
1332: *[RTD]      SKIPPED  = SKIPPED  + 1
1333:                               PSKIPPED = PSKIPPED + 1
1334:                               IJBLSK = IJBLSK + 1
1335:                            END IF
1336:                         ELSE
1337:                            NOTROT = NOTROT + 1
1338:                            PSKIPPED = PSKIPPED + 1
1339:                            IJBLSK = IJBLSK + 1
1340:                         END IF
1341: *
1342:                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
1343:      +                      THEN
1344:                            SVA( p ) = AAPP
1345:                            NOTROT = 0
1346:                            GO TO 2011
1347:                         END IF
1348:                         IF( ( i.LE.SWBAND ) .AND.
1349:      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1350:                            AAPP = -AAPP
1351:                            NOTROT = 0
1352:                            GO TO 2203
1353:                         END IF
1354: *
1355:  2200                CONTINUE
1356: *        end of the q-loop
1357:  2203                CONTINUE
1358: *
1359:                      SVA( p ) = AAPP
1360: *
1361:                   ELSE
1362: *
1363:                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1364:      +                   MIN0( jgl+KBL-1, N ) - jgl + 1
1365:                      IF( AAPP.LT.ZERO )NOTROT = 0
1366: *
1367:                   END IF
1368: *
1369:  2100          CONTINUE
1370: *     end of the p-loop
1371:  2010       CONTINUE
1372: *     end of the jbc-loop
1373:  2011       CONTINUE
1374: *2011 bailed out of the jbc-loop
1375:             DO 2012 p = igl, MIN0( igl+KBL-1, N )
1376:                SVA( p ) = ABS( SVA( p ) )
1377:  2012       CONTINUE
1378: ***
1379:  2000    CONTINUE
1380: *2000 :: end of the ibr-loop
1381: *
1382: *     .. update SVA(N)
1383:          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
1384:      +       THEN
1385:             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*WORK( N )
1386:          ELSE
1387:             T = ZERO
1388:             AAPP = ZERO
1389:             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
1390:             SVA( N ) = T*SQRT( AAPP )*WORK( N )
1391:          END IF
1392: *
1393: *     Additional steering devices
1394: *
1395:          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
1396:      +       ( ISWROT.LE.N ) ) )SWBAND = i
1397: *
1398:          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
1399:      +       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1400:             GO TO 1994
1401:          END IF
1402: *
1403:          IF( NOTROT.GE.EMPTSW )GO TO 1994
1404: *
1405:  1993 CONTINUE
1406: *     end i=1:NSWEEP loop
1407: *
1408: * #:( Reaching this point means that the procedure has not converged.
1409:       INFO = NSWEEP - 1
1410:       GO TO 1995
1411: *
1412:  1994 CONTINUE
1413: * #:) Reaching this point means numerical convergence after the i-th
1414: *     sweep.
1415: *
1416:       INFO = 0
1417: * #:) INFO = 0 confirms successful iterations.
1418:  1995 CONTINUE
1419: *
1420: *     Sort the singular values and find how many are above
1421: *     the underflow threshold.
1422: *
1423:       N2 = 0
1424:       N4 = 0
1425:       DO 5991 p = 1, N - 1
1426:          q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
1427:          IF( p.NE.q ) THEN
1428:             TEMP1 = SVA( p )
1429:             SVA( p ) = SVA( q )
1430:             SVA( q ) = TEMP1
1431:             TEMP1 = WORK( p )
1432:             WORK( p ) = WORK( q )
1433:             WORK( q ) = TEMP1
1434:             CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
1435:             IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
1436:          END IF
1437:          IF( SVA( p ).NE.ZERO ) THEN
1438:             N4 = N4 + 1
1439:             IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1
1440:          END IF
1441:  5991 CONTINUE
1442:       IF( SVA( N ).NE.ZERO ) THEN
1443:          N4 = N4 + 1
1444:          IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1
1445:       END IF
1446: *
1447: *     Normalize the left singular vectors.
1448: *
1449:       IF( LSVEC .OR. UCTOL ) THEN
1450:          DO 1998 p = 1, N2
1451:             CALL SSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
1452:  1998    CONTINUE
1453:       END IF
1454: *
1455: *     Scale the product of Jacobi rotations (assemble the fast rotations).
1456: *
1457:       IF( RSVEC ) THEN
1458:          IF( APPLV ) THEN
1459:             DO 2398 p = 1, N
1460:                CALL SSCAL( MVL, WORK( p ), V( 1, p ), 1 )
1461:  2398       CONTINUE
1462:          ELSE
1463:             DO 2399 p = 1, N
1464:                TEMP1 = ONE / SNRM2( MVL, V( 1, p ), 1 )
1465:                CALL SSCAL( MVL, TEMP1, V( 1, p ), 1 )
1466:  2399       CONTINUE
1467:          END IF
1468:       END IF
1469: *
1470: *     Undo scaling, if necessary (and possible).
1471:       IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
1472:      +    SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.
1473:      +    ( SFMIN / SCALE ) ) ) ) THEN
1474:          DO 2400 p = 1, N
1475:             SVA( p ) = SCALE*SVA( p )
1476:  2400    CONTINUE
1477:          SCALE = ONE
1478:       END IF
1479: *
1480:       WORK( 1 ) = SCALE
1481: *     The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE
1482: *     then some of the singular values may overflow or underflow and
1483: *     the spectrum is given in this factored representation.
1484: *
1485:       WORK( 2 ) = FLOAT( N4 )
1486: *     N4 is the number of computed nonzero singular values of A.
1487: *
1488:       WORK( 3 ) = FLOAT( N2 )
1489: *     N2 is the number of singular values of A greater than SFMIN.
1490: *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
1491: *     that may carry some information.
1492: *
1493:       WORK( 4 ) = FLOAT( i )
1494: *     i is the index of the last sweep before declaring convergence.
1495: *
1496:       WORK( 5 ) = MXAAPQ
1497: *     MXAAPQ is the largest absolute value of scaled pivots in the
1498: *     last sweep
1499: *
1500:       WORK( 6 ) = MXSINJ
1501: *     MXSINJ is the largest absolute value of the sines of Jacobi angles
1502: *     in the last sweep
1503: *
1504:       RETURN
1505: *     ..
1506: *     .. END OF SGESVJ
1507: *     ..
1508:       END
1509: