SRC/dgsvj0.f

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