LAPACK: SRC/dgsvj1.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
00002      +                   EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.3.0)                                    --
00005 *
00006 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
00007 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
00008 *     November 2010
00009 *
00010 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00011 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00012 *
00013 * This routine is also part of SIGMA (version 1.23, October 23. 2008.)
00014 * SIGMA is a library of algorithms for highly accurate algorithms for
00015 * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the
00016 * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
00017 *
00018       IMPLICIT           NONE
00019 *     ..
00020 *     .. Scalar Arguments ..
00021       DOUBLE PRECISION   EPS, SFMIN, TOL
00022       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
00023       CHARACTER*1        JOBV
00024 *     ..
00025 *     .. Array Arguments ..
00026       DOUBLE PRECISION   A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
00027      +                   WORK( LWORK )
00028 *     ..
00029 *
00030 *  Purpose
00031 *  =======
00032 *
00033 *  DGSVJ1 is called from SGESVJ as a pre-processor and that is its main
00034 *  purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
00035 *  it targets only particular pivots and it does not check convergence
00036 *  (stopping criterion). Few tunning parameters (marked by [TP]) are
00037 *  available for the implementer.
00038 *
00039 *  Further Details
00040 *  ~~~~~~~~~~~~~~~
00041 *  DGSVJ1 applies few sweeps of Jacobi rotations in the column space of
00042 *  the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
00043 *  off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
00044 *  block-entries (tiles) of the (1,2) off-diagonal block are marked by the
00045 *  [x]'s in the following scheme:
00046 *
00047 *     | *   *   * [x] [x] [x]|
00048 *     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
00049 *     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
00050 *     |[x] [x] [x] *   *   * |
00051 *     |[x] [x] [x] *   *   * |
00052 *     |[x] [x] [x] *   *   * |
00053 *
00054 *  In terms of the columns of A, the first N1 columns are rotated 'against'
00055 *  the remaining N-N1 columns, trying to increase the angle between the
00056 *  corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
00057 *  tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
00058 *  The number of sweeps is given in NSWEEP and the orthogonality threshold
00059 *  is given in TOL.
00060 *
00061 *  Contributors
00062 *  ~~~~~~~~~~~~
00063 *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
00064 *
00065 *  Arguments
00066 *  =========
00067 *
00068 *  JOBV    (input) CHARACTER*1
00069 *          Specifies whether the output from this procedure is used
00070 *          to compute the matrix V:
00071 *          = 'V': the product of the Jacobi rotations is accumulated
00072 *                 by postmulyiplying the N-by-N array V.
00073 *                (See the description of V.)
00074 *          = 'A': the product of the Jacobi rotations is accumulated
00075 *                 by postmulyiplying the MV-by-N array V.
00076 *                (See the descriptions of MV and V.)
00077 *          = 'N': the Jacobi rotations are not accumulated.
00078 *
00079 *  M       (input) INTEGER
00080 *          The number of rows of the input matrix A.  M >= 0.
00081 *
00082 *  N       (input) INTEGER
00083 *          The number of columns of the input matrix A.
00084 *          M >= N >= 0.
00085 *
00086 *  N1      (input) INTEGER
00087 *          N1 specifies the 2 x 2 block partition, the first N1 columns are
00088 *          rotated 'against' the remaining N-N1 columns of A.
00089 *
00090 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00091 *          On entry, M-by-N matrix A, such that A*diag(D) represents
00092 *          the input matrix.
00093 *          On exit,
00094 *          A_onexit * D_onexit represents the input matrix A*diag(D)
00095 *          post-multiplied by a sequence of Jacobi rotations, where the
00096 *          rotation threshold and the total number of sweeps are given in
00097 *          TOL and NSWEEP, respectively.
00098 *          (See the descriptions of N1, D, TOL and NSWEEP.)
00099 *
00100 *  LDA     (input) INTEGER
00101 *          The leading dimension of the array A.  LDA >= max(1,M).
00102 *
00103 *  D       (input/workspace/output) DOUBLE PRECISION array, dimension (N)
00104 *          The array D accumulates the scaling factors from the fast scaled
00105 *          Jacobi rotations.
00106 *          On entry, A*diag(D) represents the input matrix.
00107 *          On exit, A_onexit*diag(D_onexit) represents the input matrix
00108 *          post-multiplied by a sequence of Jacobi rotations, where the
00109 *          rotation threshold and the total number of sweeps are given in
00110 *          TOL and NSWEEP, respectively.
00111 *          (See the descriptions of N1, A, TOL and NSWEEP.)
00112 *
00113 *  SVA     (input/workspace/output) DOUBLE PRECISION array, dimension (N)
00114 *          On entry, SVA contains the Euclidean norms of the columns of
00115 *          the matrix A*diag(D).
00116 *          On exit, SVA contains the Euclidean norms of the columns of
00117 *          the matrix onexit*diag(D_onexit).
00118 *
00119 *  MV      (input) INTEGER
00120 *          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
00121 *                           sequence of Jacobi rotations.
00122 *          If JOBV = 'N',   then MV is not referenced.
00123 *
00124 *  V       (input/output) DOUBLE PRECISION array, dimension (LDV,N)
00125 *          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
00126 *                           sequence of Jacobi rotations.
00127 *          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
00128 *                           sequence of Jacobi rotations.
00129 *          If JOBV = 'N',   then V is not referenced.
00130 *
00131 *  LDV     (input) INTEGER
00132 *          The leading dimension of the array V,  LDV >= 1.
00133 *          If JOBV = 'V', LDV .GE. N.
00134 *          If JOBV = 'A', LDV .GE. MV.
00135 *
00136 *  EPS     (input) DOUBLE PRECISION
00137 *          EPS = DLAMCH('Epsilon')
00138 *
00139 *  SFMIN   (input) DOUBLE PRECISION
00140 *          SFMIN = DLAMCH('Safe Minimum')
00141 *
00142 *  TOL     (input) DOUBLE PRECISION
00143 *          TOL is the threshold for Jacobi rotations. For a pair
00144 *          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
00145 *          applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
00146 *
00147 *  NSWEEP  (input) INTEGER
00148 *          NSWEEP is the number of sweeps of Jacobi rotations to be
00149 *          performed.
00150 *
00151 *  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
00152 *
00153 *  LWORK   (input) INTEGER
00154 *          LWORK is the dimension of WORK. LWORK .GE. M.
00155 *
00156 *  INFO    (output) INTEGER
00157 *          = 0 : successful exit.
00158 *          < 0 : if INFO = -i, then the i-th argument had an illegal value
00159 *
00160 *  =====================================================================
00161 *
00162 *     .. Local Parameters ..
00163       DOUBLE PRECISION   ZERO, HALF, ONE, TWO
00164       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
00165      +                   TWO = 2.0D0 )
00166 *     ..
00167 *     .. Local Scalars ..
00168       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
00169      +                   BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
00170      +                   ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
00171      +                   TEMP1, THETA, THSIGN
00172       INTEGER            BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
00173      +                   ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
00174      +                   p, PSKIPPED, q, ROWSKIP, SWBAND
00175       LOGICAL            APPLV, ROTOK, RSVEC
00176 *     ..
00177 *     .. Local Arrays ..
00178       DOUBLE PRECISION   FASTR( 5 )
00179 *     ..
00180 *     .. Intrinsic Functions ..
00181       INTRINSIC          DABS, DMAX1, DBLE, MIN0, DSIGN, DSQRT
00182 *     ..
00183 *     .. External Functions ..
00184       DOUBLE PRECISION   DDOT, DNRM2
00185       INTEGER            IDAMAX
00186       LOGICAL            LSAME
00187       EXTERNAL           IDAMAX, LSAME, DDOT, DNRM2
00188 *     ..
00189 *     .. External Subroutines ..
00190       EXTERNAL           DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP
00191 *     ..
00192 *     .. Executable Statements ..
00193 *
00194 *     Test the input parameters.
00195 *
00196       APPLV = LSAME( JOBV, 'A' )
00197       RSVEC = LSAME( JOBV, 'V' )
00198       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00199          INFO = -1
00200       ELSE IF( M.LT.0 ) THEN
00201          INFO = -2
00202       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
00203          INFO = -3
00204       ELSE IF( N1.LT.0 ) THEN
00205          INFO = -4
00206       ELSE IF( LDA.LT.M ) THEN
00207          INFO = -6
00208       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
00209          INFO = -9
00210       ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
00211      &         ( APPLV.AND.( LDV.LT.MV ) )  ) THEN
00212          INFO = -11
00213       ELSE IF( TOL.LE.EPS ) THEN
00214          INFO = -14
00215       ELSE IF( NSWEEP.LT.0 ) THEN
00216          INFO = -15
00217       ELSE IF( LWORK.LT.M ) THEN
00218          INFO = -17
00219       ELSE
00220          INFO = 0
00221       END IF
00222 *
00223 *     #:(
00224       IF( INFO.NE.0 ) THEN
00225          CALL XERBLA( 'DGSVJ1', -INFO )
00226          RETURN
00227       END IF
00228 *
00229       IF( RSVEC ) THEN
00230          MVL = N
00231       ELSE IF( APPLV ) THEN
00232          MVL = MV
00233       END IF
00234       RSVEC = RSVEC .OR. APPLV
00235 
00236       ROOTEPS = DSQRT( EPS )
00237       ROOTSFMIN = DSQRT( SFMIN )
00238       SMALL = SFMIN / EPS
00239       BIG = ONE / SFMIN
00240       ROOTBIG = ONE / ROOTSFMIN
00241       LARGE = BIG / DSQRT( DBLE( M*N ) )
00242       BIGTHETA = ONE / ROOTEPS
00243       ROOTTOL = DSQRT( TOL )
00244 *
00245 *     .. Initialize the right singular vector matrix ..
00246 *
00247 *     RSVEC = LSAME( JOBV, 'Y' )
00248 *
00249       EMPTSW = N1*( N-N1 )
00250       NOTROT = 0
00251       FASTR( 1 ) = ZERO
00252 *
00253 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
00254 *
00255       KBL = MIN0( 8, N )
00256       NBLR = N1 / KBL
00257       IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
00258 
00259 *     .. the tiling is nblr-by-nblc [tiles]
00260 
00261       NBLC = ( N-N1 ) / KBL
00262       IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
00263       BLSKIP = ( KBL**2 ) + 1
00264 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
00265 
00266       ROWSKIP = MIN0( 5, KBL )
00267 *[TP] ROWSKIP is a tuning parameter.
00268       SWBAND = 0
00269 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
00270 *     if SGESVJ is used as a computational routine in the preconditioned
00271 *     Jacobi SVD algorithm SGESVJ.
00272 *
00273 *
00274 *     | *   *   * [x] [x] [x]|
00275 *     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
00276 *     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
00277 *     |[x] [x] [x] *   *   * |
00278 *     |[x] [x] [x] *   *   * |
00279 *     |[x] [x] [x] *   *   * |
00280 *
00281 *
00282       DO 1993 i = 1, NSWEEP
00283 *     .. go go go ...
00284 *
00285          MXAAPQ = ZERO
00286          MXSINJ = ZERO
00287          ISWROT = 0
00288 *
00289          NOTROT = 0
00290          PSKIPPED = 0
00291 *
00292          DO 2000 ibr = 1, NBLR
00293 
00294             igl = ( ibr-1 )*KBL + 1
00295 *
00296 *
00297 *........................................................
00298 * ... go to the off diagonal blocks
00299 
00300             igl = ( ibr-1 )*KBL + 1
00301 
00302             DO 2010 jbc = 1, NBLC
00303 
00304                jgl = N1 + ( jbc-1 )*KBL + 1
00305 
00306 *        doing the block at ( ibr, jbc )
00307 
00308                IJBLSK = 0
00309                DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
00310 
00311                   AAPP = SVA( p )
00312 
00313                   IF( AAPP.GT.ZERO ) THEN
00314 
00315                      PSKIPPED = 0
00316 
00317                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
00318 *
00319                         AAQQ = SVA( q )
00320 
00321                         IF( AAQQ.GT.ZERO ) THEN
00322                            AAPP0 = AAPP
00323 *
00324 *     .. M x 2 Jacobi SVD ..
00325 *
00326 *        .. Safe Gram matrix computation ..
00327 *
00328                            IF( AAQQ.GE.ONE ) THEN
00329                               IF( AAPP.GE.AAQQ ) THEN
00330                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
00331                               ELSE
00332                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
00333                               END IF
00334                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
00335                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
00336      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00337      +                                  / AAPP
00338                               ELSE
00339                                  CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
00340                                  CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
00341      +                                        M, 1, WORK, LDA, IERR )
00342                                  AAPQ = DDOT( M, WORK, 1, A( 1, q ),
00343      +                                  1 )*D( q ) / AAQQ
00344                               END IF
00345                            ELSE
00346                               IF( AAPP.GE.AAQQ ) THEN
00347                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
00348                               ELSE
00349                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
00350                               END IF
00351                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
00352                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
00353      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00354      +                                  / AAPP
00355                               ELSE
00356                                  CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
00357                                  CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
00358      +                                        M, 1, WORK, LDA, IERR )
00359                                  AAPQ = DDOT( M, WORK, 1, A( 1, p ),
00360      +                                  1 )*D( p ) / AAPP
00361                               END IF
00362                            END IF
00363 
00364                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
00365 
00366 *        TO rotate or NOT to rotate, THAT is the question ...
00367 *
00368                            IF( DABS( AAPQ ).GT.TOL ) THEN
00369                               NOTROT = 0
00370 *           ROTATED  = ROTATED + 1
00371                               PSKIPPED = 0
00372                               ISWROT = ISWROT + 1
00373 *
00374                               IF( ROTOK ) THEN
00375 *
00376                                  AQOAP = AAQQ / AAPP
00377                                  APOAQ = AAPP / AAQQ
00378                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /
00379      +                                   AAPQ
00380                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
00381 
00382                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
00383                                     T = HALF / THETA
00384                                     FASTR( 3 ) = T*D( p ) / D( q )
00385                                     FASTR( 4 ) = -T*D( q ) / D( p )
00386                                     CALL DROTM( M, A( 1, p ), 1,
00387      +                                          A( 1, q ), 1, FASTR )
00388                                     IF( RSVEC )CALL DROTM( MVL,
00389      +                                              V( 1, p ), 1,
00390      +                                              V( 1, q ), 1,
00391      +                                              FASTR )
00392                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
00393      +                                         ONE+T*APOAQ*AAPQ ) )
00394                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
00395      +                                     ONE-T*AQOAP*AAPQ ) )
00396                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
00397                                  ELSE
00398 *
00399 *                 .. choose correct signum for THETA and rotate
00400 *
00401                                     THSIGN = -DSIGN( ONE, AAPQ )
00402                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
00403                                     T = ONE / ( THETA+THSIGN*
00404      +                                  DSQRT( ONE+THETA*THETA ) )
00405                                     CS = DSQRT( ONE / ( ONE+T*T ) )
00406                                     SN = T*CS
00407                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
00408                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
00409      +                                         ONE+T*APOAQ*AAPQ ) )
00410                                     AAPP = AAPP*DSQRT( DMAX1( ZERO, 
00411      +                                    ONE-T*AQOAP*AAPQ ) )
00412 
00413                                     APOAQ = D( p ) / D( q )
00414                                     AQOAP = D( q ) / D( p )
00415                                     IF( D( p ).GE.ONE ) THEN
00416 *
00417                                        IF( D( q ).GE.ONE ) THEN
00418                                           FASTR( 3 ) = T*APOAQ
00419                                           FASTR( 4 ) = -T*AQOAP
00420                                           D( p ) = D( p )*CS
00421                                           D( q ) = D( q )*CS
00422                                           CALL DROTM( M, A( 1, p ), 1,
00423      +                                                A( 1, q ), 1,
00424      +                                                FASTR )
00425                                           IF( RSVEC )CALL DROTM( MVL,
00426      +                                        V( 1, p ), 1, V( 1, q ),
00427      +                                        1, FASTR )
00428                                        ELSE
00429                                           CALL DAXPY( M, -T*AQOAP,
00430      +                                                A( 1, q ), 1,
00431      +                                                A( 1, p ), 1 )
00432                                           CALL DAXPY( M, CS*SN*APOAQ,
00433      +                                                A( 1, p ), 1,
00434      +                                                A( 1, q ), 1 )
00435                                           IF( RSVEC ) THEN
00436                                              CALL DAXPY( MVL, -T*AQOAP,
00437      +                                                   V( 1, q ), 1,
00438      +                                                   V( 1, p ), 1 )
00439                                              CALL DAXPY( MVL,
00440      +                                                   CS*SN*APOAQ,
00441      +                                                   V( 1, p ), 1,
00442      +                                                   V( 1, q ), 1 )
00443                                           END IF
00444                                           D( p ) = D( p )*CS
00445                                           D( q ) = D( q ) / CS
00446                                        END IF
00447                                     ELSE
00448                                        IF( D( q ).GE.ONE ) THEN
00449                                           CALL DAXPY( M, T*APOAQ,
00450      +                                                A( 1, p ), 1,
00451      +                                                A( 1, q ), 1 )
00452                                           CALL DAXPY( M, -CS*SN*AQOAP,
00453      +                                                A( 1, q ), 1,
00454      +                                                A( 1, p ), 1 )
00455                                           IF( RSVEC ) THEN
00456                                              CALL DAXPY( MVL, T*APOAQ,
00457      +                                                   V( 1, p ), 1,
00458      +                                                   V( 1, q ), 1 )
00459                                              CALL DAXPY( MVL,
00460      +                                                   -CS*SN*AQOAP,
00461      +                                                   V( 1, q ), 1,
00462      +                                                   V( 1, p ), 1 )
00463                                           END IF
00464                                           D( p ) = D( p ) / CS
00465                                           D( q ) = D( q )*CS
00466                                        ELSE
00467                                           IF( D( p ).GE.D( q ) ) THEN
00468                                              CALL DAXPY( M, -T*AQOAP,
00469      +                                                   A( 1, q ), 1,
00470      +                                                   A( 1, p ), 1 )
00471                                              CALL DAXPY( M, CS*SN*APOAQ,
00472      +                                                   A( 1, p ), 1,
00473      +                                                   A( 1, q ), 1 )
00474                                              D( p ) = D( p )*CS
00475                                              D( q ) = D( q ) / CS
00476                                              IF( RSVEC ) THEN
00477                                                 CALL DAXPY( MVL,
00478      +                                               -T*AQOAP,
00479      +                                               V( 1, q ), 1,
00480      +                                               V( 1, p ), 1 )
00481                                                 CALL DAXPY( MVL,
00482      +                                               CS*SN*APOAQ,
00483      +                                               V( 1, p ), 1,
00484      +                                               V( 1, q ), 1 )
00485                                              END IF
00486                                           ELSE
00487                                              CALL DAXPY( M, T*APOAQ,
00488      +                                                   A( 1, p ), 1,
00489      +                                                   A( 1, q ), 1 )
00490                                              CALL DAXPY( M,
00491      +                                                   -CS*SN*AQOAP,
00492      +                                                   A( 1, q ), 1,
00493      +                                                   A( 1, p ), 1 )
00494                                              D( p ) = D( p ) / CS
00495                                              D( q ) = D( q )*CS
00496                                              IF( RSVEC ) THEN
00497                                                 CALL DAXPY( MVL,
00498      +                                               T*APOAQ, V( 1, p ),
00499      +                                               1, V( 1, q ), 1 )
00500                                                 CALL DAXPY( MVL,
00501      +                                               -CS*SN*AQOAP,
00502      +                                               V( 1, q ), 1,
00503      +                                               V( 1, p ), 1 )
00504                                              END IF
00505                                           END IF
00506                                        END IF
00507                                     END IF
00508                                  END IF
00509 
00510                               ELSE
00511                                  IF( AAPP.GT.AAQQ ) THEN
00512                                     CALL DCOPY( M, A( 1, p ), 1, WORK,
00513      +                                          1 )
00514                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
00515      +                                           M, 1, WORK, LDA, IERR )
00516                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
00517      +                                           M, 1, A( 1, q ), LDA,
00518      +                                           IERR )
00519                                     TEMP1 = -AAPQ*D( p ) / D( q )
00520                                     CALL DAXPY( M, TEMP1, WORK, 1,
00521      +                                          A( 1, q ), 1 )
00522                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
00523      +                                           M, 1, A( 1, q ), LDA,
00524      +                                           IERR )
00525                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
00526      +                                         ONE-AAPQ*AAPQ ) )
00527                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
00528                                  ELSE
00529                                     CALL DCOPY( M, A( 1, q ), 1, WORK,
00530      +                                          1 )
00531                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
00532      +                                           M, 1, WORK, LDA, IERR )
00533                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
00534      +                                           M, 1, A( 1, p ), LDA,
00535      +                                           IERR )
00536                                     TEMP1 = -AAPQ*D( q ) / D( p )
00537                                     CALL DAXPY( M, TEMP1, WORK, 1,
00538      +                                          A( 1, p ), 1 )
00539                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,
00540      +                                           M, 1, A( 1, p ), LDA,
00541      +                                           IERR )
00542                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
00543      +                                         ONE-AAPQ*AAPQ ) )
00544                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
00545                                  END IF
00546                               END IF
00547 *           END IF ROTOK THEN ... ELSE
00548 *
00549 *           In the case of cancellation in updating SVA(q)
00550 *           .. recompute SVA(q)
00551                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
00552      +                            THEN
00553                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
00554      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
00555                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
00556      +                                         D( q )
00557                                  ELSE
00558                                     T = ZERO
00559                                     AAQQ = ONE
00560                                     CALL DLASSQ( M, A( 1, q ), 1, T,
00561      +                                           AAQQ )
00562                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
00563                                  END IF
00564                               END IF
00565                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
00566                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
00567      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
00568                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
00569      +                                     D( p )
00570                                  ELSE
00571                                     T = ZERO
00572                                     AAPP = ONE
00573                                     CALL DLASSQ( M, A( 1, p ), 1, T,
00574      +                                           AAPP )
00575                                     AAPP = T*DSQRT( AAPP )*D( p )
00576                                  END IF
00577                                  SVA( p ) = AAPP
00578                               END IF
00579 *              end of OK rotation
00580                            ELSE
00581                               NOTROT = NOTROT + 1
00582 *           SKIPPED  = SKIPPED  + 1
00583                               PSKIPPED = PSKIPPED + 1
00584                               IJBLSK = IJBLSK + 1
00585                            END IF
00586                         ELSE
00587                            NOTROT = NOTROT + 1
00588                            PSKIPPED = PSKIPPED + 1
00589                            IJBLSK = IJBLSK + 1
00590                         END IF
00591 
00592 *      IF ( NOTROT .GE. EMPTSW )  GO TO 2011
00593                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
00594      +                      THEN
00595                            SVA( p ) = AAPP
00596                            NOTROT = 0
00597                            GO TO 2011
00598                         END IF
00599                         IF( ( i.LE.SWBAND ) .AND.
00600      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
00601                            AAPP = -AAPP
00602                            NOTROT = 0
00603                            GO TO 2203
00604                         END IF
00605 
00606 *
00607  2200                CONTINUE
00608 *        end of the q-loop
00609  2203                CONTINUE
00610 
00611                      SVA( p ) = AAPP
00612 *
00613                   ELSE
00614                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
00615      +                   MIN0( jgl+KBL-1, N ) - jgl + 1
00616                      IF( AAPP.LT.ZERO )NOTROT = 0
00617 ***      IF ( NOTROT .GE. EMPTSW )  GO TO 2011
00618                   END IF
00619 
00620  2100          CONTINUE
00621 *     end of the p-loop
00622  2010       CONTINUE
00623 *     end of the jbc-loop
00624  2011       CONTINUE
00625 *2011 bailed out of the jbc-loop
00626             DO 2012 p = igl, MIN0( igl+KBL-1, N )
00627                SVA( p ) = DABS( SVA( p ) )
00628  2012       CONTINUE
00629 ***   IF ( NOTROT .GE. EMPTSW ) GO TO 1994
00630  2000    CONTINUE
00631 *2000 :: end of the ibr-loop
00632 *
00633 *     .. update SVA(N)
00634          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
00635      +       THEN
00636             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
00637          ELSE
00638             T = ZERO
00639             AAPP = ONE
00640             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
00641             SVA( N ) = T*DSQRT( AAPP )*D( N )
00642          END IF
00643 *
00644 *     Additional steering devices
00645 *
00646          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
00647      +       ( ISWROT.LE.N ) ) )SWBAND = i
00648 
00649          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
00650      +       ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
00651             GO TO 1994
00652          END IF
00653 
00654 *
00655          IF( NOTROT.GE.EMPTSW )GO TO 1994
00656 
00657  1993 CONTINUE
00658 *     end i=1:NSWEEP loop
00659 * #:) Reaching this point means that the procedure has completed the given
00660 *     number of sweeps.
00661       INFO = NSWEEP - 1
00662       GO TO 1995
00663  1994 CONTINUE
00664 * #:) Reaching this point means that during the i-th sweep all pivots were
00665 *     below the given threshold, causing early exit.
00666 
00667       INFO = 0
00668 * #:) INFO = 0 confirms successful iterations.
00669  1995 CONTINUE
00670 *
00671 *     Sort the vector D
00672 *
00673       DO 5991 p = 1, N - 1
00674          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
00675          IF( p.NE.q ) THEN
00676             TEMP1 = SVA( p )
00677             SVA( p ) = SVA( q )
00678             SVA( q ) = TEMP1
00679             TEMP1 = D( p )
00680             D( p ) = D( q )
00681             D( q ) = TEMP1
00682             CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
00683             IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
00684          END IF
00685  5991 CONTINUE
00686 *
00687       RETURN
00688 *     ..
00689 *     .. END OF DGSVJ1
00690 *     ..
00691       END