SRC/sgsvj0.f

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