SRC/sgesvj.f

   1 *> \brief \b SGESVJ
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGESVJ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgesvj.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgesvj.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgesvj.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
  22 *                          LDV, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
  26 *       CHARACTER*1        JOBA, JOBU, JOBV
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               A( LDA, * ), SVA( N ), V( LDV, * ),
  30 *      $                   WORK( LWORK )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SGESVJ computes the singular value decomposition (SVD) of a real
  40 *> M-by-N matrix A, where M >= N. The SVD of A is written as
  41 *>                                    [++]   [xx]   [x0]   [xx]
  42 *>              A = U * SIGMA * V^t,  [++] = [xx] * [ox] * [xx]
  43 *>                                    [++]   [xx]
  44 *> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
  45 *> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
  46 *> of SIGMA are the singular values of A. The columns of U and V are the
  47 *> left and the right singular vectors of A, respectively.
  48 *> SGESVJ can sometimes compute tiny singular values and their singular vectors much
  49 *> more accurately than other SVD routines, see below under Further Details.
  50 *> \endverbatim
  51 *
  52 *  Arguments:
  53 *  ==========
  54 *
  55 *> \param[in] JOBA
  56 *> \verbatim
  57 *>          JOBA is CHARACTER* 1
  58 *>          Specifies the structure of A.
  59 *>          = 'L': The input matrix A is lower triangular;
  60 *>          = 'U': The input matrix A is upper triangular;
  61 *>          = 'G': The input matrix A is general M-by-N matrix, M >= N.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] JOBU
  65 *> \verbatim
  66 *>          JOBU is CHARACTER*1
  67 *>          Specifies whether to compute the left singular vectors
  68 *>          (columns of U):
  69 *>          = 'U': The left singular vectors corresponding to the nonzero
  70 *>                 singular values are computed and returned in the leading
  71 *>                 columns of A. See more details in the description of A.
  72 *>                 The default numerical orthogonality threshold is set to
  73 *>                 approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E').
  74 *>          = 'C': Analogous to JOBU='U', except that user can control the
  75 *>                 level of numerical orthogonality of the computed left
  76 *>                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
  77 *>                 CTOL is given on input in the array WORK.
  78 *>                 No CTOL smaller than ONE is allowed. CTOL greater
  79 *>                 than 1 / EPS is meaningless. The option 'C'
  80 *>                 can be used if M*EPS is satisfactory orthogonality
  81 *>                 of the computed left singular vectors, so CTOL=M could
  82 *>                 save few sweeps of Jacobi rotations.
  83 *>                 See the descriptions of A and WORK(1).
  84 *>          = 'N': The matrix U is not computed. However, see the
  85 *>                 description of A.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] JOBV
  89 *> \verbatim
  90 *>          JOBV is CHARACTER*1
  91 *>          Specifies whether to compute the right singular vectors, that
  92 *>          is, the matrix V:
  93 *>          = 'V' : the matrix V is computed and returned in the array V
  94 *>          = 'A' : the Jacobi rotations are applied to the MV-by-N
  95 *>                  array V. In other words, the right singular vector
  96 *>                  matrix V is not computed explicitly; instead it is
  97 *>                  applied to an MV-by-N matrix initially stored in the
  98 *>                  first MV rows of V.
  99 *>          = 'N' : the matrix V is not computed and the array V is not
 100 *>                  referenced
 101 *> \endverbatim
 102 *>
 103 *> \param[in] M
 104 *> \verbatim
 105 *>          M is INTEGER
 106 *>          The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] N
 110 *> \verbatim
 111 *>          N is INTEGER
 112 *>          The number of columns of the input matrix A.
 113 *>          M >= N >= 0.
 114 *> \endverbatim
 115 *>
 116 *> \param[in,out] A
 117 *> \verbatim
 118 *>          A is REAL array, dimension (LDA,N)
 119 *>          On entry, the M-by-N matrix A.
 120 *>          On exit,
 121 *>          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C':
 122 *>                 If INFO .EQ. 0 :
 123 *>                 RANKA orthonormal columns of U are returned in the
 124 *>                 leading RANKA columns of the array A. Here RANKA <= N
 125 *>                 is the number of computed singular values of A that are
 126 *>                 above the underflow threshold SLAMCH('S'). The singular
 127 *>                 vectors corresponding to underflowed or zero singular
 128 *>                 values are not computed. The value of RANKA is returned
 129 *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
 130 *>                 descriptions of SVA and WORK. The computed columns of U
 131 *>                 are mutually numerically orthogonal up to approximately
 132 *>                 TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
 133 *>                 see the description of JOBU.
 134 *>                 If INFO .GT. 0,
 135 *>                 the procedure SGESVJ did not converge in the given number
 136 *>                 of iterations (sweeps). In that case, the computed
 137 *>                 columns of U may not be orthogonal up to TOL. The output
 138 *>                 U (stored in A), SIGMA (given by the computed singular
 139 *>                 values in SVA(1:N)) and V is still a decomposition of the
 140 *>                 input matrix A in the sense that the residual
 141 *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
 142 *>          If JOBU .EQ. 'N':
 143 *>                 If INFO .EQ. 0 :
 144 *>                 Note that the left singular vectors are 'for free' in the
 145 *>                 one-sided Jacobi SVD algorithm. However, if only the
 146 *>                 singular values are needed, the level of numerical
 147 *>                 orthogonality of U is not an issue and iterations are
 148 *>                 stopped when the columns of the iterated matrix are
 149 *>                 numerically orthogonal up to approximately M*EPS. Thus,
 150 *>                 on exit, A contains the columns of U scaled with the
 151 *>                 corresponding singular values.
 152 *>                 If INFO .GT. 0 :
 153 *>                 the procedure SGESVJ did not converge in the given number
 154 *>                 of iterations (sweeps).
 155 *> \endverbatim
 156 *>
 157 *> \param[in] LDA
 158 *> \verbatim
 159 *>          LDA is INTEGER
 160 *>          The leading dimension of the array A.  LDA >= max(1,M).
 161 *> \endverbatim
 162 *>
 163 *> \param[out] SVA
 164 *> \verbatim
 165 *>          SVA is REAL array, dimension (N)
 166 *>          On exit,
 167 *>          If INFO .EQ. 0 :
 168 *>          depending on the value SCALE = WORK(1), we have:
 169 *>                 If SCALE .EQ. ONE:
 170 *>                 SVA(1:N) contains the computed singular values of A.
 171 *>                 During the computation SVA contains the Euclidean column
 172 *>                 norms of the iterated matrices in the array A.
 173 *>                 If SCALE .NE. ONE:
 174 *>                 The singular values of A are SCALE*SVA(1:N), and this
 175 *>                 factored representation is due to the fact that some of the
 176 *>                 singular values of A might underflow or overflow.
 177 *>
 178 *>          If INFO .GT. 0 :
 179 *>          the procedure SGESVJ did not converge in the given number of
 180 *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
 181 *> \endverbatim
 182 *>
 183 *> \param[in] MV
 184 *> \verbatim
 185 *>          MV is INTEGER
 186 *>          If JOBV .EQ. 'A', then the product of Jacobi rotations in SGESVJ
 187 *>          is applied to the first MV rows of V. See the description of JOBV.
 188 *> \endverbatim
 189 *>
 190 *> \param[in,out] V
 191 *> \verbatim
 192 *>          V is REAL array, dimension (LDV,N)
 193 *>          If JOBV = 'V', then V contains on exit the N-by-N matrix of
 194 *>                         the right singular vectors;
 195 *>          If JOBV = 'A', then V contains the product of the computed right
 196 *>                         singular vector matrix and the initial matrix in
 197 *>                         the array V.
 198 *>          If JOBV = 'N', then V is not referenced.
 199 *> \endverbatim
 200 *>
 201 *> \param[in] LDV
 202 *> \verbatim
 203 *>          LDV is INTEGER
 204 *>          The leading dimension of the array V, LDV .GE. 1.
 205 *>          If JOBV .EQ. 'V', then LDV .GE. max(1,N).
 206 *>          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
 207 *> \endverbatim
 208 *>
 209 *> \param[in,out] WORK
 210 *> \verbatim
 211 *>          WORK is REAL array, dimension max(4,M+N).
 212 *>          On entry,
 213 *>          If JOBU .EQ. 'C' :
 214 *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
 215 *>                    The process stops if all columns of A are mutually
 216 *>                    orthogonal up to CTOL*EPS, EPS=SLAMCH('E').
 217 *>                    It is required that CTOL >= ONE, i.e. it is not
 218 *>                    allowed to force the routine to obtain orthogonality
 219 *>                    below EPSILON.
 220 *>          On exit,
 221 *>          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
 222 *>                    are the computed singular vcalues of A.
 223 *>                    (See description of SVA().)
 224 *>          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
 225 *>                    singular values.
 226 *>          WORK(3) = NINT(WORK(3)) is the number of the computed singular
 227 *>                    values that are larger than the underflow threshold.
 228 *>          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
 229 *>                    rotations needed for numerical convergence.
 230 *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
 231 *>                    This is useful information in cases when SGESVJ did
 232 *>                    not converge, as it can be used to estimate whether
 233 *>                    the output is stil useful and for post festum analysis.
 234 *>          WORK(6) = the largest absolute value over all sines of the
 235 *>                    Jacobi rotation angles in the last sweep. It can be
 236 *>                    useful for a post festum analysis.
 237 *> \endverbatim
 238 *>
 239 *> \param[in] LWORK
 240 *> \verbatim
 241 *>          LWORK is INTEGER
 242 *>         length of WORK, WORK >= MAX(6,M+N)
 243 *> \endverbatim
 244 *>
 245 *> \param[out] INFO
 246 *> \verbatim
 247 *>          INFO is INTEGER
 248 *>          = 0 : successful exit.
 249 *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
 250 *>          > 0 : SGESVJ did not converge in the maximal allowed number (30)
 251 *>                of sweeps. The output may still be useful. See the
 252 *>                description of WORK.
 253 *> \endverbatim
 254 *
 255 *  Authors:
 256 *  ========
 257 *
 258 *> \author Univ. of Tennessee
 259 *> \author Univ. of California Berkeley
 260 *> \author Univ. of Colorado Denver
 261 *> \author NAG Ltd.
 262 *
 263 *> \date November 2015
 264 *
 265 *> \ingroup realGEcomputational
 266 *
 267 *> \par Further Details:
 268 *  =====================
 269 *>
 270 *> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
 271 *> rotations. The rotations are implemented as fast scaled rotations of
 272 *> Anda and Park [1]. In the case of underflow of the Jacobi angle, a
 273 *> modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
 274 *> column interchanges of de Rijk [2]. The relative accuracy of the computed
 275 *> singular values and the accuracy of the computed singular vectors (in
 276 *> angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
 277 *> The condition number that determines the accuracy in the full rank case
 278 *> is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
 279 *> spectral condition number. The best performance of this Jacobi SVD
 280 *> procedure is achieved if used in an  accelerated version of Drmac and
 281 *> Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
 282 *> Some tunning parameters (marked with [TP]) are available for the
 283 *> implementer. \n
 284 *> The computational range for the nonzero singular values is the  machine
 285 *> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
 286 *> denormalized singular values can be computed with the corresponding
 287 *> gradual loss of accurate digits.
 288 *>
 289 *> \par Contributors:
 290 *  ==================
 291 *>
 292 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
 293 *>
 294 *> \par References:
 295 *  ================
 296 *>
 297 *> [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. \n
 298 *>    SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. \n\n
 299 *> [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
 300 *>    singular value decomposition on a vector computer. \n
 301 *>    SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. \n\n
 302 *> [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. \n
 303 *> [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
 304 *>    value computation in floating point arithmetic. \n
 305 *>    SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. \n\n
 306 *> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. \n
 307 *>    SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. \n
 308 *>    LAPACK Working note 169. \n\n
 309 *> [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. \n
 310 *>    SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. \n
 311 *>    LAPACK Working note 170. \n\n
 312 *> [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
 313 *>    QSVD, (H,K)-SVD computations.\n
 314 *>    Department of Mathematics, University of Zagreb, 2008.
 315 *>
 316 *> \par Bugs, Examples and Comments:
 317 *  =================================
 318 *>
 319 *> Please report all bugs and send interesting test examples and comments to
 320 *> drmac@math.hr. Thank you.
 321 *
 322 *  =====================================================================
 323       SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
 324      $                   LDV, WORK, LWORK, INFO )
 325 *
 326 *  -- LAPACK computational routine (version 3.6.0) --
 327 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 328 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 329 *     November 2015
 330 *
 331 *     .. Scalar Arguments ..
 332       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
 333       CHARACTER*1        JOBA, JOBU, JOBV
 334 *     ..
 335 *     .. Array Arguments ..
 336       REAL               A( LDA, * ), SVA( N ), V( LDV, * ),
 337      $                   WORK( LWORK )
 338 *     ..
 339 *
 340 *  =====================================================================
 341 *
 342 *     .. Local Parameters ..
 343       REAL               ZERO, HALF, ONE
 344       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0)
 345       INTEGER            NSWEEP
 346       PARAMETER          ( NSWEEP = 30 )
 347 *     ..
 348 *     .. Local Scalars ..
 349       REAL               AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
 350      $                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
 351      $                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
 352      $                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
 353      $                   THSIGN, TOL
 354       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
 355      $                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
 356      $                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
 357      $                   SWBAND
 358       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
 359      $                   RSVEC, UCTOL, UPPER
 360 *     ..
 361 *     .. Local Arrays ..
 362       REAL               FASTR( 5 )
 363 *     ..
 364 *     .. Intrinsic Functions ..
 365       INTRINSIC          ABS, MAX, MIN, FLOAT, SIGN, SQRT
 366 *     ..
 367 *     .. External Functions ..
 368 *     ..
 369 *     from BLAS
 370       REAL               SDOT, SNRM2
 371       EXTERNAL           SDOT, SNRM2
 372       INTEGER            ISAMAX
 373       EXTERNAL           ISAMAX
 374 *     from LAPACK
 375       REAL               SLAMCH
 376       EXTERNAL           SLAMCH
 377       LOGICAL            LSAME
 378       EXTERNAL           LSAME
 379 *     ..
 380 *     .. External Subroutines ..
 381 *     ..
 382 *     from BLAS
 383       EXTERNAL           SAXPY, SCOPY, SROTM, SSCAL, SSWAP
 384 *     from LAPACK
 385       EXTERNAL           SLASCL, SLASET, SLASSQ, XERBLA
 386 *
 387       EXTERNAL           SGSVJ0, SGSVJ1
 388 *     ..
 389 *     .. Executable Statements ..
 390 *
 391 *     Test the input arguments
 392 *
 393       LSVEC = LSAME( JOBU, 'U' )
 394       UCTOL = LSAME( JOBU, 'C' )
 395       RSVEC = LSAME( JOBV, 'V' )
 396       APPLV = LSAME( JOBV, 'A' )
 397       UPPER = LSAME( JOBA, 'U' )
 398       LOWER = LSAME( JOBA, 'L' )
 399 *
 400       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
 401          INFO = -1
 402       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
 403          INFO = -2
 404       ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 405          INFO = -3
 406       ELSE IF( M.LT.0 ) THEN
 407          INFO = -4
 408       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
 409          INFO = -5
 410       ELSE IF( LDA.LT.M ) THEN
 411          INFO = -7
 412       ELSE IF( MV.LT.0 ) THEN
 413          INFO = -9
 414       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
 415      $         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
 416          INFO = -11
 417       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
 418          INFO = -12
 419       ELSE IF( LWORK.LT.MAX( M+N, 6 ) ) THEN
 420          INFO = -13
 421       ELSE
 422          INFO = 0
 423       END IF
 424 *
 425 *     #:(
 426       IF( INFO.NE.0 ) THEN
 427          CALL XERBLA( 'SGESVJ', -INFO )
 428          RETURN
 429       END IF
 430 *
 431 * #:) Quick return for void matrix
 432 *
 433       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
 434 *
 435 *     Set numerical parameters
 436 *     The stopping criterion for Jacobi rotations is
 437 *
 438 *     max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
 439 *
 440 *     where EPS is the round-off and CTOL is defined as follows:
 441 *
 442       IF( UCTOL ) THEN
 443 *        ... user controlled
 444          CTOL = WORK( 1 )
 445       ELSE
 446 *        ... default
 447          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
 448             CTOL = SQRT( FLOAT( M ) )
 449          ELSE
 450             CTOL = FLOAT( M )
 451          END IF
 452       END IF
 453 *     ... and the machine dependent parameters are
 454 *[!]  (Make sure that SLAMCH() works properly on the target machine.)
 455 *
 456       EPSLN = SLAMCH( 'Epsilon' )
 457       ROOTEPS = SQRT( EPSLN )
 458       SFMIN = SLAMCH( 'SafeMinimum' )
 459       ROOTSFMIN = SQRT( SFMIN )
 460       SMALL = SFMIN / EPSLN
 461       BIG = SLAMCH( 'Overflow' )
 462 *     BIG         = ONE    / SFMIN
 463       ROOTBIG = ONE / ROOTSFMIN
 464       LARGE = BIG / SQRT( FLOAT( M*N ) )
 465       BIGTHETA = ONE / ROOTEPS
 466 *
 467       TOL = CTOL*EPSLN
 468       ROOTTOL = SQRT( TOL )
 469 *
 470       IF( FLOAT( M )*EPSLN.GE.ONE ) THEN
 471          INFO = -4
 472          CALL XERBLA( 'SGESVJ', -INFO )
 473          RETURN
 474       END IF
 475 *
 476 *     Initialize the right singular vector matrix.
 477 *
 478       IF( RSVEC ) THEN
 479          MVL = N
 480          CALL SLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
 481       ELSE IF( APPLV ) THEN
 482          MVL = MV
 483       END IF
 484       RSVEC = RSVEC .OR. APPLV
 485 *
 486 *     Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
 487 *(!)  If necessary, scale A to protect the largest singular value
 488 *     from overflow. It is possible that saving the largest singular
 489 *     value destroys the information about the small ones.
 490 *     This initial scaling is almost minimal in the sense that the
 491 *     goal is to make sure that no column norm overflows, and that
 492 *     SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 493 *     in A are detected, the procedure returns with INFO=-6.
 494 *
 495       SKL = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
 496       NOSCALE = .TRUE.
 497       GOSCALE = .TRUE.
 498 *
 499       IF( LOWER ) THEN
 500 *        the input matrix is M-by-N lower triangular (trapezoidal)
 501          DO 1874 p = 1, N
 502             AAPP = ZERO
 503             AAQQ = ONE
 504             CALL SLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
 505             IF( AAPP.GT.BIG ) THEN
 506                INFO = -6
 507                CALL XERBLA( 'SGESVJ', -INFO )
 508                RETURN
 509             END IF
 510             AAQQ = SQRT( AAQQ )
 511             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
 512                SVA( p ) = AAPP*AAQQ
 513             ELSE
 514                NOSCALE = .FALSE.
 515                SVA( p ) = AAPP*( AAQQ*SKL )
 516                IF( GOSCALE ) THEN
 517                   GOSCALE = .FALSE.
 518                   DO 1873 q = 1, p - 1
 519                      SVA( q ) = SVA( q )*SKL
 520  1873             CONTINUE
 521                END IF
 522             END IF
 523  1874    CONTINUE
 524       ELSE IF( UPPER ) THEN
 525 *        the input matrix is M-by-N upper triangular (trapezoidal)
 526          DO 2874 p = 1, N
 527             AAPP = ZERO
 528             AAQQ = ONE
 529             CALL SLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
 530             IF( AAPP.GT.BIG ) THEN
 531                INFO = -6
 532                CALL XERBLA( 'SGESVJ', -INFO )
 533                RETURN
 534             END IF
 535             AAQQ = SQRT( AAQQ )
 536             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
 537                SVA( p ) = AAPP*AAQQ
 538             ELSE
 539                NOSCALE = .FALSE.
 540                SVA( p ) = AAPP*( AAQQ*SKL )
 541                IF( GOSCALE ) THEN
 542                   GOSCALE = .FALSE.
 543                   DO 2873 q = 1, p - 1
 544                      SVA( q ) = SVA( q )*SKL
 545  2873             CONTINUE
 546                END IF
 547             END IF
 548  2874    CONTINUE
 549       ELSE
 550 *        the input matrix is M-by-N general dense
 551          DO 3874 p = 1, N
 552             AAPP = ZERO
 553             AAQQ = ONE
 554             CALL SLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
 555             IF( AAPP.GT.BIG ) THEN
 556                INFO = -6
 557                CALL XERBLA( 'SGESVJ', -INFO )
 558                RETURN
 559             END IF
 560             AAQQ = SQRT( AAQQ )
 561             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
 562                SVA( p ) = AAPP*AAQQ
 563             ELSE
 564                NOSCALE = .FALSE.
 565                SVA( p ) = AAPP*( AAQQ*SKL )
 566                IF( GOSCALE ) THEN
 567                   GOSCALE = .FALSE.
 568                   DO 3873 q = 1, p - 1
 569                      SVA( q ) = SVA( q )*SKL
 570  3873             CONTINUE
 571                END IF
 572             END IF
 573  3874    CONTINUE
 574       END IF
 575 *
 576       IF( NOSCALE )SKL = ONE
 577 *
 578 *     Move the smaller part of the spectrum from the underflow threshold
 579 *(!)  Start by determining the position of the nonzero entries of the
 580 *     array SVA() relative to ( SFMIN, BIG ).
 581 *
 582       AAPP = ZERO
 583       AAQQ = BIG
 584       DO 4781 p = 1, N
 585          IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) )
 586          AAPP = MAX( AAPP, SVA( p ) )
 587  4781 CONTINUE
 588 *
 589 * #:) Quick return for zero matrix
 590 *
 591       IF( AAPP.EQ.ZERO ) THEN
 592          IF( LSVEC )CALL SLASET( 'G', M, N, ZERO, ONE, A, LDA )
 593          WORK( 1 ) = ONE
 594          WORK( 2 ) = ZERO
 595          WORK( 3 ) = ZERO
 596          WORK( 4 ) = ZERO
 597          WORK( 5 ) = ZERO
 598          WORK( 6 ) = ZERO
 599          RETURN
 600       END IF
 601 *
 602 * #:) Quick return for one-column matrix
 603 *
 604       IF( N.EQ.1 ) THEN
 605          IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
 606      $                           A( 1, 1 ), LDA, IERR )
 607          WORK( 1 ) = ONE / SKL
 608          IF( SVA( 1 ).GE.SFMIN ) THEN
 609             WORK( 2 ) = ONE
 610          ELSE
 611             WORK( 2 ) = ZERO
 612          END IF
 613          WORK( 3 ) = ZERO
 614          WORK( 4 ) = ZERO
 615          WORK( 5 ) = ZERO
 616          WORK( 6 ) = ZERO
 617          RETURN
 618       END IF
 619 *
 620 *     Protect small singular values from underflow, and try to
 621 *     avoid underflows/overflows in computing Jacobi rotations.
 622 *
 623       SN = SQRT( SFMIN / EPSLN )
 624       TEMP1 = SQRT( BIG / FLOAT( N ) )
 625       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
 626      $    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
 627          TEMP1 = MIN( BIG, TEMP1 / AAPP )
 628 *         AAQQ  = AAQQ*TEMP1
 629 *         AAPP  = AAPP*TEMP1
 630       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
 631          TEMP1 = MIN( SN / AAQQ, BIG / ( AAPP*SQRT( FLOAT( N ) ) ) )
 632 *         AAQQ  = AAQQ*TEMP1
 633 *         AAPP  = AAPP*TEMP1
 634       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
 635          TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP )
 636 *         AAQQ  = AAQQ*TEMP1
 637 *         AAPP  = AAPP*TEMP1
 638       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
 639          TEMP1 = MIN( SN / AAQQ, BIG / ( SQRT( FLOAT( N ) )*AAPP ) )
 640 *         AAQQ  = AAQQ*TEMP1
 641 *         AAPP  = AAPP*TEMP1
 642       ELSE
 643          TEMP1 = ONE
 644       END IF
 645 *
 646 *     Scale, if necessary
 647 *
 648       IF( TEMP1.NE.ONE ) THEN
 649          CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
 650       END IF
 651       SKL = TEMP1*SKL
 652       IF( SKL.NE.ONE ) THEN
 653          CALL SLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
 654          SKL = ONE / SKL
 655       END IF
 656 *
 657 *     Row-cyclic Jacobi SVD algorithm with column pivoting
 658 *
 659       EMPTSW = ( N*( N-1 ) ) / 2
 660       NOTROT = 0
 661       FASTR( 1 ) = ZERO
 662 *
 663 *     A is represented in factored form A = A * diag(WORK), where diag(WORK)
 664 *     is initialized to identity. WORK is updated during fast scaled
 665 *     rotations.
 666 *
 667       DO 1868 q = 1, N
 668          WORK( q ) = ONE
 669  1868 CONTINUE
 670 *
 671 *
 672       SWBAND = 3
 673 *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
 674 *     if SGESVJ is used as a computational routine in the preconditioned
 675 *     Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
 676 *     works on pivots inside a band-like region around the diagonal.
 677 *     The boundaries are determined dynamically, based on the number of
 678 *     pivots above a threshold.
 679 *
 680       KBL = MIN( 8, N )
 681 *[TP] KBL is a tuning parameter that defines the tile size in the
 682 *     tiling of the p-q loops of pivot pairs. In general, an optimal
 683 *     value of KBL depends on the matrix dimensions and on the
 684 *     parameters of the computer's memory.
 685 *
 686       NBL = N / KBL
 687       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
 688 *
 689       BLSKIP = KBL**2
 690 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 691 *
 692       ROWSKIP = MIN( 5, KBL )
 693 *[TP] ROWSKIP is a tuning parameter.
 694 *
 695       LKAHEAD = 1
 696 *[TP] LKAHEAD is a tuning parameter.
 697 *
 698 *     Quasi block transformations, using the lower (upper) triangular
 699 *     structure of the input matrix. The quasi-block-cycling usually
 700 *     invokes cubic convergence. Big part of this cycle is done inside
 701 *     canonical subspaces of dimensions less than M.
 702 *
 703       IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN
 704 *[TP] The number of partition levels and the actual partition are
 705 *     tuning parameters.
 706          N4 = N / 4
 707          N2 = N / 2
 708          N34 = 3*N4
 709          IF( APPLV ) THEN
 710             q = 0
 711          ELSE
 712             q = 1
 713          END IF
 714 *
 715          IF( LOWER ) THEN
 716 *
 717 *     This works very well on lower triangular matrices, in particular
 718 *     in the framework of the preconditioned Jacobi SVD (xGEJSV).
 719 *     The idea is simple:
 720 *     [+ 0 0 0]   Note that Jacobi transformations of [0 0]
 721 *     [+ + 0 0]                                       [0 0]
 722 *     [+ + x 0]   actually work on [x 0]              [x 0]
 723 *     [+ + x x]                    [x x].             [x x]
 724 *
 725             CALL SGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
 726      $                   WORK( N34+1 ), SVA( N34+1 ), MVL,
 727      $                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
 728      $                   2, WORK( N+1 ), LWORK-N, IERR )
 729 *
 730             CALL SGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
 731      $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
 732      $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
 733      $                   WORK( N+1 ), LWORK-N, IERR )
 734 *
 735             CALL SGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
 736      $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
 737      $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
 738      $                   WORK( N+1 ), LWORK-N, IERR )
 739 *
 740             CALL SGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
 741      $                   WORK( N4+1 ), SVA( N4+1 ), MVL,
 742      $                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
 743      $                   WORK( N+1 ), LWORK-N, IERR )
 744 *
 745             CALL SGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
 746      $                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
 747      $                   IERR )
 748 *
 749             CALL SGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
 750      $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
 751      $                   LWORK-N, IERR )
 752 *
 753 *
 754          ELSE IF( UPPER ) THEN
 755 *
 756 *
 757             CALL SGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
 758      $                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
 759      $                   IERR )
 760 *
 761             CALL SGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
 762      $                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
 763      $                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
 764      $                   IERR )
 765 *
 766             CALL SGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
 767      $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
 768      $                   LWORK-N, IERR )
 769 *
 770             CALL SGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
 771      $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
 772      $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
 773      $                   WORK( N+1 ), LWORK-N, IERR )
 774
 775          END IF
 776 *
 777       END IF
 778 *
 779 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 780 *
 781       DO 1993 i = 1, NSWEEP
 782 *
 783 *     .. go go go ...
 784 *
 785          MXAAPQ = ZERO
 786          MXSINJ = ZERO
 787          ISWROT = 0
 788 *
 789          NOTROT = 0
 790          PSKIPPED = 0
 791 *
 792 *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
 793 *     1 <= p < q <= N. This is the first step toward a blocked implementation
 794 *     of the rotations. New implementation, based on block transformations,
 795 *     is under development.
 796 *
 797          DO 2000 ibr = 1, NBL
 798 *
 799             igl = ( ibr-1 )*KBL + 1
 800 *
 801             DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
 802 *
 803                igl = igl + ir1*KBL
 804 *
 805                DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
 806 *
 807 *     .. de Rijk's pivoting
 808 *
 809                   q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
 810                   IF( p.NE.q ) THEN
 811                      CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
 812                      IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1,
 813      $                                      V( 1, q ), 1 )
 814                      TEMP1 = SVA( p )
 815                      SVA( p ) = SVA( q )
 816                      SVA( q ) = TEMP1
 817                      TEMP1 = WORK( p )
 818                      WORK( p ) = WORK( q )
 819                      WORK( q ) = TEMP1
 820                   END IF
 821 *
 822                   IF( ir1.EQ.0 ) THEN
 823 *
 824 *        Column norms are periodically updated by explicit
 825 *        norm computation.
 826 *        Caveat:
 827 *        Unfortunately, some BLAS implementations compute SNRM2(M,A(1,p),1)
 828 *        as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
 829 *        overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
 830 *        underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
 831 *        Hence, SNRM2 cannot be trusted, not even in the case when
 832 *        the true norm is far from the under(over)flow boundaries.
 833 *        If properly implemented SNRM2 is available, the IF-THEN-ELSE
 834 *        below should read "AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)".
 835 *
 836                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
 837      $                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
 838                         SVA( p ) = SNRM2( M, A( 1, p ), 1 )*WORK( p )
 839                      ELSE
 840                         TEMP1 = ZERO
 841                         AAPP = ONE
 842                         CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
 843                         SVA( p ) = TEMP1*SQRT( AAPP )*WORK( p )
 844                      END IF
 845                      AAPP = SVA( p )
 846                   ELSE
 847                      AAPP = SVA( p )
 848                   END IF
 849 *
 850                   IF( AAPP.GT.ZERO ) THEN
 851 *
 852                      PSKIPPED = 0
 853 *
 854                      DO 2002 q = p + 1, MIN( igl+KBL-1, N )
 855 *
 856                         AAQQ = SVA( q )
 857 *
 858                         IF( AAQQ.GT.ZERO ) THEN
 859 *
 860                            AAPP0 = AAPP
 861                            IF( AAQQ.GE.ONE ) THEN
 862                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
 863                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
 864                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
 865      $                                  q ), 1 )*WORK( p )*WORK( q ) /
 866      $                                  AAQQ ) / AAPP
 867                               ELSE
 868                                  CALL SCOPY( M, A( 1, p ), 1,
 869      $                                       WORK( N+1 ), 1 )
 870                                  CALL SLASCL( 'G', 0, 0, AAPP,
 871      $                                        WORK( p ), M, 1,
 872      $                                        WORK( N+1 ), LDA, IERR )
 873                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
 874      $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
 875                               END IF
 876                            ELSE
 877                               ROTOK = AAPP.LE.( AAQQ / SMALL )
 878                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
 879                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
 880      $                                  q ), 1 )*WORK( p )*WORK( q ) /
 881      $                                  AAQQ ) / AAPP
 882                               ELSE
 883                                  CALL SCOPY( M, A( 1, q ), 1,
 884      $                                       WORK( N+1 ), 1 )
 885                                  CALL SLASCL( 'G', 0, 0, AAQQ,
 886      $                                        WORK( q ), M, 1,
 887      $                                        WORK( N+1 ), LDA, IERR )
 888                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
 889      $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
 890                               END IF
 891                            END IF
 892 *
 893                            MXAAPQ = MAX( MXAAPQ, ABS( AAPQ ) )
 894 *
 895 *        TO rotate or NOT to rotate, THAT is the question ...
 896 *
 897                            IF( ABS( AAPQ ).GT.TOL ) THEN
 898 *
 899 *           .. rotate
 900 *[RTD]      ROTATED = ROTATED + ONE
 901 *
 902                               IF( ir1.EQ.0 ) THEN
 903                                  NOTROT = 0
 904                                  PSKIPPED = 0
 905                                  ISWROT = ISWROT + 1
 906                               END IF
 907 *
 908                               IF( ROTOK ) THEN
 909 *
 910                                  AQOAP = AAQQ / AAPP
 911                                  APOAQ = AAPP / AAQQ
 912                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
 913 *
 914                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
 915 *
 916                                     T = HALF / THETA
 917                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
 918                                     FASTR( 4 ) = -T*WORK( q ) /
 919      $                                           WORK( p )
 920                                     CALL SROTM( M, A( 1, p ), 1,
 921      $                                          A( 1, q ), 1, FASTR )
 922                                     IF( RSVEC )CALL SROTM( MVL,
 923      $                                              V( 1, p ), 1,
 924      $                                              V( 1, q ), 1,
 925      $                                              FASTR )
 926                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
 927      $                                         ONE+T*APOAQ*AAPQ ) )
 928                                     AAPP = AAPP*SQRT( MAX( ZERO,
 929      $                                         ONE-T*AQOAP*AAPQ ) )
 930                                     MXSINJ = MAX( MXSINJ, ABS( T ) )
 931 *
 932                                  ELSE
 933 *
 934 *                 .. choose correct signum for THETA and rotate
 935 *
 936                                     THSIGN = -SIGN( ONE, AAPQ )
 937                                     T = ONE / ( THETA+THSIGN*
 938      $                                  SQRT( ONE+THETA*THETA ) )
 939                                     CS = SQRT( ONE / ( ONE+T*T ) )
 940                                     SN = T*CS
 941 *
 942                                     MXSINJ = MAX( MXSINJ, ABS( SN ) )
 943                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
 944      $                                         ONE+T*APOAQ*AAPQ ) )
 945                                     AAPP = AAPP*SQRT( MAX( ZERO,
 946      $                                     ONE-T*AQOAP*AAPQ ) )
 947 *
 948                                     APOAQ = WORK( p ) / WORK( q )
 949                                     AQOAP = WORK( q ) / WORK( p )
 950                                     IF( WORK( p ).GE.ONE ) THEN
 951                                        IF( WORK( q ).GE.ONE ) THEN
 952                                           FASTR( 3 ) = T*APOAQ
 953                                           FASTR( 4 ) = -T*AQOAP
 954                                           WORK( p ) = WORK( p )*CS
 955                                           WORK( q ) = WORK( q )*CS
 956                                           CALL SROTM( M, A( 1, p ), 1,
 957      $                                                A( 1, q ), 1,
 958      $                                                FASTR )
 959                                           IF( RSVEC )CALL SROTM( MVL,
 960      $                                        V( 1, p ), 1, V( 1, q ),
 961      $                                        1, FASTR )
 962                                        ELSE
 963                                           CALL SAXPY( M, -T*AQOAP,
 964      $                                                A( 1, q ), 1,
 965      $                                                A( 1, p ), 1 )
 966                                           CALL SAXPY( M, CS*SN*APOAQ,
 967      $                                                A( 1, p ), 1,
 968      $                                                A( 1, q ), 1 )
 969                                           WORK( p ) = WORK( p )*CS
 970                                           WORK( q ) = WORK( q ) / CS
 971                                           IF( RSVEC ) THEN
 972                                              CALL SAXPY( MVL, -T*AQOAP,
 973      $                                                   V( 1, q ), 1,
 974      $                                                   V( 1, p ), 1 )
 975                                              CALL SAXPY( MVL,
 976      $                                                   CS*SN*APOAQ,
 977      $                                                   V( 1, p ), 1,
 978      $                                                   V( 1, q ), 1 )
 979                                           END IF
 980                                        END IF
 981                                     ELSE
 982                                        IF( WORK( q ).GE.ONE ) THEN
 983                                           CALL SAXPY( M, T*APOAQ,
 984      $                                                A( 1, p ), 1,
 985      $                                                A( 1, q ), 1 )
 986                                           CALL SAXPY( M, -CS*SN*AQOAP,
 987      $                                                A( 1, q ), 1,
 988      $                                                A( 1, p ), 1 )
 989                                           WORK( p ) = WORK( p ) / CS
 990                                           WORK( q ) = WORK( q )*CS
 991                                           IF( RSVEC ) THEN
 992                                              CALL SAXPY( MVL, T*APOAQ,
 993      $                                                   V( 1, p ), 1,
 994      $                                                   V( 1, q ), 1 )
 995                                              CALL SAXPY( MVL,
 996      $                                                   -CS*SN*AQOAP,
 997      $                                                   V( 1, q ), 1,
 998      $                                                   V( 1, p ), 1 )
 999                                           END IF
1000                                        ELSE
1001                                           IF( WORK( p ).GE.WORK( q ) )
1002      $                                        THEN
1003                                              CALL SAXPY( M, -T*AQOAP,
1004      $                                                   A( 1, q ), 1,
1005      $                                                   A( 1, p ), 1 )
1006                                              CALL SAXPY( M, CS*SN*APOAQ,
1007      $                                                   A( 1, p ), 1,
1008      $                                                   A( 1, q ), 1 )
1009                                              WORK( p ) = WORK( p )*CS
1010                                              WORK( q ) = WORK( q ) / CS
1011                                              IF( RSVEC ) THEN
1012                                                 CALL SAXPY( MVL,
1013      $                                               -T*AQOAP,
1014      $                                               V( 1, q ), 1,
1015      $                                               V( 1, p ), 1 )
1016                                                 CALL SAXPY( MVL,
1017      $                                               CS*SN*APOAQ,
1018      $                                               V( 1, p ), 1,
1019      $                                               V( 1, q ), 1 )
1020                                              END IF
1021                                           ELSE
1022                                              CALL SAXPY( M, T*APOAQ,
1023      $                                                   A( 1, p ), 1,
1024      $                                                   A( 1, q ), 1 )
1025                                              CALL SAXPY( M,
1026      $                                                   -CS*SN*AQOAP,
1027      $                                                   A( 1, q ), 1,
1028      $                                                   A( 1, p ), 1 )
1029                                              WORK( p ) = WORK( p ) / CS
1030                                              WORK( q ) = WORK( q )*CS
1031                                              IF( RSVEC ) THEN
1032                                                 CALL SAXPY( MVL,
1033      $                                               T*APOAQ, V( 1, p ),
1034      $                                               1, V( 1, q ), 1 )
1035                                                 CALL SAXPY( MVL,
1036      $                                               -CS*SN*AQOAP,
1037      $                                               V( 1, q ), 1,
1038      $                                               V( 1, p ), 1 )
1039                                              END IF
1040                                           END IF
1041                                        END IF
1042                                     END IF
1043                                  END IF
1044 *
1045                               ELSE
1046 *              .. have to use modified Gram-Schmidt like transformation
1047                                  CALL SCOPY( M, A( 1, p ), 1,
1048      $                                       WORK( N+1 ), 1 )
1049                                  CALL SLASCL( 'G', 0, 0, AAPP, ONE, M,
1050      $                                        1, WORK( N+1 ), LDA,
1051      $                                        IERR )
1052                                  CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M,
1053      $                                        1, A( 1, q ), LDA, IERR )
1054                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )
1055                                  CALL SAXPY( M, TEMP1, WORK( N+1 ), 1,
1056      $                                       A( 1, q ), 1 )
1057                                  CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M,
1058      $                                        1, A( 1, q ), LDA, IERR )
1059                                  SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1060      $                                      ONE-AAPQ*AAPQ ) )
1061                                  MXSINJ = MAX( MXSINJ, SFMIN )
1062                               END IF
1063 *           END IF ROTOK THEN ... ELSE
1064 *
1065 *           In the case of cancellation in updating SVA(q), SVA(p)
1066 *           recompute SVA(q), SVA(p).
1067 *
1068                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1069      $                            THEN
1070                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
1071      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
1072                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
1073      $                                         WORK( q )
1074                                  ELSE
1075                                     T = ZERO
1076                                     AAQQ = ONE
1077                                     CALL SLASSQ( M, A( 1, q ), 1, T,
1078      $                                           AAQQ )
1079                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
1080                                  END IF
1081                               END IF
1082                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
1083                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
1084      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1085                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
1086      $                                     WORK( p )
1087                                  ELSE
1088                                     T = ZERO
1089                                     AAPP = ONE
1090                                     CALL SLASSQ( M, A( 1, p ), 1, T,
1091      $                                           AAPP )
1092                                     AAPP = T*SQRT( AAPP )*WORK( p )
1093                                  END IF
1094                                  SVA( p ) = AAPP
1095                               END IF
1096 *
1097                            ELSE
1098 *        A(:,p) and A(:,q) already numerically orthogonal
1099                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1100 *[RTD]      SKIPPED  = SKIPPED  + 1
1101                               PSKIPPED = PSKIPPED + 1
1102                            END IF
1103                         ELSE
1104 *        A(:,q) is zero column
1105                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
1106                            PSKIPPED = PSKIPPED + 1
1107                         END IF
1108 *
1109                         IF( ( i.LE.SWBAND ) .AND.
1110      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1111                            IF( ir1.EQ.0 )AAPP = -AAPP
1112                            NOTROT = 0
1113                            GO TO 2103
1114                         END IF
1115 *
1116  2002                CONTINUE
1117 *     END q-LOOP
1118 *
1119  2103                CONTINUE
1120 *     bailed out of q-loop
1121 *
1122                      SVA( p ) = AAPP
1123 *
1124                   ELSE
1125                      SVA( p ) = AAPP
1126                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
1127      $                   NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p
1128                   END IF
1129 *
1130  2001          CONTINUE
1131 *     end of the p-loop
1132 *     end of doing the block ( ibr, ibr )
1133  1002       CONTINUE
1134 *     end of ir1-loop
1135 *
1136 * ... go to the off diagonal blocks
1137 *
1138             igl = ( ibr-1 )*KBL + 1
1139 *
1140             DO 2010 jbc = ibr + 1, NBL
1141 *
1142                jgl = ( jbc-1 )*KBL + 1
1143 *
1144 *        doing the block at ( ibr, jbc )
1145 *
1146                IJBLSK = 0
1147                DO 2100 p = igl, MIN( igl+KBL-1, N )
1148 *
1149                   AAPP = SVA( p )
1150                   IF( AAPP.GT.ZERO ) THEN
1151 *
1152                      PSKIPPED = 0
1153 *
1154                      DO 2200 q = jgl, MIN( jgl+KBL-1, N )
1155 *
1156                         AAQQ = SVA( q )
1157                         IF( AAQQ.GT.ZERO ) THEN
1158                            AAPP0 = AAPP
1159 *
1160 *     .. M x 2 Jacobi SVD ..
1161 *
1162 *        Safe Gram matrix computation
1163 *
1164                            IF( AAQQ.GE.ONE ) THEN
1165                               IF( AAPP.GE.AAQQ ) THEN
1166                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
1167                               ELSE
1168                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
1169                               END IF
1170                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
1171                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
1172      $                                  q ), 1 )*WORK( p )*WORK( q ) /
1173      $                                  AAQQ ) / AAPP
1174                               ELSE
1175                                  CALL SCOPY( M, A( 1, p ), 1,
1176      $                                       WORK( N+1 ), 1 )
1177                                  CALL SLASCL( 'G', 0, 0, AAPP,
1178      $                                        WORK( p ), M, 1,
1179      $                                        WORK( N+1 ), LDA, IERR )
1180                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
1181      $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
1182                               END IF
1183                            ELSE
1184                               IF( AAPP.GE.AAQQ ) THEN
1185                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
1186                               ELSE
1187                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
1188                               END IF
1189                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
1190                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
1191      $                                  q ), 1 )*WORK( p )*WORK( q ) /
1192      $                                  AAQQ ) / AAPP
1193                               ELSE
1194                                  CALL SCOPY( M, A( 1, q ), 1,
1195      $                                       WORK( N+1 ), 1 )
1196                                  CALL SLASCL( 'G', 0, 0, AAQQ,
1197      $                                        WORK( q ), M, 1,
1198      $                                        WORK( N+1 ), LDA, IERR )
1199                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
1200      $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
1201                               END IF
1202                            END IF
1203 *
1204                            MXAAPQ = MAX( MXAAPQ, ABS( AAPQ ) )
1205 *
1206 *        TO rotate or NOT to rotate, THAT is the question ...
1207 *
1208                            IF( ABS( AAPQ ).GT.TOL ) THEN
1209                               NOTROT = 0
1210 *[RTD]      ROTATED  = ROTATED + 1
1211                               PSKIPPED = 0
1212                               ISWROT = ISWROT + 1
1213 *
1214                               IF( ROTOK ) THEN
1215 *
1216                                  AQOAP = AAQQ / AAPP
1217                                  APOAQ = AAPP / AAQQ
1218                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
1219                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
1220 *
1221                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
1222                                     T = HALF / THETA
1223                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
1224                                     FASTR( 4 ) = -T*WORK( q ) /
1225      $                                           WORK( p )
1226                                     CALL SROTM( M, A( 1, p ), 1,
1227      $                                          A( 1, q ), 1, FASTR )
1228                                     IF( RSVEC )CALL SROTM( MVL,
1229      $                                              V( 1, p ), 1,
1230      $                                              V( 1, q ), 1,
1231      $                                              FASTR )
1232                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1233      $                                         ONE+T*APOAQ*AAPQ ) )
1234                                     AAPP = AAPP*SQRT( MAX( ZERO,
1235      $                                     ONE-T*AQOAP*AAPQ ) )
1236                                     MXSINJ = MAX( MXSINJ, ABS( T ) )
1237                                  ELSE
1238 *
1239 *                 .. choose correct signum for THETA and rotate
1240 *
1241                                     THSIGN = -SIGN( ONE, AAPQ )
1242                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
1243                                     T = ONE / ( THETA+THSIGN*
1244      $                                  SQRT( ONE+THETA*THETA ) )
1245                                     CS = SQRT( ONE / ( ONE+T*T ) )
1246                                     SN = T*CS
1247                                     MXSINJ = MAX( MXSINJ, ABS( SN ) )
1248                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1249      $                                         ONE+T*APOAQ*AAPQ ) )
1250                                     AAPP = AAPP*SQRT( MAX( ZERO,
1251      $                                         ONE-T*AQOAP*AAPQ ) )
1252 *
1253                                     APOAQ = WORK( p ) / WORK( q )
1254                                     AQOAP = WORK( q ) / WORK( p )
1255                                     IF( WORK( p ).GE.ONE ) THEN
1256 *
1257                                        IF( WORK( q ).GE.ONE ) THEN
1258                                           FASTR( 3 ) = T*APOAQ
1259                                           FASTR( 4 ) = -T*AQOAP
1260                                           WORK( p ) = WORK( p )*CS
1261                                           WORK( q ) = WORK( q )*CS
1262                                           CALL SROTM( M, A( 1, p ), 1,
1263      $                                                A( 1, q ), 1,
1264      $                                                FASTR )
1265                                           IF( RSVEC )CALL SROTM( MVL,
1266      $                                        V( 1, p ), 1, V( 1, q ),
1267      $                                        1, FASTR )
1268                                        ELSE
1269                                           CALL SAXPY( M, -T*AQOAP,
1270      $                                                A( 1, q ), 1,
1271      $                                                A( 1, p ), 1 )
1272                                           CALL SAXPY( M, CS*SN*APOAQ,
1273      $                                                A( 1, p ), 1,
1274      $                                                A( 1, q ), 1 )
1275                                           IF( RSVEC ) THEN
1276                                              CALL SAXPY( MVL, -T*AQOAP,
1277      $                                                   V( 1, q ), 1,
1278      $                                                   V( 1, p ), 1 )
1279                                              CALL SAXPY( MVL,
1280      $                                                   CS*SN*APOAQ,
1281      $                                                   V( 1, p ), 1,
1282      $                                                   V( 1, q ), 1 )
1283                                           END IF
1284                                           WORK( p ) = WORK( p )*CS
1285                                           WORK( q ) = WORK( q ) / CS
1286                                        END IF
1287                                     ELSE
1288                                        IF( WORK( q ).GE.ONE ) THEN
1289                                           CALL SAXPY( M, T*APOAQ,
1290      $                                                A( 1, p ), 1,
1291      $                                                A( 1, q ), 1 )
1292                                           CALL SAXPY( M, -CS*SN*AQOAP,
1293      $                                                A( 1, q ), 1,
1294      $                                                A( 1, p ), 1 )
1295                                           IF( RSVEC ) THEN
1296                                              CALL SAXPY( MVL, T*APOAQ,
1297      $                                                   V( 1, p ), 1,
1298      $                                                   V( 1, q ), 1 )
1299                                              CALL SAXPY( MVL,
1300      $                                                   -CS*SN*AQOAP,
1301      $                                                   V( 1, q ), 1,
1302      $                                                   V( 1, p ), 1 )
1303                                           END IF
1304                                           WORK( p ) = WORK( p ) / CS
1305                                           WORK( q ) = WORK( q )*CS
1306                                        ELSE
1307                                           IF( WORK( p ).GE.WORK( q ) )
1308      $                                        THEN
1309                                              CALL SAXPY( M, -T*AQOAP,
1310      $                                                   A( 1, q ), 1,
1311      $                                                   A( 1, p ), 1 )
1312                                              CALL SAXPY( M, CS*SN*APOAQ,
1313      $                                                   A( 1, p ), 1,
1314      $                                                   A( 1, q ), 1 )
1315                                              WORK( p ) = WORK( p )*CS
1316                                              WORK( q ) = WORK( q ) / CS
1317                                              IF( RSVEC ) THEN
1318                                                 CALL SAXPY( MVL,
1319      $                                               -T*AQOAP,
1320      $                                               V( 1, q ), 1,
1321      $                                               V( 1, p ), 1 )
1322                                                 CALL SAXPY( MVL,
1323      $                                               CS*SN*APOAQ,
1324      $                                               V( 1, p ), 1,
1325      $                                               V( 1, q ), 1 )
1326                                              END IF
1327                                           ELSE
1328                                              CALL SAXPY( M, T*APOAQ,
1329      $                                                   A( 1, p ), 1,
1330      $                                                   A( 1, q ), 1 )
1331                                              CALL SAXPY( M,
1332      $                                                   -CS*SN*AQOAP,
1333      $                                                   A( 1, q ), 1,
1334      $                                                   A( 1, p ), 1 )
1335                                              WORK( p ) = WORK( p ) / CS
1336                                              WORK( q ) = WORK( q )*CS
1337                                              IF( RSVEC ) THEN
1338                                                 CALL SAXPY( MVL,
1339      $                                               T*APOAQ, V( 1, p ),
1340      $                                               1, V( 1, q ), 1 )
1341                                                 CALL SAXPY( MVL,
1342      $                                               -CS*SN*AQOAP,
1343      $                                               V( 1, q ), 1,
1344      $                                               V( 1, p ), 1 )
1345                                              END IF
1346                                           END IF
1347                                        END IF
1348                                     END IF
1349                                  END IF
1350 *
1351                               ELSE
1352                                  IF( AAPP.GT.AAQQ ) THEN
1353                                     CALL SCOPY( M, A( 1, p ), 1,
1354      $                                          WORK( N+1 ), 1 )
1355                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
1356      $                                           M, 1, WORK( N+1 ), LDA,
1357      $                                           IERR )
1358                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
1359      $                                           M, 1, A( 1, q ), LDA,
1360      $                                           IERR )
1361                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )
1362                                     CALL SAXPY( M, TEMP1, WORK( N+1 ),
1363      $                                          1, A( 1, q ), 1 )
1364                                     CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
1365      $                                           M, 1, A( 1, q ), LDA,
1366      $                                           IERR )
1367                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1368      $                                         ONE-AAPQ*AAPQ ) )
1369                                     MXSINJ = MAX( MXSINJ, SFMIN )
1370                                  ELSE
1371                                     CALL SCOPY( M, A( 1, q ), 1,
1372      $                                          WORK( N+1 ), 1 )
1373                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
1374      $                                           M, 1, WORK( N+1 ), LDA,
1375      $                                           IERR )
1376                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
1377      $                                           M, 1, A( 1, p ), LDA,
1378      $                                           IERR )
1379                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )
1380                                     CALL SAXPY( M, TEMP1, WORK( N+1 ),
1381      $                                          1, A( 1, p ), 1 )
1382                                     CALL SLASCL( 'G', 0, 0, ONE, AAPP,
1383      $                                           M, 1, A( 1, p ), LDA,
1384      $                                           IERR )
1385                                     SVA( p ) = AAPP*SQRT( MAX( ZERO,
1386      $                                         ONE-AAPQ*AAPQ ) )
1387                                     MXSINJ = MAX( MXSINJ, SFMIN )
1388                                  END IF
1389                               END IF
1390 *           END IF ROTOK THEN ... ELSE
1391 *
1392 *           In the case of cancellation in updating SVA(q)
1393 *           .. recompute SVA(q)
1394                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1395      $                            THEN
1396                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
1397      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
1398                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
1399      $                                         WORK( q )
1400                                  ELSE
1401                                     T = ZERO
1402                                     AAQQ = ONE
1403                                     CALL SLASSQ( M, A( 1, q ), 1, T,
1404      $                                           AAQQ )
1405                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
1406                                  END IF
1407                               END IF
1408                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
1409                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
1410      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1411                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
1412      $                                     WORK( p )
1413                                  ELSE
1414                                     T = ZERO
1415                                     AAPP = ONE
1416                                     CALL SLASSQ( M, A( 1, p ), 1, T,
1417      $                                           AAPP )
1418                                     AAPP = T*SQRT( AAPP )*WORK( p )
1419                                  END IF
1420                                  SVA( p ) = AAPP
1421                               END IF
1422 *              end of OK rotation
1423                            ELSE
1424                               NOTROT = NOTROT + 1
1425 *[RTD]      SKIPPED  = SKIPPED  + 1
1426                               PSKIPPED = PSKIPPED + 1
1427                               IJBLSK = IJBLSK + 1
1428                            END IF
1429                         ELSE
1430                            NOTROT = NOTROT + 1
1431                            PSKIPPED = PSKIPPED + 1
1432                            IJBLSK = IJBLSK + 1
1433                         END IF
1434 *
1435                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
1436      $                      THEN
1437                            SVA( p ) = AAPP
1438                            NOTROT = 0
1439                            GO TO 2011
1440                         END IF
1441                         IF( ( i.LE.SWBAND ) .AND.
1442      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1443                            AAPP = -AAPP
1444                            NOTROT = 0
1445                            GO TO 2203
1446                         END IF
1447 *
1448  2200                CONTINUE
1449 *        end of the q-loop
1450  2203                CONTINUE
1451 *
1452                      SVA( p ) = AAPP
1453 *
1454                   ELSE
1455 *
1456                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1457      $                   MIN( jgl+KBL-1, N ) - jgl + 1
1458                      IF( AAPP.LT.ZERO )NOTROT = 0
1459 *
1460                   END IF
1461 *
1462  2100          CONTINUE
1463 *     end of the p-loop
1464  2010       CONTINUE
1465 *     end of the jbc-loop
1466  2011       CONTINUE
1467 *2011 bailed out of the jbc-loop
1468             DO 2012 p = igl, MIN( igl+KBL-1, N )
1469                SVA( p ) = ABS( SVA( p ) )
1470  2012       CONTINUE
1471 ***
1472  2000    CONTINUE
1473 *2000 :: end of the ibr-loop
1474 *
1475 *     .. update SVA(N)
1476          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
1477      $       THEN
1478             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*WORK( N )
1479          ELSE
1480             T = ZERO
1481             AAPP = ONE
1482             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
1483             SVA( N ) = T*SQRT( AAPP )*WORK( N )
1484          END IF
1485 *
1486 *     Additional steering devices
1487 *
1488          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
1489      $       ( ISWROT.LE.N ) ) )SWBAND = i
1490 *
1491          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
1492      $       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1493             GO TO 1994
1494          END IF
1495 *
1496          IF( NOTROT.GE.EMPTSW )GO TO 1994
1497 *
1498  1993 CONTINUE
1499 *     end i=1:NSWEEP loop
1500 *
1501 * #:( Reaching this point means that the procedure has not converged.
1502       INFO = NSWEEP - 1
1503       GO TO 1995
1504 *
1505  1994 CONTINUE
1506 * #:) Reaching this point means numerical convergence after the i-th
1507 *     sweep.
1508 *
1509       INFO = 0
1510 * #:) INFO = 0 confirms successful iterations.
1511  1995 CONTINUE
1512 *
1513 *     Sort the singular values and find how many are above
1514 *     the underflow threshold.
1515 *
1516       N2 = 0
1517       N4 = 0
1518       DO 5991 p = 1, N - 1
1519          q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
1520          IF( p.NE.q ) THEN
1521             TEMP1 = SVA( p )
1522             SVA( p ) = SVA( q )
1523             SVA( q ) = TEMP1
1524             TEMP1 = WORK( p )
1525             WORK( p ) = WORK( q )
1526             WORK( q ) = TEMP1
1527             CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
1528             IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
1529          END IF
1530          IF( SVA( p ).NE.ZERO ) THEN
1531             N4 = N4 + 1
1532             IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
1533          END IF
1534  5991 CONTINUE
1535       IF( SVA( N ).NE.ZERO ) THEN
1536          N4 = N4 + 1
1537          IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
1538       END IF
1539 *
1540 *     Normalize the left singular vectors.
1541 *
1542       IF( LSVEC .OR. UCTOL ) THEN
1543          DO 1998 p = 1, N2
1544             CALL SSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
1545  1998    CONTINUE
1546       END IF
1547 *
1548 *     Scale the product of Jacobi rotations (assemble the fast rotations).
1549 *
1550       IF( RSVEC ) THEN
1551          IF( APPLV ) THEN
1552             DO 2398 p = 1, N
1553                CALL SSCAL( MVL, WORK( p ), V( 1, p ), 1 )
1554  2398       CONTINUE
1555          ELSE
1556             DO 2399 p = 1, N
1557                TEMP1 = ONE / SNRM2( MVL, V( 1, p ), 1 )
1558                CALL SSCAL( MVL, TEMP1, V( 1, p ), 1 )
1559  2399       CONTINUE
1560          END IF
1561       END IF
1562 *
1563 *     Undo scaling, if necessary (and possible).
1564       IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL ) ) )
1565      $    .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT.
1566      $    ( SFMIN / SKL ) ) ) ) THEN
1567          DO 2400 p = 1, N
1568             SVA( P ) = SKL*SVA( P )
1569  2400    CONTINUE
1570          SKL = ONE
1571       END IF
1572 *
1573       WORK( 1 ) = SKL
1574 *     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
1575 *     then some of the singular values may overflow or underflow and
1576 *     the spectrum is given in this factored representation.
1577 *
1578       WORK( 2 ) = FLOAT( N4 )
1579 *     N4 is the number of computed nonzero singular values of A.
1580 *
1581       WORK( 3 ) = FLOAT( N2 )
1582 *     N2 is the number of singular values of A greater than SFMIN.
1583 *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
1584 *     that may carry some information.
1585 *
1586       WORK( 4 ) = FLOAT( i )
1587 *     i is the index of the last sweep before declaring convergence.
1588 *
1589       WORK( 5 ) = MXAAPQ
1590 *     MXAAPQ is the largest absolute value of scaled pivots in the
1591 *     last sweep
1592 *
1593       WORK( 6 ) = MXSINJ
1594 *     MXSINJ is the largest absolute value of the sines of Jacobi angles
1595 *     in the last sweep
1596 *
1597       RETURN
1598 *     ..
1599 *     .. END OF SGESVJ
1600 *     ..
1601       END