SRC/zgesvj.f

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