SRC/dgsvj1.f

   1 *> \brief \b DGSVJ1 pre-processor for the routine dgesvj, applies Jacobi rotations targeting only particular pivots.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGSVJ1 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj1.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj1.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj1.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
  22 *                          EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       DOUBLE PRECISION   EPS, SFMIN, TOL
  26 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
  27 *       CHARACTER*1        JOBV
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       DOUBLE PRECISION   A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
  31 *      $                   WORK( LWORK )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DGSVJ1 is called from DGESVJ as a pre-processor and that is its main
  41 *> purpose. It applies Jacobi rotations in the same way as DGESVJ does, but
  42 *> it targets only particular pivots and it does not check convergence
  43 *> (stopping criterion). Few tunning parameters (marked by [TP]) are
  44 *> available for the implementer.
  45 *>
  46 *> Further Details
  47 *> ~~~~~~~~~~~~~~~
  48 *> DGSVJ1 applies few sweeps of Jacobi rotations in the column space of
  49 *> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
  50 *> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
  51 *> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
  52 *> [x]'s in the following scheme:
  53 *>
  54 *>    | *  *  * [x] [x] [x]|
  55 *>    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
  56 *>    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
  57 *>    |[x] [x] [x] *  *  * |
  58 *>    |[x] [x] [x] *  *  * |
  59 *>    |[x] [x] [x] *  *  * |
  60 *>
  61 *> In terms of the columns of A, the first N1 columns are rotated 'against'
  62 *> the remaining N-N1 columns, trying to increase the angle between the
  63 *> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
  64 *> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
  65 *> The number of sweeps is given in NSWEEP and the orthogonality threshold
  66 *> is given in TOL.
  67 *> \endverbatim
  68 *
  69 *  Arguments:
  70 *  ==========
  71 *
  72 *> \param[in] JOBV
  73 *> \verbatim
  74 *>          JOBV is CHARACTER*1
  75 *>          Specifies whether the output from this procedure is used
  76 *>          to compute the matrix V:
  77 *>          = 'V': the product of the Jacobi rotations is accumulated
  78 *>                 by postmulyiplying the N-by-N array V.
  79 *>                (See the description of V.)
  80 *>          = 'A': the product of the Jacobi rotations is accumulated
  81 *>                 by postmulyiplying the MV-by-N array V.
  82 *>                (See the descriptions of MV and V.)
  83 *>          = 'N': the Jacobi rotations are not accumulated.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] M
  87 *> \verbatim
  88 *>          M is INTEGER
  89 *>          The number of rows of the input matrix A.  M >= 0.
  90 *> \endverbatim
  91 *>
  92 *> \param[in] N
  93 *> \verbatim
  94 *>          N is INTEGER
  95 *>          The number of columns of the input matrix A.
  96 *>          M >= N >= 0.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] N1
 100 *> \verbatim
 101 *>          N1 is INTEGER
 102 *>          N1 specifies the 2 x 2 block partition, the first N1 columns are
 103 *>          rotated 'against' the remaining N-N1 columns of A.
 104 *> \endverbatim
 105 *>
 106 *> \param[in,out] A
 107 *> \verbatim
 108 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 109 *>          On entry, M-by-N matrix A, such that A*diag(D) represents
 110 *>          the input matrix.
 111 *>          On exit,
 112 *>          A_onexit * D_onexit represents the input matrix A*diag(D)
 113 *>          post-multiplied by a sequence of Jacobi rotations, where the
 114 *>          rotation threshold and the total number of sweeps are given in
 115 *>          TOL and NSWEEP, respectively.
 116 *>          (See the descriptions of N1, D, TOL and NSWEEP.)
 117 *> \endverbatim
 118 *>
 119 *> \param[in] LDA
 120 *> \verbatim
 121 *>          LDA is INTEGER
 122 *>          The leading dimension of the array A.  LDA >= max(1,M).
 123 *> \endverbatim
 124 *>
 125 *> \param[in,out] D
 126 *> \verbatim
 127 *>          D is DOUBLE PRECISION array, dimension (N)
 128 *>          The array D accumulates the scaling factors from the fast scaled
 129 *>          Jacobi rotations.
 130 *>          On entry, A*diag(D) represents the input matrix.
 131 *>          On exit, A_onexit*diag(D_onexit) represents the input matrix
 132 *>          post-multiplied by a sequence of Jacobi rotations, where the
 133 *>          rotation threshold and the total number of sweeps are given in
 134 *>          TOL and NSWEEP, respectively.
 135 *>          (See the descriptions of N1, A, TOL and NSWEEP.)
 136 *> \endverbatim
 137 *>
 138 *> \param[in,out] SVA
 139 *> \verbatim
 140 *>          SVA is DOUBLE PRECISION array, dimension (N)
 141 *>          On entry, SVA contains the Euclidean norms of the columns of
 142 *>          the matrix A*diag(D).
 143 *>          On exit, SVA contains the Euclidean norms of the columns of
 144 *>          the matrix onexit*diag(D_onexit).
 145 *> \endverbatim
 146 *>
 147 *> \param[in] MV
 148 *> \verbatim
 149 *>          MV is INTEGER
 150 *>          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
 151 *>                           sequence of Jacobi rotations.
 152 *>          If JOBV = 'N',   then MV is not referenced.
 153 *> \endverbatim
 154 *>
 155 *> \param[in,out] V
 156 *> \verbatim
 157 *>          V is DOUBLE PRECISION array, dimension (LDV,N)
 158 *>          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
 159 *>                           sequence of Jacobi rotations.
 160 *>          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
 161 *>                           sequence of Jacobi rotations.
 162 *>          If JOBV = 'N',   then V is not referenced.
 163 *> \endverbatim
 164 *>
 165 *> \param[in] LDV
 166 *> \verbatim
 167 *>          LDV is INTEGER
 168 *>          The leading dimension of the array V,  LDV >= 1.
 169 *>          If JOBV = 'V', LDV .GE. N.
 170 *>          If JOBV = 'A', LDV .GE. MV.
 171 *> \endverbatim
 172 *>
 173 *> \param[in] EPS
 174 *> \verbatim
 175 *>          EPS is DOUBLE PRECISION
 176 *>          EPS = DLAMCH('Epsilon')
 177 *> \endverbatim
 178 *>
 179 *> \param[in] SFMIN
 180 *> \verbatim
 181 *>          SFMIN is DOUBLE PRECISION
 182 *>          SFMIN = DLAMCH('Safe Minimum')
 183 *> \endverbatim
 184 *>
 185 *> \param[in] TOL
 186 *> \verbatim
 187 *>          TOL is DOUBLE PRECISION
 188 *>          TOL is the threshold for Jacobi rotations. For a pair
 189 *>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
 190 *>          applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
 191 *> \endverbatim
 192 *>
 193 *> \param[in] NSWEEP
 194 *> \verbatim
 195 *>          NSWEEP is INTEGER
 196 *>          NSWEEP is the number of sweeps of Jacobi rotations to be
 197 *>          performed.
 198 *> \endverbatim
 199 *>
 200 *> \param[out] WORK
 201 *> \verbatim
 202 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 203 *> \endverbatim
 204 *>
 205 *> \param[in] LWORK
 206 *> \verbatim
 207 *>          LWORK is INTEGER
 208 *>          LWORK is the dimension of WORK. LWORK .GE. M.
 209 *> \endverbatim
 210 *>
 211 *> \param[out] INFO
 212 *> \verbatim
 213 *>          INFO is INTEGER
 214 *>          = 0 : successful exit.
 215 *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
 216 *> \endverbatim
 217 *
 218 *  Authors:
 219 *  ========
 220 *
 221 *> \author Univ. of Tennessee
 222 *> \author Univ. of California Berkeley
 223 *> \author Univ. of Colorado Denver
 224 *> \author NAG Ltd.
 225 *
 226 *> \date June 2016
 227 *
 228 *> \ingroup doubleOTHERcomputational
 229 *
 230 *> \par Contributors:
 231 *  ==================
 232 *>
 233 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
 234 *
 235 *  =====================================================================
 236       SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
 237      $                   EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
 238 *
 239 *  -- LAPACK computational routine (version 3.6.1) --
 240 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 241 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 242 *     June 2016
 243 *
 244 *     .. Scalar Arguments ..
 245       DOUBLE PRECISION   EPS, SFMIN, TOL
 246       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
 247       CHARACTER*1        JOBV
 248 *     ..
 249 *     .. Array Arguments ..
 250       DOUBLE PRECISION   A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
 251      $                   WORK( LWORK )
 252 *     ..
 253 *
 254 *  =====================================================================
 255 *
 256 *     .. Local Parameters ..
 257       DOUBLE PRECISION   ZERO, HALF, ONE
 258       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 )
 259 *     ..
 260 *     .. Local Scalars ..
 261       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
 262      $                   BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
 263      $                   ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
 264      $                   TEMP1, THETA, THSIGN
 265       INTEGER            BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
 266      $                   ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
 267      $                   p, PSKIPPED, q, ROWSKIP, SWBAND
 268       LOGICAL            APPLV, ROTOK, RSVEC
 269 *     ..
 270 *     .. Local Arrays ..
 271       DOUBLE PRECISION   FASTR( 5 )
 272 *     ..
 273 *     .. Intrinsic Functions ..
 274       INTRINSIC          DABS, MAX, DBLE, MIN, DSIGN, DSQRT
 275 *     ..
 276 *     .. External Functions ..
 277       DOUBLE PRECISION   DDOT, DNRM2
 278       INTEGER            IDAMAX
 279       LOGICAL            LSAME
 280       EXTERNAL           IDAMAX, LSAME, DDOT, DNRM2
 281 *     ..
 282 *     .. External Subroutines ..
 283       EXTERNAL           DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP
 284 *     ..
 285 *     .. Executable Statements ..
 286 *
 287 *     Test the input parameters.
 288 *
 289       APPLV = LSAME( JOBV, 'A' )
 290       RSVEC = LSAME( JOBV, 'V' )
 291       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
 292          INFO = -1
 293       ELSE IF( M.LT.0 ) THEN
 294          INFO = -2
 295       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
 296          INFO = -3
 297       ELSE IF( N1.LT.0 ) THEN
 298          INFO = -4
 299       ELSE IF( LDA.LT.M ) THEN
 300          INFO = -6
 301       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
 302          INFO = -9
 303       ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
 304      $         ( APPLV.AND.( LDV.LT.MV ) )  ) THEN
 305          INFO = -11
 306       ELSE IF( TOL.LE.EPS ) THEN
 307          INFO = -14
 308       ELSE IF( NSWEEP.LT.0 ) THEN
 309          INFO = -15
 310       ELSE IF( LWORK.LT.M ) THEN
 311          INFO = -17
 312       ELSE
 313          INFO = 0
 314       END IF
 315 *
 316 *     #:(
 317       IF( INFO.NE.0 ) THEN
 318          CALL XERBLA( 'DGSVJ1', -INFO )
 319          RETURN
 320       END IF
 321 *
 322       IF( RSVEC ) THEN
 323          MVL = N
 324       ELSE IF( APPLV ) THEN
 325          MVL = MV
 326       END IF
 327       RSVEC = RSVEC .OR. APPLV
 328
 329       ROOTEPS = DSQRT( EPS )
 330       ROOTSFMIN = DSQRT( SFMIN )
 331       SMALL = SFMIN / EPS
 332       BIG = ONE / SFMIN
 333       ROOTBIG = ONE / ROOTSFMIN
 334       LARGE = BIG / DSQRT( DBLE( M*N ) )
 335       BIGTHETA = ONE / ROOTEPS
 336       ROOTTOL = DSQRT( TOL )
 337 *
 338 *     .. Initialize the right singular vector matrix ..
 339 *
 340 *     RSVEC = LSAME( JOBV, 'Y' )
 341 *
 342       EMPTSW = N1*( N-N1 )
 343       NOTROT = 0
 344       FASTR( 1 ) = ZERO
 345 *
 346 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 347 *
 348       KBL = MIN( 8, N )
 349       NBLR = N1 / KBL
 350       IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
 351
 352 *     .. the tiling is nblr-by-nblc [tiles]
 353
 354       NBLC = ( N-N1 ) / KBL
 355       IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
 356       BLSKIP = ( KBL**2 ) + 1
 357 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 358
 359       ROWSKIP = MIN( 5, KBL )
 360 *[TP] ROWSKIP is a tuning parameter.
 361       SWBAND = 0
 362 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
 363 *     if SGESVJ is used as a computational routine in the preconditioned
 364 *     Jacobi SVD algorithm SGESVJ.
 365 *
 366 *
 367 *     | *   *   * [x] [x] [x]|
 368 *     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
 369 *     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
 370 *     |[x] [x] [x] *   *   * |
 371 *     |[x] [x] [x] *   *   * |
 372 *     |[x] [x] [x] *   *   * |
 373 *
 374 *
 375       DO 1993 i = 1, NSWEEP
 376 *     .. go go go ...
 377 *
 378          MXAAPQ = ZERO
 379          MXSINJ = ZERO
 380          ISWROT = 0
 381 *
 382          NOTROT = 0
 383          PSKIPPED = 0
 384 *
 385          DO 2000 ibr = 1, NBLR
 386
 387             igl = ( ibr-1 )*KBL + 1
 388 *
 389 *
 390 *........................................................
 391 * ... go to the off diagonal blocks
 392
 393             igl = ( ibr-1 )*KBL + 1
 394
 395             DO 2010 jbc = 1, NBLC
 396
 397                jgl = N1 + ( jbc-1 )*KBL + 1
 398
 399 *        doing the block at ( ibr, jbc )
 400
 401                IJBLSK = 0
 402                DO 2100 p = igl, MIN( igl+KBL-1, N1 )
 403
 404                   AAPP = SVA( p )
 405
 406                   IF( AAPP.GT.ZERO ) THEN
 407
 408                      PSKIPPED = 0
 409
 410                      DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 411 *
 412                         AAQQ = SVA( q )
 413
 414                         IF( AAQQ.GT.ZERO ) THEN
 415                            AAPP0 = AAPP
 416 *
 417 *     .. M x 2 Jacobi SVD ..
 418 *
 419 *        .. Safe Gram matrix computation ..
 420 *
 421                            IF( AAQQ.GE.ONE ) THEN
 422                               IF( AAPP.GE.AAQQ ) THEN
 423                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
 424                               ELSE
 425                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
 426                               END IF
 427                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
 428                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
 429      $                                  q ), 1 )*D( p )*D( q ) / AAQQ )
 430      $                                  / AAPP
 431                               ELSE
 432                                  CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
 433                                  CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
 434      $                                        M, 1, WORK, LDA, IERR )
 435                                  AAPQ = DDOT( M, WORK, 1, A( 1, q ),
 436      $                                  1 )*D( q ) / AAQQ
 437                               END IF
 438                            ELSE
 439                               IF( AAPP.GE.AAQQ ) THEN
 440                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
 441                               ELSE
 442                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
 443                               END IF
 444                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
 445                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
 446      $                                  q ), 1 )*D( p )*D( q ) / AAQQ )
 447      $                                  / AAPP
 448                               ELSE
 449                                  CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
 450                                  CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
 451      $                                        M, 1, WORK, LDA, IERR )
 452                                  AAPQ = DDOT( M, WORK, 1, A( 1, p ),
 453      $                                  1 )*D( p ) / AAPP
 454                               END IF
 455                            END IF
 456
 457                            MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) )
 458
 459 *        TO rotate or NOT to rotate, THAT is the question ...
 460 *
 461                            IF( DABS( AAPQ ).GT.TOL ) THEN
 462                               NOTROT = 0
 463 *           ROTATED  = ROTATED + 1
 464                               PSKIPPED = 0
 465                               ISWROT = ISWROT + 1
 466 *
 467                               IF( ROTOK ) THEN
 468 *
 469                                  AQOAP = AAQQ / AAPP
 470                                  APOAQ = AAPP / AAQQ
 471                                  THETA = -HALF*DABS(AQOAP-APOAQ) / AAPQ
 472                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
 473
 474                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
 475                                     T = HALF / THETA
 476                                     FASTR( 3 ) = T*D( p ) / D( q )
 477                                     FASTR( 4 ) = -T*D( q ) / D( p )
 478                                     CALL DROTM( M, A( 1, p ), 1,
 479      $                                          A( 1, q ), 1, FASTR )
 480                                     IF( RSVEC )CALL DROTM( MVL,
 481      $                                              V( 1, p ), 1,
 482      $                                              V( 1, q ), 1,
 483      $                                              FASTR )
 484                                     SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
 485      $                                         ONE+T*APOAQ*AAPQ ) )
 486                                     AAPP = AAPP*DSQRT( MAX( ZERO,
 487      $                                     ONE-T*AQOAP*AAPQ ) )
 488                                     MXSINJ = MAX( MXSINJ, DABS( T ) )
 489                                  ELSE
 490 *
 491 *                 .. choose correct signum for THETA and rotate
 492 *
 493                                     THSIGN = -DSIGN( ONE, AAPQ )
 494                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
 495                                     T = ONE / ( THETA+THSIGN*
 496      $                                  DSQRT( ONE+THETA*THETA ) )
 497                                     CS = DSQRT( ONE / ( ONE+T*T ) )
 498                                     SN = T*CS
 499                                     MXSINJ = MAX( MXSINJ, DABS( SN ) )
 500                                     SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
 501      $                                         ONE+T*APOAQ*AAPQ ) )
 502                                     AAPP = AAPP*DSQRT( MAX( ZERO,
 503      $                                    ONE-T*AQOAP*AAPQ ) )
 504
 505                                     APOAQ = D( p ) / D( q )
 506                                     AQOAP = D( q ) / D( p )
 507                                     IF( D( p ).GE.ONE ) THEN
 508 *
 509                                        IF( D( q ).GE.ONE ) THEN
 510                                           FASTR( 3 ) = T*APOAQ
 511                                           FASTR( 4 ) = -T*AQOAP
 512                                           D( p ) = D( p )*CS
 513                                           D( q ) = D( q )*CS
 514                                           CALL DROTM( M, A( 1, p ), 1,
 515      $                                                A( 1, q ), 1,
 516      $                                                FASTR )
 517                                           IF( RSVEC )CALL DROTM( MVL,
 518      $                                        V( 1, p ), 1, V( 1, q ),
 519      $                                        1, FASTR )
 520                                        ELSE
 521                                           CALL DAXPY( M, -T*AQOAP,
 522      $                                                A( 1, q ), 1,
 523      $                                                A( 1, p ), 1 )
 524                                           CALL DAXPY( M, CS*SN*APOAQ,
 525      $                                                A( 1, p ), 1,
 526      $                                                A( 1, q ), 1 )
 527                                           IF( RSVEC ) THEN
 528                                              CALL DAXPY( MVL, -T*AQOAP,
 529      $                                                   V( 1, q ), 1,
 530      $                                                   V( 1, p ), 1 )
 531                                              CALL DAXPY( MVL,
 532      $                                                   CS*SN*APOAQ,
 533      $                                                   V( 1, p ), 1,
 534      $                                                   V( 1, q ), 1 )
 535                                           END IF
 536                                           D( p ) = D( p )*CS
 537                                           D( q ) = D( q ) / CS
 538                                        END IF
 539                                     ELSE
 540                                        IF( D( q ).GE.ONE ) THEN
 541                                           CALL DAXPY( M, T*APOAQ,
 542      $                                                A( 1, p ), 1,
 543      $                                                A( 1, q ), 1 )
 544                                           CALL DAXPY( M, -CS*SN*AQOAP,
 545      $                                                A( 1, q ), 1,
 546      $                                                A( 1, p ), 1 )
 547                                           IF( RSVEC ) THEN
 548                                              CALL DAXPY( MVL, T*APOAQ,
 549      $                                                   V( 1, p ), 1,
 550      $                                                   V( 1, q ), 1 )
 551                                              CALL DAXPY( MVL,
 552      $                                                   -CS*SN*AQOAP,
 553      $                                                   V( 1, q ), 1,
 554      $                                                   V( 1, p ), 1 )
 555                                           END IF
 556                                           D( p ) = D( p ) / CS
 557                                           D( q ) = D( q )*CS
 558                                        ELSE
 559                                           IF( D( p ).GE.D( q ) ) THEN
 560                                              CALL DAXPY( M, -T*AQOAP,
 561      $                                                   A( 1, q ), 1,
 562      $                                                   A( 1, p ), 1 )
 563                                              CALL DAXPY( M, CS*SN*APOAQ,
 564      $                                                   A( 1, p ), 1,
 565      $                                                   A( 1, q ), 1 )
 566                                              D( p ) = D( p )*CS
 567                                              D( q ) = D( q ) / CS
 568                                              IF( RSVEC ) THEN
 569                                                 CALL DAXPY( MVL,
 570      $                                               -T*AQOAP,
 571      $                                               V( 1, q ), 1,
 572      $                                               V( 1, p ), 1 )
 573                                                 CALL DAXPY( MVL,
 574      $                                               CS*SN*APOAQ,
 575      $                                               V( 1, p ), 1,
 576      $                                               V( 1, q ), 1 )
 577                                              END IF
 578                                           ELSE
 579                                              CALL DAXPY( M, T*APOAQ,
 580      $                                                   A( 1, p ), 1,
 581      $                                                   A( 1, q ), 1 )
 582                                              CALL DAXPY( M,
 583      $                                                   -CS*SN*AQOAP,
 584      $                                                   A( 1, q ), 1,
 585      $                                                   A( 1, p ), 1 )
 586                                              D( p ) = D( p ) / CS
 587                                              D( q ) = D( q )*CS
 588                                              IF( RSVEC ) THEN
 589                                                 CALL DAXPY( MVL,
 590      $                                               T*APOAQ, V( 1, p ),
 591      $                                               1, V( 1, q ), 1 )
 592                                                 CALL DAXPY( MVL,
 593      $                                               -CS*SN*AQOAP,
 594      $                                               V( 1, q ), 1,
 595      $                                               V( 1, p ), 1 )
 596                                              END IF
 597                                           END IF
 598                                        END IF
 599                                     END IF
 600                                  END IF
 601
 602                               ELSE
 603                                  IF( AAPP.GT.AAQQ ) THEN
 604                                     CALL DCOPY( M, A( 1, p ), 1, WORK,
 605      $                                          1 )
 606                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
 607      $                                           M, 1, WORK, LDA, IERR )
 608                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
 609      $                                           M, 1, A( 1, q ), LDA,
 610      $                                           IERR )
 611                                     TEMP1 = -AAPQ*D( p ) / D( q )
 612                                     CALL DAXPY( M, TEMP1, WORK, 1,
 613      $                                          A( 1, q ), 1 )
 614                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
 615      $                                           M, 1, A( 1, q ), LDA,
 616      $                                           IERR )
 617                                     SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
 618      $                                         ONE-AAPQ*AAPQ ) )
 619                                     MXSINJ = MAX( MXSINJ, SFMIN )
 620                                  ELSE
 621                                     CALL DCOPY( M, A( 1, q ), 1, WORK,
 622      $                                          1 )
 623                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
 624      $                                           M, 1, WORK, LDA, IERR )
 625                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
 626      $                                           M, 1, A( 1, p ), LDA,
 627      $                                           IERR )
 628                                     TEMP1 = -AAPQ*D( q ) / D( p )
 629                                     CALL DAXPY( M, TEMP1, WORK, 1,
 630      $                                          A( 1, p ), 1 )
 631                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,
 632      $                                           M, 1, A( 1, p ), LDA,
 633      $                                           IERR )
 634                                     SVA( p ) = AAPP*DSQRT( MAX( ZERO,
 635      $                                         ONE-AAPQ*AAPQ ) )
 636                                     MXSINJ = MAX( MXSINJ, SFMIN )
 637                                  END IF
 638                               END IF
 639 *           END IF ROTOK THEN ... ELSE
 640 *
 641 *           In the case of cancellation in updating SVA(q)
 642 *           .. recompute SVA(q)
 643                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
 644      $                            THEN
 645                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
 646      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
 647                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
 648      $                                         D( q )
 649                                  ELSE
 650                                     T = ZERO
 651                                     AAQQ = ONE
 652                                     CALL DLASSQ( M, A( 1, q ), 1, T,
 653      $                                           AAQQ )
 654                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
 655                                  END IF
 656                               END IF
 657                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
 658                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
 659      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
 660                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
 661      $                                     D( p )
 662                                  ELSE
 663                                     T = ZERO
 664                                     AAPP = ONE
 665                                     CALL DLASSQ( M, A( 1, p ), 1, T,
 666      $                                           AAPP )
 667                                     AAPP = T*DSQRT( AAPP )*D( p )
 668                                  END IF
 669                                  SVA( p ) = AAPP
 670                               END IF
 671 *              end of OK rotation
 672                            ELSE
 673                               NOTROT = NOTROT + 1
 674 *           SKIPPED  = SKIPPED  + 1
 675                               PSKIPPED = PSKIPPED + 1
 676                               IJBLSK = IJBLSK + 1
 677                            END IF
 678                         ELSE
 679                            NOTROT = NOTROT + 1
 680                            PSKIPPED = PSKIPPED + 1
 681                            IJBLSK = IJBLSK + 1
 682                         END IF
 683
 684 *      IF ( NOTROT .GE. EMPTSW )  GO TO 2011
 685                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
 686      $                      THEN
 687                            SVA( p ) = AAPP
 688                            NOTROT = 0
 689                            GO TO 2011
 690                         END IF
 691                         IF( ( i.LE.SWBAND ) .AND.
 692      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
 693                            AAPP = -AAPP
 694                            NOTROT = 0
 695                            GO TO 2203
 696                         END IF
 697
 698 *
 699  2200                CONTINUE
 700 *        end of the q-loop
 701  2203                CONTINUE
 702
 703                      SVA( p ) = AAPP
 704 *
 705                   ELSE
 706                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
 707      $                   MIN( jgl+KBL-1, N ) - jgl + 1
 708                      IF( AAPP.LT.ZERO )NOTROT = 0
 709 ***      IF ( NOTROT .GE. EMPTSW )  GO TO 2011
 710                   END IF
 711
 712  2100          CONTINUE
 713 *     end of the p-loop
 714  2010       CONTINUE
 715 *     end of the jbc-loop
 716  2011       CONTINUE
 717 *2011 bailed out of the jbc-loop
 718             DO 2012 p = igl, MIN( igl+KBL-1, N )
 719                SVA( p ) = DABS( SVA( p ) )
 720  2012       CONTINUE
 721 ***   IF ( NOTROT .GE. EMPTSW ) GO TO 1994
 722  2000    CONTINUE
 723 *2000 :: end of the ibr-loop
 724 *
 725 *     .. update SVA(N)
 726          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
 727      $       THEN
 728             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
 729          ELSE
 730             T = ZERO
 731             AAPP = ONE
 732             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
 733             SVA( N ) = T*DSQRT( AAPP )*D( N )
 734          END IF
 735 *
 736 *     Additional steering devices
 737 *
 738          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
 739      $       ( ISWROT.LE.N ) ) )SWBAND = i
 740
 741          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
 742      $       ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
 743             GO TO 1994
 744          END IF
 745
 746 *
 747          IF( NOTROT.GE.EMPTSW )GO TO 1994
 748
 749  1993 CONTINUE
 750 *     end i=1:NSWEEP loop
 751 * #:) Reaching this point means that the procedure has completed the given
 752 *     number of sweeps.
 753       INFO = NSWEEP - 1
 754       GO TO 1995
 755  1994 CONTINUE
 756 * #:) Reaching this point means that during the i-th sweep all pivots were
 757 *     below the given threshold, causing early exit.
 758
 759       INFO = 0
 760 * #:) INFO = 0 confirms successful iterations.
 761  1995 CONTINUE
 762 *
 763 *     Sort the vector D
 764 *
 765       DO 5991 p = 1, N - 1
 766          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
 767          IF( p.NE.q ) THEN
 768             TEMP1 = SVA( p )
 769             SVA( p ) = SVA( q )
 770             SVA( q ) = TEMP1
 771             TEMP1 = D( p )
 772             D( p ) = D( q )
 773             D( q ) = TEMP1
 774             CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
 775             IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
 776          END IF
 777  5991 CONTINUE
 778 *
 779       RETURN
 780 *     ..
 781 *     .. END OF DGSVJ1
 782 *     ..
 783       END