SRC/cgesvj.f

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