SRC/sgsvj0.f

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