SRC/cgsvj0.f

   1 *> \brief \b CGSVJ0 pre-processor for the routine cgesvj.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGSVJ0 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgsvj0.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgsvj0.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgsvj0.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
  22 *                          SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
  26 *       REAL               EPS, SFMIN, TOL
  27 *       CHARACTER*1        JOBV
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       COMPLEX            A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
  31 *       REAL               SVA( N )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> CGSVJ0 is called from CGESVJ as a pre-processor and that is its main
  41 *> purpose. It applies Jacobi rotations in the same way as CGESVJ does, but
  42 *> it does not check convergence (stopping criterion). Few tuning
  43 *> parameters (marked by [TP]) are available for the implementer.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] JOBV
  50 *> \verbatim
  51 *>          JOBV is CHARACTER*1
  52 *>          Specifies whether the output from this procedure is used
  53 *>          to compute the matrix V:
  54 *>          = 'V': the product of the Jacobi rotations is accumulated
  55 *>                 by postmulyiplying the N-by-N array V.
  56 *>                (See the description of V.)
  57 *>          = 'A': the product of the Jacobi rotations is accumulated
  58 *>                 by postmulyiplying the MV-by-N array V.
  59 *>                (See the descriptions of MV and V.)
  60 *>          = 'N': the Jacobi rotations are not accumulated.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] M
  64 *> \verbatim
  65 *>          M is INTEGER
  66 *>          The number of rows of the input matrix A.  M >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] N
  70 *> \verbatim
  71 *>          N is INTEGER
  72 *>          The number of columns of the input matrix A.
  73 *>          M >= N >= 0.
  74 *> \endverbatim
  75 *>
  76 *> \param[in,out] A
  77 *> \verbatim
  78 *>          A is COMPLEX array, dimension (LDA,N)
  79 *>          On entry, M-by-N matrix A, such that A*diag(D) represents
  80 *>          the input matrix.
  81 *>          On exit,
  82 *>          A_onexit * diag(D_onexit) represents the input matrix A*diag(D)
  83 *>          post-multiplied by a sequence of Jacobi rotations, where the
  84 *>          rotation threshold and the total number of sweeps are given in
  85 *>          TOL and NSWEEP, respectively.
  86 *>          (See the descriptions of D, TOL and NSWEEP.)
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDA
  90 *> \verbatim
  91 *>          LDA is INTEGER
  92 *>          The leading dimension of the array A.  LDA >= max(1,M).
  93 *> \endverbatim
  94 *>
  95 *> \param[in,out] D
  96 *> \verbatim
  97 *>          D is COMPLEX array, dimension (N)
  98 *>          The array D accumulates the scaling factors from the complex scaled
  99 *>          Jacobi rotations.
 100 *>          On entry, A*diag(D) represents the input matrix.
 101 *>          On exit, A_onexit*diag(D_onexit) represents the input matrix
 102 *>          post-multiplied by a sequence of Jacobi rotations, where the
 103 *>          rotation threshold and the total number of sweeps are given in
 104 *>          TOL and NSWEEP, respectively.
 105 *>          (See the descriptions of A, TOL and NSWEEP.)
 106 *> \endverbatim
 107 *>
 108 *> \param[in,out] SVA
 109 *> \verbatim
 110 *>          SVA is REAL array, dimension (N)
 111 *>          On entry, SVA contains the Euclidean norms of the columns of
 112 *>          the matrix A*diag(D).
 113 *>          On exit, SVA contains the Euclidean norms of the columns of
 114 *>          the matrix A_onexit*diag(D_onexit).
 115 *> \endverbatim
 116 *>
 117 *> \param[in] MV
 118 *> \verbatim
 119 *>          MV is INTEGER
 120 *>          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
 121 *>                           sequence of Jacobi rotations.
 122 *>          If JOBV = 'N',   then MV is not referenced.
 123 *> \endverbatim
 124 *>
 125 *> \param[in,out] V
 126 *> \verbatim
 127 *>          V is COMPLEX array, dimension (LDV,N)
 128 *>          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
 129 *>                           sequence of Jacobi rotations.
 130 *>          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
 131 *>                           sequence of Jacobi rotations.
 132 *>          If JOBV = 'N',   then V is not referenced.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] LDV
 136 *> \verbatim
 137 *>          LDV is INTEGER
 138 *>          The leading dimension of the array V,  LDV >= 1.
 139 *>          If JOBV = 'V', LDV .GE. N.
 140 *>          If JOBV = 'A', LDV .GE. MV.
 141 *> \endverbatim
 142 *>
 143 *> \param[in] EPS
 144 *> \verbatim
 145 *>          EPS is REAL
 146 *>          EPS = SLAMCH('Epsilon')
 147 *> \endverbatim
 148 *>
 149 *> \param[in] SFMIN
 150 *> \verbatim
 151 *>          SFMIN is REAL
 152 *>          SFMIN = SLAMCH('Safe Minimum')
 153 *> \endverbatim
 154 *>
 155 *> \param[in] TOL
 156 *> \verbatim
 157 *>          TOL is REAL
 158 *>          TOL is the threshold for Jacobi rotations. For a pair
 159 *>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
 160 *>          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
 161 *> \endverbatim
 162 *>
 163 *> \param[in] NSWEEP
 164 *> \verbatim
 165 *>          NSWEEP is INTEGER
 166 *>          NSWEEP is the number of sweeps of Jacobi rotations to be
 167 *>          performed.
 168 *> \endverbatim
 169 *>
 170 *> \param[out] WORK
 171 *> \verbatim
 172 *>          WORK is COMPLEX array, dimension LWORK.
 173 *> \endverbatim
 174 *>
 175 *> \param[in] LWORK
 176 *> \verbatim
 177 *>          LWORK is INTEGER
 178 *>          LWORK is the dimension of WORK. LWORK .GE. M.
 179 *> \endverbatim
 180 *>
 181 *> \param[out] INFO
 182 *> \verbatim
 183 *>          INFO is INTEGER
 184 *>          = 0 : successful exit.
 185 *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
 186 *> \endverbatim
 187 *
 188 *  Authors:
 189 *  ========
 190 *
 191 *> \author Univ. of Tennessee
 192 *> \author Univ. of California Berkeley
 193 *> \author Univ. of Colorado Denver
 194 *> \author NAG Ltd.
 195 *
 196 *> \date June 2016
 197 *
 198 *> \ingroup complexOTHERcomputational
 199 *
 200 *> \par Further Details:
 201 *  =====================
 202 *>
 203 *> CGSVJ0 is used just to enable CGESVJ to call a simplified version of
 204 *> itself to work on a submatrix of the original matrix.
 205 *>
 206 *> \par Contributors:
 207 *  ==================
 208 *>
 209 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
 210 *>
 211 *> \par Bugs, Examples and Comments:
 212 *  =================================
 213 *>
 214 *> Please report all bugs and send interesting test examples and comments to
 215 *> drmac@math.hr. Thank you.
 216 *
 217 *  =====================================================================
 218       SUBROUTINE CGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
 219      $                   SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 220 *
 221 *  -- LAPACK computational routine (version 3.6.1) --
 222 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 223 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 224 *     June 2016
 225 *
 226       IMPLICIT NONE
 227 *     .. Scalar Arguments ..
 228       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
 229       REAL               EPS, SFMIN, TOL
 230       CHARACTER*1        JOBV
 231 *     ..
 232 *     .. Array Arguments ..
 233       COMPLEX            A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
 234       REAL               SVA( N )
 235 *     ..
 236 *
 237 *  =====================================================================
 238 *
 239 *     .. Local Parameters ..
 240       REAL               ZERO, HALF, ONE
 241       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0)
 242       COMPLEX      CZERO,                  CONE
 243       PARAMETER  ( CZERO = (0.0E0, 0.0E0), CONE = (1.0E0, 0.0E0) )
 244 *     ..
 245 *     .. Local Scalars ..
 246       COMPLEX            AAPQ, OMPQ
 247       REAL               AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
 248      $                   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
 249      $                   ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
 250      $                   THSIGN
 251       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
 252      $                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
 253      $                   NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
 254       LOGICAL            APPLV, ROTOK, RSVEC
 255 *     ..
 256 *     ..
 257 *     .. Intrinsic Functions ..
 258       INTRINSIC ABS, AMAX1, CONJG, FLOAT, MIN0, SIGN, SQRT
 259 *     ..
 260 *     .. External Functions ..
 261       REAL               SCNRM2
 262       COMPLEX            CDOTC
 263       INTEGER            ISAMAX
 264       LOGICAL            LSAME
 265       EXTERNAL           ISAMAX, LSAME, CDOTC, SCNRM2
 266 *     ..
 267 *     ..
 268 *     .. External Subroutines ..
 269 *     ..
 270 *     from BLAS
 271       EXTERNAL           CCOPY, CROT, CSSCAL, CSWAP
 272 *     from LAPACK
 273       EXTERNAL           CLASCL, CLASSQ, XERBLA
 274 *     ..
 275 *     .. Executable Statements ..
 276 *
 277 *     Test the input parameters.
 278 *
 279       APPLV = LSAME( JOBV, 'A' )
 280       RSVEC = LSAME( JOBV, 'V' )
 281       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 282          INFO = -1
 283       ELSE IF( M.LT.0 ) THEN
 284          INFO = -2
 285       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
 286          INFO = -3
 287       ELSE IF( LDA.LT.M ) THEN
 288          INFO = -5
 289       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
 290          INFO = -8
 291       ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
 292      $         ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
 293          INFO = -10
 294       ELSE IF( TOL.LE.EPS ) THEN
 295          INFO = -13
 296       ELSE IF( NSWEEP.LT.0 ) THEN
 297          INFO = -14
 298       ELSE IF( LWORK.LT.M ) THEN
 299          INFO = -16
 300       ELSE
 301          INFO = 0
 302       END IF
 303 *
 304 *     #:(
 305       IF( INFO.NE.0 ) THEN
 306          CALL XERBLA( 'CGSVJ0', -INFO )
 307          RETURN
 308       END IF
 309 *
 310       IF( RSVEC ) THEN
 311          MVL = N
 312       ELSE IF( APPLV ) THEN
 313          MVL = MV
 314       END IF
 315       RSVEC = RSVEC .OR. APPLV
 316
 317       ROOTEPS = SQRT( EPS )
 318       ROOTSFMIN = SQRT( SFMIN )
 319       SMALL = SFMIN / EPS
 320       BIG = ONE / SFMIN
 321       ROOTBIG = ONE / ROOTSFMIN
 322       BIGTHETA = ONE / ROOTEPS
 323       ROOTTOL = SQRT( TOL )
 324 *
 325 *     .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
 326 *
 327       EMPTSW = ( N*( N-1 ) ) / 2
 328       NOTROT = 0
 329 *
 330 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 331 *
 332
 333       SWBAND = 0
 334 *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
 335 *     if CGESVJ is used as a computational routine in the preconditioned
 336 *     Jacobi SVD algorithm CGEJSV. For sweeps i=1:SWBAND the procedure
 337 *     works on pivots inside a band-like region around the diagonal.
 338 *     The boundaries are determined dynamically, based on the number of
 339 *     pivots above a threshold.
 340 *
 341       KBL = MIN0( 8, N )
 342 *[TP] KBL is a tuning parameter that defines the tile size in the
 343 *     tiling of the p-q loops of pivot pairs. In general, an optimal
 344 *     value of KBL depends on the matrix dimensions and on the
 345 *     parameters of the computer's memory.
 346 *
 347       NBL = N / KBL
 348       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
 349 *
 350       BLSKIP = KBL**2
 351 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 352 *
 353       ROWSKIP = MIN0( 5, KBL )
 354 *[TP] ROWSKIP is a tuning parameter.
 355 *
 356       LKAHEAD = 1
 357 *[TP] LKAHEAD is a tuning parameter.
 358 *
 359 *     Quasi block transformations, using the lower (upper) triangular
 360 *     structure of the input matrix. The quasi-block-cycling usually
 361 *     invokes cubic convergence. Big part of this cycle is done inside
 362 *     canonical subspaces of dimensions less than M.
 363 *
 364 *
 365 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 366 *
 367       DO 1993 i = 1, NSWEEP
 368 *
 369 *     .. go go go ...
 370 *
 371          MXAAPQ = ZERO
 372          MXSINJ = ZERO
 373          ISWROT = 0
 374 *
 375          NOTROT = 0
 376          PSKIPPED = 0
 377 *
 378 *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
 379 *     1 <= p < q <= N. This is the first step toward a blocked implementation
 380 *     of the rotations. New implementation, based on block transformations,
 381 *     is under development.
 382 *
 383          DO 2000 ibr = 1, NBL
 384 *
 385             igl = ( ibr-1 )*KBL + 1
 386 *
 387             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
 388 *
 389                igl = igl + ir1*KBL
 390 *
 391                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
 392 *
 393 *     .. de Rijk's pivoting
 394 *
 395                   q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
 396                   IF( p.NE.q ) THEN
 397                      CALL CSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
 398                      IF( RSVEC )CALL CSWAP( MVL, V( 1, p ), 1,
 399      $                                           V( 1, q ), 1 )
 400                      TEMP1 = SVA( p )
 401                      SVA( p ) = SVA( q )
 402                      SVA( q ) = TEMP1
 403                      AAPQ = D(p)
 404                      D(p) = D(q)
 405                      D(q) = AAPQ
 406                   END IF
 407 *
 408                   IF( ir1.EQ.0 ) THEN
 409 *
 410 *        Column norms are periodically updated by explicit
 411 *        norm computation.
 412 *        Caveat:
 413 *        Unfortunately, some BLAS implementations compute SNCRM2(M,A(1,p),1)
 414 *        as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
 415 *        overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
 416 *        underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
 417 *        Hence, SCNRM2 cannot be trusted, not even in the case when
 418 *        the true norm is far from the under(over)flow boundaries.
 419 *        If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF
 420 *        below should be replaced with "AAPP = SCNRM2( M, A(1,p), 1 )".
 421 *
 422                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
 423      $                    ( SVA( p ).GT.ROOTSFMIN ) ) THEN
 424                         SVA( p ) = SCNRM2( M, A( 1, p ), 1 )
 425                      ELSE
 426                         TEMP1 = ZERO
 427                         AAPP = ONE
 428                         CALL CLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
 429                         SVA( p ) = TEMP1*SQRT( AAPP )
 430                      END IF
 431                      AAPP = SVA( p )
 432                   ELSE
 433                      AAPP = SVA( p )
 434                   END IF
 435 *
 436                   IF( AAPP.GT.ZERO ) THEN
 437 *
 438                      PSKIPPED = 0
 439 *
 440                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
 441 *
 442                         AAQQ = SVA( q )
 443 *
 444                         IF( AAQQ.GT.ZERO ) THEN
 445 *
 446                            AAPP0 = AAPP
 447                            IF( AAQQ.GE.ONE ) THEN
 448                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
 449                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
 450                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
 451      $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
 452                               ELSE
 453                                  CALL CCOPY( M, A( 1, p ), 1,
 454      $                                        WORK, 1 )
 455                                  CALL CLASCL( 'G', 0, 0, AAPP, ONE,
 456      $                                M, 1, WORK, LDA, IERR )
 457                                  AAPQ = CDOTC( M, WORK, 1,
 458      $                                   A( 1, q ), 1 ) / AAQQ
 459                               END IF
 460                            ELSE
 461                               ROTOK = AAPP.LE.( AAQQ / SMALL )
 462                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
 463                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
 464      $                                    A( 1, q ), 1 ) / AAQQ ) / AAPP
 465                               ELSE
 466                                  CALL CCOPY( M, A( 1, q ), 1,
 467      $                                        WORK, 1 )
 468                                  CALL CLASCL( 'G', 0, 0, AAQQ,
 469      $                                         ONE, M, 1,
 470      $                                         WORK, LDA, IERR )
 471                                  AAPQ = CDOTC( M, A( 1, p ), 1,
 472      $                                   WORK, 1 ) / AAPP
 473                               END IF
 474                            END IF
 475 *
 476                            OMPQ = AAPQ / ABS(AAPQ)
 477 *                           AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
 478                            AAPQ1  = -ABS(AAPQ)
 479                            MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
 480 *
 481 *        TO rotate or NOT to rotate, THAT is the question ...
 482 *
 483                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
 484 *
 485 *           .. rotate
 486 *[RTD]      ROTATED = ROTATED + ONE
 487 *
 488                               IF( ir1.EQ.0 ) THEN
 489                                  NOTROT = 0
 490                                  PSKIPPED = 0
 491                                  ISWROT = ISWROT + 1
 492                               END IF
 493 *
 494                               IF( ROTOK ) THEN
 495 *
 496                                  AQOAP = AAQQ / AAPP
 497                                  APOAQ = AAPP / AAQQ
 498                                  THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1
 499 *
 500                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
 501 *
 502                                     T  = HALF / THETA
 503                                     CS = ONE
 504
 505                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
 506      $                                          CS, CONJG(OMPQ)*T )
 507                                     IF ( RSVEC ) THEN
 508                                         CALL CROT( MVL, V(1,p), 1,
 509      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
 510                                     END IF
 511
 512                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
 513      $                                          ONE+T*APOAQ*AAPQ1 ) )
 514                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
 515      $                                          ONE-T*AQOAP*AAPQ1 ) )
 516                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
 517 *
 518                                  ELSE
 519 *
 520 *                 .. choose correct signum for THETA and rotate
 521 *
 522                                     THSIGN = -SIGN( ONE, AAPQ1 )
 523                                     T = ONE / ( THETA+THSIGN*
 524      $                                   SQRT( ONE+THETA*THETA ) )
 525                                     CS = SQRT( ONE / ( ONE+T*T ) )
 526                                     SN = T*CS
 527 *
 528                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
 529                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
 530      $                                          ONE+T*APOAQ*AAPQ1 ) )
 531                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
 532      $                                      ONE-T*AQOAP*AAPQ1 ) )
 533 *
 534                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
 535      $                                          CS, CONJG(OMPQ)*SN )
 536                                     IF ( RSVEC ) THEN
 537                                         CALL CROT( MVL, V(1,p), 1,
 538      $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
 539                                     END IF
 540                                  END IF
 541                                  D(p) = -D(q) * OMPQ
 542 *
 543                                  ELSE
 544 *              .. have to use modified Gram-Schmidt like transformation
 545                                  CALL CCOPY( M, A( 1, p ), 1,
 546      $                                       WORK, 1 )
 547                                  CALL CLASCL( 'G', 0, 0, AAPP, ONE, M,
 548      $                                        1, WORK, LDA,
 549      $                                        IERR )
 550                                  CALL CLASCL( 'G', 0, 0, AAQQ, ONE, M,
 551      $                                        1, A( 1, q ), LDA, IERR )
 552                                  CALL CAXPY( M, -AAPQ, WORK, 1,
 553      $                                       A( 1, q ), 1 )
 554                                  CALL CLASCL( 'G', 0, 0, ONE, AAQQ, M,
 555      $                                        1, A( 1, q ), LDA, IERR )
 556                                  SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
 557      $                                      ONE-AAPQ1*AAPQ1 ) )
 558                                  MXSINJ = AMAX1( MXSINJ, SFMIN )
 559                               END IF
 560 *           END IF ROTOK THEN ... ELSE
 561 *
 562 *           In the case of cancellation in updating SVA(q), SVA(p)
 563 *           recompute SVA(q), SVA(p).
 564 *
 565                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
 566      $                            THEN
 567                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
 568      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
 569                                     SVA( q ) = SCNRM2( M, A( 1, q ), 1 )
 570                                  ELSE
 571                                     T = ZERO
 572                                     AAQQ = ONE
 573                                     CALL CLASSQ( M, A( 1, q ), 1, T,
 574      $                                           AAQQ )
 575                                     SVA( q ) = T*SQRT( AAQQ )
 576                                  END IF
 577                               END IF
 578                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
 579                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
 580      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
 581                                     AAPP = SCNRM2( M, A( 1, p ), 1 )
 582                                  ELSE
 583                                     T = ZERO
 584                                     AAPP = ONE
 585                                     CALL CLASSQ( M, A( 1, p ), 1, T,
 586      $                                           AAPP )
 587                                     AAPP = T*SQRT( AAPP )
 588                                  END IF
 589                                  SVA( p ) = AAPP
 590                               END IF
 591 *
 592                            ELSE
 593 *        A(:,p) and A(:,q) already numerically orthogonal
 594                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
 595 *[RTD]      SKIPPED  = SKIPPED  + 1
 596                               PSKIPPED = PSKIPPED + 1
 597                            END IF
 598                         ELSE
 599 *        A(:,q) is zero column
 600                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
 601                            PSKIPPED = PSKIPPED + 1
 602                         END IF
 603 *
 604                         IF( ( i.LE.SWBAND ) .AND.
 605      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
 606                            IF( ir1.EQ.0 )AAPP = -AAPP
 607                            NOTROT = 0
 608                            GO TO 2103
 609                         END IF
 610 *
 611  2002                CONTINUE
 612 *     END q-LOOP
 613 *
 614  2103                CONTINUE
 615 *     bailed out of q-loop
 616 *
 617                      SVA( p ) = AAPP
 618 *
 619                   ELSE
 620                      SVA( p ) = AAPP
 621                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
 622      $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
 623                   END IF
 624 *
 625  2001          CONTINUE
 626 *     end of the p-loop
 627 *     end of doing the block ( ibr, ibr )
 628  1002       CONTINUE
 629 *     end of ir1-loop
 630 *
 631 * ... go to the off diagonal blocks
 632 *
 633             igl = ( ibr-1 )*KBL + 1
 634 *
 635             DO 2010 jbc = ibr + 1, NBL
 636 *
 637                jgl = ( jbc-1 )*KBL + 1
 638 *
 639 *        doing the block at ( ibr, jbc )
 640 *
 641                IJBLSK = 0
 642                DO 2100 p = igl, MIN0( igl+KBL-1, N )
 643 *
 644                   AAPP = SVA( p )
 645                   IF( AAPP.GT.ZERO ) THEN
 646 *
 647                      PSKIPPED = 0
 648 *
 649                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
 650 *
 651                         AAQQ = SVA( q )
 652                         IF( AAQQ.GT.ZERO ) THEN
 653                            AAPP0 = AAPP
 654 *
 655 *     .. M x 2 Jacobi SVD ..
 656 *
 657 *        Safe Gram matrix computation
 658 *
 659                            IF( AAQQ.GE.ONE ) THEN
 660                               IF( AAPP.GE.AAQQ ) THEN
 661                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
 662                               ELSE
 663                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
 664                               END IF
 665                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
 666                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
 667      $                                  A( 1, q ), 1 ) / AAQQ ) / AAPP
 668                               ELSE
 669                                  CALL CCOPY( M, A( 1, p ), 1,
 670      $                                       WORK, 1 )
 671                                  CALL CLASCL( 'G', 0, 0, AAPP,
 672      $                                        ONE, M, 1,
 673      $                                        WORK, LDA, IERR )
 674                                  AAPQ = CDOTC( M, WORK, 1,
 675      $                                  A( 1, q ), 1 ) / AAQQ
 676                               END IF
 677                            ELSE
 678                               IF( AAPP.GE.AAQQ ) THEN
 679                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
 680                               ELSE
 681                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
 682                               END IF
 683                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
 684                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
 685      $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
 686                               ELSE
 687                                  CALL CCOPY( M, A( 1, q ), 1,
 688      $                                       WORK, 1 )
 689                                  CALL CLASCL( 'G', 0, 0, AAQQ,
 690      $                                        ONE, M, 1,
 691      $                                        WORK, LDA, IERR )
 692                                  AAPQ = CDOTC( M, A( 1, p ), 1,
 693      $                                  WORK, 1 ) / AAPP
 694                               END IF
 695                            END IF
 696 *
 697                            OMPQ = AAPQ / ABS(AAPQ)
 698 *                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
 699                            AAPQ1  = -ABS(AAPQ)
 700                            MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
 701 *
 702 *        TO rotate or NOT to rotate, THAT is the question ...
 703 *
 704                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
 705                               NOTROT = 0
 706 *[RTD]      ROTATED  = ROTATED + 1
 707                               PSKIPPED = 0
 708                               ISWROT = ISWROT + 1
 709 *
 710                               IF( ROTOK ) THEN
 711 *
 712                                  AQOAP = AAQQ / AAPP
 713                                  APOAQ = AAPP / AAQQ
 714                                  THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
 715                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
 716 *
 717                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
 718                                     T  = HALF / THETA
 719                                     CS = ONE
 720                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
 721      $                                          CS, CONJG(OMPQ)*T )
 722                                     IF( RSVEC ) THEN
 723                                         CALL CROT( MVL, V(1,p), 1,
 724      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
 725                                     END IF
 726                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
 727      $                                         ONE+T*APOAQ*AAPQ1 ) )
 728                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
 729      $                                     ONE-T*AQOAP*AAPQ1 ) )
 730                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
 731                                  ELSE
 732 *
 733 *                 .. choose correct signum for THETA and rotate
 734 *
 735                                     THSIGN = -SIGN( ONE, AAPQ1 )
 736                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
 737                                     T = ONE / ( THETA+THSIGN*
 738      $                                  SQRT( ONE+THETA*THETA ) )
 739                                     CS = SQRT( ONE / ( ONE+T*T ) )
 740                                     SN = T*CS
 741                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
 742                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
 743      $                                         ONE+T*APOAQ*AAPQ1 ) )
 744                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
 745      $                                         ONE-T*AQOAP*AAPQ1 ) )
 746 *
 747                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
 748      $                                          CS, CONJG(OMPQ)*SN )
 749                                     IF( RSVEC ) THEN
 750                                         CALL CROT( MVL, V(1,p), 1,
 751      $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
 752                                     END IF
 753                                  END IF
 754                                  D(p) = -D(q) * OMPQ
 755 *
 756                               ELSE
 757 *              .. have to use modified Gram-Schmidt like transformation
 758                                IF( AAPP.GT.AAQQ ) THEN
 759                                     CALL CCOPY( M, A( 1, p ), 1,
 760      $                                          WORK, 1 )
 761                                     CALL CLASCL( 'G', 0, 0, AAPP, ONE,
 762      $                                           M, 1, WORK,LDA,
 763      $                                           IERR )
 764                                     CALL CLASCL( 'G', 0, 0, AAQQ, ONE,
 765      $                                           M, 1, A( 1, q ), LDA,
 766      $                                           IERR )
 767                                     CALL CAXPY( M, -AAPQ, WORK,
 768      $                                          1, A( 1, q ), 1 )
 769                                     CALL CLASCL( 'G', 0, 0, ONE, AAQQ,
 770      $                                           M, 1, A( 1, q ), LDA,
 771      $                                           IERR )
 772                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
 773      $                                         ONE-AAPQ1*AAPQ1 ) )
 774                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
 775                                ELSE
 776                                    CALL CCOPY( M, A( 1, q ), 1,
 777      $                                          WORK, 1 )
 778                                     CALL CLASCL( 'G', 0, 0, AAQQ, ONE,
 779      $                                           M, 1, WORK,LDA,
 780      $                                           IERR )
 781                                     CALL CLASCL( 'G', 0, 0, AAPP, ONE,
 782      $                                           M, 1, A( 1, p ), LDA,
 783      $                                           IERR )
 784                                     CALL CAXPY( M, -CONJG(AAPQ),
 785      $                                   WORK, 1, A( 1, p ), 1 )
 786                                     CALL CLASCL( 'G', 0, 0, ONE, AAPP,
 787      $                                           M, 1, A( 1, p ), LDA,
 788      $                                           IERR )
 789                                     SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
 790      $                                         ONE-AAPQ1*AAPQ1 ) )
 791                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
 792                                END IF
 793                               END IF
 794 *           END IF ROTOK THEN ... ELSE
 795 *
 796 *           In the case of cancellation in updating SVA(q), SVA(p)
 797 *           .. recompute SVA(q), SVA(p)
 798                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
 799      $                            THEN
 800                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
 801      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
 802                                     SVA( q ) = SCNRM2( M, A( 1, q ), 1)
 803                                   ELSE
 804                                     T = ZERO
 805                                     AAQQ = ONE
 806                                     CALL CLASSQ( M, A( 1, q ), 1, T,
 807      $                                           AAQQ )
 808                                     SVA( q ) = T*SQRT( AAQQ )
 809                                  END IF
 810                               END IF
 811                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
 812                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
 813      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
 814                                     AAPP = SCNRM2( M, A( 1, p ), 1 )
 815                                  ELSE
 816                                     T = ZERO
 817                                     AAPP = ONE
 818                                     CALL CLASSQ( M, A( 1, p ), 1, T,
 819      $                                           AAPP )
 820                                     AAPP = T*SQRT( AAPP )
 821                                  END IF
 822                                  SVA( p ) = AAPP
 823                               END IF
 824 *              end of OK rotation
 825                            ELSE
 826                               NOTROT = NOTROT + 1
 827 *[RTD]      SKIPPED  = SKIPPED  + 1
 828                               PSKIPPED = PSKIPPED + 1
 829                               IJBLSK = IJBLSK + 1
 830                            END IF
 831                         ELSE
 832                            NOTROT = NOTROT + 1
 833                            PSKIPPED = PSKIPPED + 1
 834                            IJBLSK = IJBLSK + 1
 835                         END IF
 836 *
 837                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
 838      $                      THEN
 839                            SVA( p ) = AAPP
 840                            NOTROT = 0
 841                            GO TO 2011
 842                         END IF
 843                         IF( ( i.LE.SWBAND ) .AND.
 844      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
 845                            AAPP = -AAPP
 846                            NOTROT = 0
 847                            GO TO 2203
 848                         END IF
 849 *
 850  2200                CONTINUE
 851 *        end of the q-loop
 852  2203                CONTINUE
 853 *
 854                      SVA( p ) = AAPP
 855 *
 856                   ELSE
 857 *
 858                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
 859      $                   MIN0( jgl+KBL-1, N ) - jgl + 1
 860                      IF( AAPP.LT.ZERO )NOTROT = 0
 861 *
 862                   END IF
 863 *
 864  2100          CONTINUE
 865 *     end of the p-loop
 866  2010       CONTINUE
 867 *     end of the jbc-loop
 868  2011       CONTINUE
 869 *2011 bailed out of the jbc-loop
 870             DO 2012 p = igl, MIN0( igl+KBL-1, N )
 871                SVA( p ) = ABS( SVA( p ) )
 872  2012       CONTINUE
 873 ***
 874  2000    CONTINUE
 875 *2000 :: end of the ibr-loop
 876 *
 877 *     .. update SVA(N)
 878          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
 879      $       THEN
 880             SVA( N ) = SCNRM2( M, A( 1, N ), 1 )
 881          ELSE
 882             T = ZERO
 883             AAPP = ONE
 884             CALL CLASSQ( M, A( 1, N ), 1, T, AAPP )
 885             SVA( N ) = T*SQRT( AAPP )
 886          END IF
 887 *
 888 *     Additional steering devices
 889 *
 890          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
 891      $       ( ISWROT.LE.N ) ) )SWBAND = i
 892 *
 893          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
 894      $       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
 895             GO TO 1994
 896          END IF
 897 *
 898          IF( NOTROT.GE.EMPTSW )GO TO 1994
 899 *
 900  1993 CONTINUE
 901 *     end i=1:NSWEEP loop
 902 *
 903 * #:( Reaching this point means that the procedure has not converged.
 904       INFO = NSWEEP - 1
 905       GO TO 1995
 906 *
 907  1994 CONTINUE
 908 * #:) Reaching this point means numerical convergence after the i-th
 909 *     sweep.
 910 *
 911       INFO = 0
 912 * #:) INFO = 0 confirms successful iterations.
 913  1995 CONTINUE
 914 *
 915 *     Sort the vector SVA() of column norms.
 916       DO 5991 p = 1, N - 1
 917          q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
 918          IF( p.NE.q ) THEN
 919             TEMP1 = SVA( p )
 920             SVA( p ) = SVA( q )
 921             SVA( q ) = TEMP1
 922             AAPQ = D( p )
 923             D( p ) = D( q )
 924             D( q ) = AAPQ
 925             CALL CSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
 926             IF( RSVEC )CALL CSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
 927          END IF
 928  5991 CONTINUE
 929 *
 930       RETURN
 931 *     ..
 932 *     .. END OF CGSVJ0
 933 *     ..
 934       END