LAPACK: SRC/sgsvj0.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
00002      +                   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       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
00022       REAL               EPS, SFMIN, TOL
00023       CHARACTER*1        JOBV
00024 *     ..
00025 *     .. Array Arguments ..
00026       REAL               A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
00027      +                   WORK( LWORK )
00028 *     ..
00029 *
00030 *  Purpose
00031 *  =======
00032 *
00033 *  SGSVJ0 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 does not check convergence (stopping criterion). Few tuning
00036 *  parameters (marked by [TP]) are available for the implementer.
00037 *
00038 *  Further Details
00039 *  ~~~~~~~~~~~~~~~
00040 *  SGSVJ0 is used just to enable SGESVJ to call a simplified version of
00041 *  itself to work on a submatrix of the original matrix.
00042 *
00043 *  Contributors
00044 *  ~~~~~~~~~~~~
00045 *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
00046 *
00047 *  Bugs, Examples and Comments
00048 *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
00049 *  Please report all bugs and send interesting test examples and comments to
00050 *  drmac@math.hr. Thank you.
00051 *
00052 *  Arguments
00053 *  =========
00054 *
00055 *  JOBV    (input) CHARACTER*1
00056 *          Specifies whether the output from this procedure is used
00057 *          to compute the matrix V:
00058 *          = 'V': the product of the Jacobi rotations is accumulated
00059 *                 by postmulyiplying the N-by-N array V.
00060 *                (See the description of V.)
00061 *          = 'A': the product of the Jacobi rotations is accumulated
00062 *                 by postmulyiplying the MV-by-N array V.
00063 *                (See the descriptions of MV and V.)
00064 *          = 'N': the Jacobi rotations are not accumulated.
00065 *
00066 *  M       (input) INTEGER
00067 *          The number of rows of the input matrix A.  M >= 0.
00068 *
00069 *  N       (input) INTEGER
00070 *          The number of columns of the input matrix A.
00071 *          M >= N >= 0.
00072 *
00073 *  A       (input/output) REAL array, dimension (LDA,N)
00074 *          On entry, M-by-N matrix A, such that A*diag(D) represents
00075 *          the input matrix.
00076 *          On exit,
00077 *          A_onexit * D_onexit represents the input matrix A*diag(D)
00078 *          post-multiplied by a sequence of Jacobi rotations, where the
00079 *          rotation threshold and the total number of sweeps are given in
00080 *          TOL and NSWEEP, respectively.
00081 *          (See the descriptions of D, TOL and NSWEEP.)
00082 *
00083 *  LDA     (input) INTEGER
00084 *          The leading dimension of the array A.  LDA >= max(1,M).
00085 *
00086 *  D       (input/workspace/output) REAL array, dimension (N)
00087 *          The array D accumulates the scaling factors from the fast scaled
00088 *          Jacobi rotations.
00089 *          On entry, A*diag(D) represents the input matrix.
00090 *          On exit, A_onexit*diag(D_onexit) represents the input matrix
00091 *          post-multiplied by a sequence of Jacobi rotations, where the
00092 *          rotation threshold and the total number of sweeps are given in
00093 *          TOL and NSWEEP, respectively.
00094 *          (See the descriptions of A, TOL and NSWEEP.)
00095 *
00096 *  SVA     (input/workspace/output) REAL array, dimension (N)
00097 *          On entry, SVA contains the Euclidean norms of the columns of
00098 *          the matrix A*diag(D).
00099 *          On exit, SVA contains the Euclidean norms of the columns of
00100 *          the matrix onexit*diag(D_onexit).
00101 *
00102 *  MV      (input) INTEGER
00103 *          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
00104 *                           sequence of Jacobi rotations.
00105 *          If JOBV = 'N',   then MV is not referenced.
00106 *
00107 *  V       (input/output) REAL array, dimension (LDV,N)
00108 *          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
00109 *                           sequence of Jacobi rotations.
00110 *          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
00111 *                           sequence of Jacobi rotations.
00112 *          If JOBV = 'N',   then V is not referenced.
00113 *
00114 *  LDV     (input) INTEGER
00115 *          The leading dimension of the array V,  LDV >= 1.
00116 *          If JOBV = 'V', LDV .GE. N.
00117 *          If JOBV = 'A', LDV .GE. MV.
00118 *
00119 *  EPS     (input) INTEGER
00120 *          EPS = SLAMCH('Epsilon')
00121 *
00122 *  SFMIN   (input) INTEGER
00123 *          SFMIN = SLAMCH('Safe Minimum')
00124 *
00125 *  TOL     (input) REAL
00126 *          TOL is the threshold for Jacobi rotations. For a pair
00127 *          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
00128 *          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
00129 *
00130 *  NSWEEP  (input) INTEGER
00131 *          NSWEEP is the number of sweeps of Jacobi rotations to be
00132 *          performed.
00133 *
00134 *  WORK    (workspace) REAL array, dimension LWORK.
00135 *
00136 *  LWORK   (input) INTEGER
00137 *          LWORK is the dimension of WORK. LWORK .GE. M.
00138 *
00139 *  INFO    (output) INTEGER
00140 *          = 0 : successful exit.
00141 *          < 0 : if INFO = -i, then the i-th argument had an illegal value
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Local Parameters ..
00146       REAL               ZERO, HALF, ONE, TWO
00147       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
00148      +                   TWO = 2.0E0 )
00149 *     ..
00150 *     .. Local Scalars ..
00151       REAL               AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
00152      +                   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
00153      +                   ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
00154      +                   THSIGN
00155       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
00156      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
00157      +                   NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
00158       LOGICAL            APPLV, ROTOK, RSVEC
00159 *     ..
00160 *     .. Local Arrays ..
00161       REAL               FASTR( 5 )
00162 *     ..
00163 *     .. Intrinsic Functions ..
00164       INTRINSIC          ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT
00165 *     ..
00166 *     .. External Functions ..
00167       REAL               SDOT, SNRM2
00168       INTEGER            ISAMAX
00169       LOGICAL            LSAME
00170       EXTERNAL           ISAMAX, LSAME, SDOT, SNRM2
00171 *     ..
00172 *     .. External Subroutines ..
00173       EXTERNAL           SAXPY, SCOPY, SLASCL, SLASSQ, SROTM, SSWAP
00174 *     ..
00175 *     .. Executable Statements ..
00176 *
00177 *     Test the input parameters.
00178 *
00179       APPLV = LSAME( JOBV, 'A' )
00180       RSVEC = LSAME( JOBV, 'V' )
00181       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00182          INFO = -1
00183       ELSE IF( M.LT.0 ) THEN
00184          INFO = -2
00185       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
00186          INFO = -3
00187       ELSE IF( LDA.LT.M ) THEN
00188          INFO = -5
00189       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
00190          INFO = -8
00191       ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
00192      &         ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
00193          INFO = -10
00194       ELSE IF( TOL.LE.EPS ) THEN
00195          INFO = -13
00196       ELSE IF( NSWEEP.LT.0 ) THEN
00197          INFO = -14
00198       ELSE IF( LWORK.LT.M ) THEN
00199          INFO = -16
00200       ELSE
00201          INFO = 0
00202       END IF
00203 *
00204 *     #:(
00205       IF( INFO.NE.0 ) THEN
00206          CALL XERBLA( 'SGSVJ0', -INFO )
00207          RETURN
00208       END IF
00209 *
00210       IF( RSVEC ) THEN
00211          MVL = N
00212       ELSE IF( APPLV ) THEN
00213          MVL = MV
00214       END IF
00215       RSVEC = RSVEC .OR. APPLV
00216 
00217       ROOTEPS = SQRT( EPS )
00218       ROOTSFMIN = SQRT( SFMIN )
00219       SMALL = SFMIN / EPS
00220       BIG = ONE / SFMIN
00221       ROOTBIG = ONE / ROOTSFMIN
00222       BIGTHETA = ONE / ROOTEPS
00223       ROOTTOL = SQRT( TOL )
00224 *
00225 *
00226 *     .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
00227 *
00228       EMPTSW = ( N*( N-1 ) ) / 2
00229       NOTROT = 0
00230       FASTR( 1 ) = ZERO
00231 *
00232 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
00233 *
00234 
00235       SWBAND = 0
00236 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
00237 *     if SGESVJ is used as a computational routine in the preconditioned
00238 *     Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
00239 *     ......
00240 
00241       KBL = MIN0( 8, N )
00242 *[TP] KBL is a tuning parameter that defines the tile size in the
00243 *     tiling of the p-q loops of pivot pairs. In general, an optimal
00244 *     value of KBL depends on the matrix dimensions and on the
00245 *     parameters of the computer's memory.
00246 *
00247       NBL = N / KBL
00248       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
00249 
00250       BLSKIP = ( KBL**2 ) + 1
00251 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
00252 
00253       ROWSKIP = MIN0( 5, KBL )
00254 *[TP] ROWSKIP is a tuning parameter.
00255 
00256       LKAHEAD = 1
00257 *[TP] LKAHEAD is a tuning parameter.
00258       SWBAND = 0
00259       PSKIPPED = 0
00260 *
00261       DO 1993 i = 1, NSWEEP
00262 *     .. go go go ...
00263 *
00264          MXAAPQ = ZERO
00265          MXSINJ = ZERO
00266          ISWROT = 0
00267 *
00268          NOTROT = 0
00269          PSKIPPED = 0
00270 *
00271          DO 2000 ibr = 1, NBL
00272 
00273             igl = ( ibr-1 )*KBL + 1
00274 *
00275             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
00276 *
00277                igl = igl + ir1*KBL
00278 *
00279                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
00280 
00281 *     .. de Rijk's pivoting
00282                   q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
00283                   IF( p.NE.q ) THEN
00284                      CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
00285                      IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1,
00286      +                                      V( 1, q ), 1 )
00287                      TEMP1 = SVA( p )
00288                      SVA( p ) = SVA( q )
00289                      SVA( q ) = TEMP1
00290                      TEMP1 = D( p )
00291                      D( p ) = D( q )
00292                      D( q ) = TEMP1
00293                   END IF
00294 *
00295                   IF( ir1.EQ.0 ) THEN
00296 *
00297 *        Column norms are periodically updated by explicit
00298 *        norm computation.
00299 *        Caveat:
00300 *        Some BLAS implementations compute SNRM2(M,A(1,p),1)
00301 *        as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may result in
00302 *        overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and
00303 *        undeflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
00304 *        Hence, SNRM2 cannot be trusted, not even in the case when
00305 *        the true norm is far from the under(over)flow boundaries.
00306 *        If properly implemented SNRM2 is available, the IF-THEN-ELSE
00307 *        below should read "AAPP = SNRM2( M, A(1,p), 1 ) * D(p)".
00308 *
00309                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
00310      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
00311                         SVA( p ) = SNRM2( M, A( 1, p ), 1 )*D( p )
00312                      ELSE
00313                         TEMP1 = ZERO
00314                         AAPP = ONE
00315                         CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
00316                         SVA( p ) = TEMP1*SQRT( AAPP )*D( p )
00317                      END IF
00318                      AAPP = SVA( p )
00319                   ELSE
00320                      AAPP = SVA( p )
00321                   END IF
00322 
00323 *
00324                   IF( AAPP.GT.ZERO ) THEN
00325 *
00326                      PSKIPPED = 0
00327 *
00328                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
00329 *
00330                         AAQQ = SVA( q )
00331 
00332                         IF( AAQQ.GT.ZERO ) THEN
00333 *
00334                            AAPP0 = AAPP
00335                            IF( AAQQ.GE.ONE ) THEN
00336                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
00337                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
00338                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00339      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00340      +                                  / AAPP
00341                               ELSE
00342                                  CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
00343                                  CALL SLASCL( 'G', 0, 0, AAPP, D( p ),
00344      +                                        M, 1, WORK, LDA, IERR )
00345                                  AAPQ = SDOT( M, WORK, 1, A( 1, q ),
00346      +                                  1 )*D( q ) / AAQQ
00347                               END IF
00348                            ELSE
00349                               ROTOK = AAPP.LE.( AAQQ / SMALL )
00350                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
00351                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00352      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00353      +                                  / AAPP
00354                               ELSE
00355                                  CALL SCOPY( M, A( 1, q ), 1, WORK, 1 )
00356                                  CALL SLASCL( 'G', 0, 0, AAQQ, D( q ),
00357      +                                        M, 1, WORK, LDA, IERR )
00358                                  AAPQ = SDOT( M, WORK, 1, A( 1, p ),
00359      +                                  1 )*D( p ) / AAPP
00360                               END IF
00361                            END IF
00362 *
00363                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
00364 *
00365 *        TO rotate or NOT to rotate, THAT is the question ...
00366 *
00367                            IF( ABS( AAPQ ).GT.TOL ) THEN
00368 *
00369 *           .. rotate
00370 *           ROTATED = ROTATED + ONE
00371 *
00372                               IF( ir1.EQ.0 ) THEN
00373                                  NOTROT = 0
00374                                  PSKIPPED = 0
00375                                  ISWROT = ISWROT + 1
00376                               END IF
00377 *
00378                               IF( ROTOK ) THEN
00379 *
00380                                  AQOAP = AAQQ / AAPP
00381                                  APOAQ = AAPP / AAQQ
00382                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
00383 *
00384                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
00385 *
00386                                     T = HALF / THETA
00387                                     FASTR( 3 ) = T*D( p ) / D( q )
00388                                     FASTR( 4 ) = -T*D( q ) / D( p )
00389                                     CALL SROTM( M, A( 1, p ), 1,
00390      +                                          A( 1, q ), 1, FASTR )
00391                                     IF( RSVEC )CALL SROTM( MVL,
00392      +                                              V( 1, p ), 1,
00393      +                                              V( 1, q ), 1,
00394      +                                              FASTR )
00395                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00396      +                                         ONE+T*APOAQ*AAPQ ) )
00397                                     AAPP = AAPP*SQRT( AMAX1( ZERO, 
00398      +                                         ONE-T*AQOAP*AAPQ ) )
00399                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
00400 *
00401                                  ELSE
00402 *
00403 *                 .. choose correct signum for THETA and rotate
00404 *
00405                                     THSIGN = -SIGN( ONE, AAPQ )
00406                                     T = ONE / ( THETA+THSIGN*
00407      +                                  SQRT( ONE+THETA*THETA ) )
00408                                     CS = SQRT( ONE / ( ONE+T*T ) )
00409                                     SN = T*CS
00410 *
00411                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
00412                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00413      +                                         ONE+T*APOAQ*AAPQ ) )
00414                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
00415      +                                     ONE-T*AQOAP*AAPQ ) )
00416 *
00417                                     APOAQ = D( p ) / D( q )
00418                                     AQOAP = D( q ) / D( p )
00419                                     IF( D( p ).GE.ONE ) THEN
00420                                        IF( D( q ).GE.ONE ) THEN
00421                                           FASTR( 3 ) = T*APOAQ
00422                                           FASTR( 4 ) = -T*AQOAP
00423                                           D( p ) = D( p )*CS
00424                                           D( q ) = D( q )*CS
00425                                           CALL SROTM( M, A( 1, p ), 1,
00426      +                                                A( 1, q ), 1,
00427      +                                                FASTR )
00428                                           IF( RSVEC )CALL SROTM( MVL,
00429      +                                        V( 1, p ), 1, V( 1, q ),
00430      +                                        1, FASTR )
00431                                        ELSE
00432                                           CALL SAXPY( M, -T*AQOAP,
00433      +                                                A( 1, q ), 1,
00434      +                                                A( 1, p ), 1 )
00435                                           CALL SAXPY( M, CS*SN*APOAQ,
00436      +                                                A( 1, p ), 1,
00437      +                                                A( 1, q ), 1 )
00438                                           D( p ) = D( p )*CS
00439                                           D( q ) = D( q ) / CS
00440                                           IF( RSVEC ) THEN
00441                                              CALL SAXPY( MVL, -T*AQOAP,
00442      +                                                   V( 1, q ), 1,
00443      +                                                   V( 1, p ), 1 )
00444                                              CALL SAXPY( MVL,
00445      +                                                   CS*SN*APOAQ,
00446      +                                                   V( 1, p ), 1,
00447      +                                                   V( 1, q ), 1 )
00448                                           END IF
00449                                        END IF
00450                                     ELSE
00451                                        IF( D( q ).GE.ONE ) THEN
00452                                           CALL SAXPY( M, T*APOAQ,
00453      +                                                A( 1, p ), 1,
00454      +                                                A( 1, q ), 1 )
00455                                           CALL SAXPY( M, -CS*SN*AQOAP,
00456      +                                                A( 1, q ), 1,
00457      +                                                A( 1, p ), 1 )
00458                                           D( p ) = D( p ) / CS
00459                                           D( q ) = D( q )*CS
00460                                           IF( RSVEC ) THEN
00461                                              CALL SAXPY( MVL, T*APOAQ,
00462      +                                                   V( 1, p ), 1,
00463      +                                                   V( 1, q ), 1 )
00464                                              CALL SAXPY( MVL,
00465      +                                                   -CS*SN*AQOAP,
00466      +                                                   V( 1, q ), 1,
00467      +                                                   V( 1, p ), 1 )
00468                                           END IF
00469                                        ELSE
00470                                           IF( D( p ).GE.D( q ) ) THEN
00471                                              CALL SAXPY( M, -T*AQOAP,
00472      +                                                   A( 1, q ), 1,
00473      +                                                   A( 1, p ), 1 )
00474                                              CALL SAXPY( M, CS*SN*APOAQ,
00475      +                                                   A( 1, p ), 1,
00476      +                                                   A( 1, q ), 1 )
00477                                              D( p ) = D( p )*CS
00478                                              D( q ) = D( q ) / CS
00479                                              IF( RSVEC ) THEN
00480                                                 CALL SAXPY( MVL,
00481      +                                               -T*AQOAP,
00482      +                                               V( 1, q ), 1,
00483      +                                               V( 1, p ), 1 )
00484                                                 CALL SAXPY( MVL,
00485      +                                               CS*SN*APOAQ,
00486      +                                               V( 1, p ), 1,
00487      +                                               V( 1, q ), 1 )
00488                                              END IF
00489                                           ELSE
00490                                              CALL SAXPY( M, T*APOAQ,
00491      +                                                   A( 1, p ), 1,
00492      +                                                   A( 1, q ), 1 )
00493                                              CALL SAXPY( M,
00494      +                                                   -CS*SN*AQOAP,
00495      +                                                   A( 1, q ), 1,
00496      +                                                   A( 1, p ), 1 )
00497                                              D( p ) = D( p ) / CS
00498                                              D( q ) = D( q )*CS
00499                                              IF( RSVEC ) THEN
00500                                                 CALL SAXPY( MVL,
00501      +                                               T*APOAQ, V( 1, p ),
00502      +                                               1, V( 1, q ), 1 )
00503                                                 CALL SAXPY( MVL,
00504      +                                               -CS*SN*AQOAP,
00505      +                                               V( 1, q ), 1,
00506      +                                               V( 1, p ), 1 )
00507                                              END IF
00508                                           END IF
00509                                        END IF
00510                                     END IF
00511                                  END IF
00512 *
00513                               ELSE
00514 *              .. have to use modified Gram-Schmidt like transformation
00515                                  CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
00516                                  CALL SLASCL( 'G', 0, 0, AAPP, ONE, M,
00517      +                                        1, WORK, LDA, IERR )
00518                                  CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M,
00519      +                                        1, A( 1, q ), LDA, IERR )
00520                                  TEMP1 = -AAPQ*D( p ) / D( q )
00521                                  CALL SAXPY( M, TEMP1, WORK, 1,
00522      +                                       A( 1, q ), 1 )
00523                                  CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M,
00524      +                                        1, A( 1, q ), LDA, IERR )
00525                                  SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00526      +                                      ONE-AAPQ*AAPQ ) )
00527                                  MXSINJ = AMAX1( MXSINJ, SFMIN )
00528                               END IF
00529 *           END IF ROTOK THEN ... ELSE
00530 *
00531 *           In the case of cancellation in updating SVA(q), SVA(p)
00532 *           recompute SVA(q), SVA(p).
00533                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
00534      +                            THEN
00535                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
00536      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
00537                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
00538      +                                         D( q )
00539                                  ELSE
00540                                     T = ZERO
00541                                     AAQQ = ONE
00542                                     CALL SLASSQ( M, A( 1, q ), 1, T,
00543      +                                           AAQQ )
00544                                     SVA( q ) = T*SQRT( AAQQ )*D( q )
00545                                  END IF
00546                               END IF
00547                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
00548                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
00549      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
00550                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
00551      +                                     D( p )
00552                                  ELSE
00553                                     T = ZERO
00554                                     AAPP = ONE
00555                                     CALL SLASSQ( M, A( 1, p ), 1, T,
00556      +                                           AAPP )
00557                                     AAPP = T*SQRT( AAPP )*D( p )
00558                                  END IF
00559                                  SVA( p ) = AAPP
00560                               END IF
00561 *
00562                            ELSE
00563 *        A(:,p) and A(:,q) already numerically orthogonal
00564                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
00565                               PSKIPPED = PSKIPPED + 1
00566                            END IF
00567                         ELSE
00568 *        A(:,q) is zero column
00569                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
00570                            PSKIPPED = PSKIPPED + 1
00571                         END IF
00572 *
00573                         IF( ( i.LE.SWBAND ) .AND.
00574      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
00575                            IF( ir1.EQ.0 )AAPP = -AAPP
00576                            NOTROT = 0
00577                            GO TO 2103
00578                         END IF
00579 *
00580  2002                CONTINUE
00581 *     END q-LOOP
00582 *
00583  2103                CONTINUE
00584 *     bailed out of q-loop
00585 
00586                      SVA( p ) = AAPP
00587 
00588                   ELSE
00589                      SVA( p ) = AAPP
00590                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
00591      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
00592                   END IF
00593 *
00594  2001          CONTINUE
00595 *     end of the p-loop
00596 *     end of doing the block ( ibr, ibr )
00597  1002       CONTINUE
00598 *     end of ir1-loop
00599 *
00600 *........................................................
00601 * ... go to the off diagonal blocks
00602 *
00603             igl = ( ibr-1 )*KBL + 1
00604 *
00605             DO 2010 jbc = ibr + 1, NBL
00606 *
00607                jgl = ( jbc-1 )*KBL + 1
00608 *
00609 *        doing the block at ( ibr, jbc )
00610 *
00611                IJBLSK = 0
00612                DO 2100 p = igl, MIN0( igl+KBL-1, N )
00613 *
00614                   AAPP = SVA( p )
00615 *
00616                   IF( AAPP.GT.ZERO ) THEN
00617 *
00618                      PSKIPPED = 0
00619 *
00620                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
00621 *
00622                         AAQQ = SVA( q )
00623 *
00624                         IF( AAQQ.GT.ZERO ) THEN
00625                            AAPP0 = AAPP
00626 *
00627 *     .. M x 2 Jacobi SVD ..
00628 *
00629 *        .. Safe Gram matrix computation ..
00630 *
00631                            IF( AAQQ.GE.ONE ) THEN
00632                               IF( AAPP.GE.AAQQ ) THEN
00633                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
00634                               ELSE
00635                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
00636                               END IF
00637                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
00638                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00639      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00640      +                                  / AAPP
00641                               ELSE
00642                                  CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
00643                                  CALL SLASCL( 'G', 0, 0, AAPP, D( p ),
00644      +                                        M, 1, WORK, LDA, IERR )
00645                                  AAPQ = SDOT( M, WORK, 1, A( 1, q ),
00646      +                                  1 )*D( q ) / AAQQ
00647                               END IF
00648                            ELSE
00649                               IF( AAPP.GE.AAQQ ) THEN
00650                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
00651                               ELSE
00652                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
00653                               END IF
00654                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
00655                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00656      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00657      +                                  / AAPP
00658                               ELSE
00659                                  CALL SCOPY( M, A( 1, q ), 1, WORK, 1 )
00660                                  CALL SLASCL( 'G', 0, 0, AAQQ, D( q ),
00661      +                                        M, 1, WORK, LDA, IERR )
00662                                  AAPQ = SDOT( M, WORK, 1, A( 1, p ),
00663      +                                  1 )*D( p ) / AAPP
00664                               END IF
00665                            END IF
00666 *
00667                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
00668 *
00669 *        TO rotate or NOT to rotate, THAT is the question ...
00670 *
00671                            IF( ABS( AAPQ ).GT.TOL ) THEN
00672                               NOTROT = 0
00673 *           ROTATED  = ROTATED + 1
00674                               PSKIPPED = 0
00675                               ISWROT = ISWROT + 1
00676 *
00677                               IF( ROTOK ) THEN
00678 *
00679                                  AQOAP = AAQQ / AAPP
00680                                  APOAQ = AAPP / AAQQ
00681                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
00682                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
00683 *
00684                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
00685                                     T = HALF / THETA
00686                                     FASTR( 3 ) = T*D( p ) / D( q )
00687                                     FASTR( 4 ) = -T*D( q ) / D( p )
00688                                     CALL SROTM( M, A( 1, p ), 1,
00689      +                                          A( 1, q ), 1, FASTR )
00690                                     IF( RSVEC )CALL SROTM( MVL,
00691      +                                              V( 1, p ), 1,
00692      +                                              V( 1, q ), 1,
00693      +                                              FASTR )
00694                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00695      +                                         ONE+T*APOAQ*AAPQ ) )
00696                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
00697      +                                     ONE-T*AQOAP*AAPQ ) )
00698                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
00699                                  ELSE
00700 *
00701 *                 .. choose correct signum for THETA and rotate
00702 *
00703                                     THSIGN = -SIGN( ONE, AAPQ )
00704                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
00705                                     T = ONE / ( THETA+THSIGN*
00706      +                                  SQRT( ONE+THETA*THETA ) )
00707                                     CS = SQRT( ONE / ( ONE+T*T ) )
00708                                     SN = T*CS
00709                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
00710                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00711      +                                         ONE+T*APOAQ*AAPQ ) )
00712                                     AAPP = AAPP*SQRT( AMAX1( ZERO, 
00713      +                                         ONE-T*AQOAP*AAPQ ) )
00714 *
00715                                     APOAQ = D( p ) / D( q )
00716                                     AQOAP = D( q ) / D( p )
00717                                     IF( D( p ).GE.ONE ) THEN
00718 *
00719                                        IF( D( q ).GE.ONE ) THEN
00720                                           FASTR( 3 ) = T*APOAQ
00721                                           FASTR( 4 ) = -T*AQOAP
00722                                           D( p ) = D( p )*CS
00723                                           D( q ) = D( q )*CS
00724                                           CALL SROTM( M, A( 1, p ), 1,
00725      +                                                A( 1, q ), 1,
00726      +                                                FASTR )
00727                                           IF( RSVEC )CALL SROTM( MVL,
00728      +                                        V( 1, p ), 1, V( 1, q ),
00729      +                                        1, FASTR )
00730                                        ELSE
00731                                           CALL SAXPY( M, -T*AQOAP,
00732      +                                                A( 1, q ), 1,
00733      +                                                A( 1, p ), 1 )
00734                                           CALL SAXPY( M, CS*SN*APOAQ,
00735      +                                                A( 1, p ), 1,
00736      +                                                A( 1, q ), 1 )
00737                                           IF( RSVEC ) THEN
00738                                              CALL SAXPY( MVL, -T*AQOAP,
00739      +                                                   V( 1, q ), 1,
00740      +                                                   V( 1, p ), 1 )
00741                                              CALL SAXPY( MVL,
00742      +                                                   CS*SN*APOAQ,
00743      +                                                   V( 1, p ), 1,
00744      +                                                   V( 1, q ), 1 )
00745                                           END IF
00746                                           D( p ) = D( p )*CS
00747                                           D( q ) = D( q ) / CS
00748                                        END IF
00749                                     ELSE
00750                                        IF( D( q ).GE.ONE ) THEN
00751                                           CALL SAXPY( M, T*APOAQ,
00752      +                                                A( 1, p ), 1,
00753      +                                                A( 1, q ), 1 )
00754                                           CALL SAXPY( M, -CS*SN*AQOAP,
00755      +                                                A( 1, q ), 1,
00756      +                                                A( 1, p ), 1 )
00757                                           IF( RSVEC ) THEN
00758                                              CALL SAXPY( MVL, T*APOAQ,
00759      +                                                   V( 1, p ), 1,
00760      +                                                   V( 1, q ), 1 )
00761                                              CALL SAXPY( MVL,
00762      +                                                   -CS*SN*AQOAP,
00763      +                                                   V( 1, q ), 1,
00764      +                                                   V( 1, p ), 1 )
00765                                           END IF
00766                                           D( p ) = D( p ) / CS
00767                                           D( q ) = D( q )*CS
00768                                        ELSE
00769                                           IF( D( p ).GE.D( q ) ) THEN
00770                                              CALL SAXPY( M, -T*AQOAP,
00771      +                                                   A( 1, q ), 1,
00772      +                                                   A( 1, p ), 1 )
00773                                              CALL SAXPY( M, CS*SN*APOAQ,
00774      +                                                   A( 1, p ), 1,
00775      +                                                   A( 1, q ), 1 )
00776                                              D( p ) = D( p )*CS
00777                                              D( q ) = D( q ) / CS
00778                                              IF( RSVEC ) THEN
00779                                                 CALL SAXPY( MVL,
00780      +                                               -T*AQOAP,
00781      +                                               V( 1, q ), 1,
00782      +                                               V( 1, p ), 1 )
00783                                                 CALL SAXPY( MVL,
00784      +                                               CS*SN*APOAQ,
00785      +                                               V( 1, p ), 1,
00786      +                                               V( 1, q ), 1 )
00787                                              END IF
00788                                           ELSE
00789                                              CALL SAXPY( M, T*APOAQ,
00790      +                                                   A( 1, p ), 1,
00791      +                                                   A( 1, q ), 1 )
00792                                              CALL SAXPY( M,
00793      +                                                   -CS*SN*AQOAP,
00794      +                                                   A( 1, q ), 1,
00795      +                                                   A( 1, p ), 1 )
00796                                              D( p ) = D( p ) / CS
00797                                              D( q ) = D( q )*CS
00798                                              IF( RSVEC ) THEN
00799                                                 CALL SAXPY( MVL,
00800      +                                               T*APOAQ, V( 1, p ),
00801      +                                               1, V( 1, q ), 1 )
00802                                                 CALL SAXPY( MVL,
00803      +                                               -CS*SN*AQOAP,
00804      +                                               V( 1, q ), 1,
00805      +                                               V( 1, p ), 1 )
00806                                              END IF
00807                                           END IF
00808                                        END IF
00809                                     END IF
00810                                  END IF
00811 *
00812                               ELSE
00813                                  IF( AAPP.GT.AAQQ ) THEN
00814                                     CALL SCOPY( M, A( 1, p ), 1, WORK,
00815      +                                          1 )
00816                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
00817      +                                           M, 1, WORK, LDA, IERR )
00818                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
00819      +                                           M, 1, A( 1, q ), LDA,
00820      +                                           IERR )
00821                                     TEMP1 = -AAPQ*D( p ) / D( q )
00822                                     CALL SAXPY( M, TEMP1, WORK, 1,
00823      +                                          A( 1, q ), 1 )
00824                                     CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
00825      +                                           M, 1, A( 1, q ), LDA,
00826      +                                           IERR )
00827                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00828      +                                         ONE-AAPQ*AAPQ ) )
00829                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
00830                                  ELSE
00831                                     CALL SCOPY( M, A( 1, q ), 1, WORK,
00832      +                                          1 )
00833                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
00834      +                                           M, 1, WORK, LDA, IERR )
00835                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
00836      +                                           M, 1, A( 1, p ), LDA,
00837      +                                           IERR )
00838                                     TEMP1 = -AAPQ*D( q ) / D( p )
00839                                     CALL SAXPY( M, TEMP1, WORK, 1,
00840      +                                          A( 1, p ), 1 )
00841                                     CALL SLASCL( 'G', 0, 0, ONE, AAPP,
00842      +                                           M, 1, A( 1, p ), LDA,
00843      +                                           IERR )
00844                                     SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
00845      +                                         ONE-AAPQ*AAPQ ) )
00846                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
00847                                  END IF
00848                               END IF
00849 *           END IF ROTOK THEN ... ELSE
00850 *
00851 *           In the case of cancellation in updating SVA(q)
00852 *           .. recompute SVA(q)
00853                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
00854      +                            THEN
00855                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
00856      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
00857                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
00858      +                                         D( q )
00859                                  ELSE
00860                                     T = ZERO
00861                                     AAQQ = ONE
00862                                     CALL SLASSQ( M, A( 1, q ), 1, T,
00863      +                                           AAQQ )
00864                                     SVA( q ) = T*SQRT( AAQQ )*D( q )
00865                                  END IF
00866                               END IF
00867                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
00868                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
00869      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
00870                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
00871      +                                     D( p )
00872                                  ELSE
00873                                     T = ZERO
00874                                     AAPP = ONE
00875                                     CALL SLASSQ( M, A( 1, p ), 1, T,
00876      +                                           AAPP )
00877                                     AAPP = T*SQRT( AAPP )*D( p )
00878                                  END IF
00879                                  SVA( p ) = AAPP
00880                               END IF
00881 *              end of OK rotation
00882                            ELSE
00883                               NOTROT = NOTROT + 1
00884                               PSKIPPED = PSKIPPED + 1
00885                               IJBLSK = IJBLSK + 1
00886                            END IF
00887                         ELSE
00888                            NOTROT = NOTROT + 1
00889                            PSKIPPED = PSKIPPED + 1
00890                            IJBLSK = IJBLSK + 1
00891                         END IF
00892 *
00893                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
00894      +                      THEN
00895                            SVA( p ) = AAPP
00896                            NOTROT = 0
00897                            GO TO 2011
00898                         END IF
00899                         IF( ( i.LE.SWBAND ) .AND.
00900      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
00901                            AAPP = -AAPP
00902                            NOTROT = 0
00903                            GO TO 2203
00904                         END IF
00905 *
00906  2200                CONTINUE
00907 *        end of the q-loop
00908  2203                CONTINUE
00909 *
00910                      SVA( p ) = AAPP
00911 *
00912                   ELSE
00913                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
00914      +                   MIN0( jgl+KBL-1, N ) - jgl + 1
00915                      IF( AAPP.LT.ZERO )NOTROT = 0
00916                   END IF
00917 
00918  2100          CONTINUE
00919 *     end of the p-loop
00920  2010       CONTINUE
00921 *     end of the jbc-loop
00922  2011       CONTINUE
00923 *2011 bailed out of the jbc-loop
00924             DO 2012 p = igl, MIN0( igl+KBL-1, N )
00925                SVA( p ) = ABS( SVA( p ) )
00926  2012       CONTINUE
00927 *
00928  2000    CONTINUE
00929 *2000 :: end of the ibr-loop
00930 *
00931 *     .. update SVA(N)
00932          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
00933      +       THEN
00934             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
00935          ELSE
00936             T = ZERO
00937             AAPP = ONE
00938             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
00939             SVA( N ) = T*SQRT( AAPP )*D( N )
00940          END IF
00941 *
00942 *     Additional steering devices
00943 *
00944          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
00945      +       ( ISWROT.LE.N ) ) )SWBAND = i
00946 *
00947          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.FLOAT( N )*TOL ) .AND.
00948      +       ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
00949             GO TO 1994
00950          END IF
00951 *
00952          IF( NOTROT.GE.EMPTSW )GO TO 1994
00953 
00954  1993 CONTINUE
00955 *     end i=1:NSWEEP loop
00956 * #:) Reaching this point means that the procedure has comleted the given
00957 *     number of iterations.
00958       INFO = NSWEEP - 1
00959       GO TO 1995
00960  1994 CONTINUE
00961 * #:) Reaching this point means that during the i-th sweep all pivots were
00962 *     below the given tolerance, causing early exit.
00963 *
00964       INFO = 0
00965 * #:) INFO = 0 confirms successful iterations.
00966  1995 CONTINUE
00967 *
00968 *     Sort the vector D.
00969       DO 5991 p = 1, N - 1
00970          q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
00971          IF( p.NE.q ) THEN
00972             TEMP1 = SVA( p )
00973             SVA( p ) = SVA( q )
00974             SVA( q ) = TEMP1
00975             TEMP1 = D( p )
00976             D( p ) = D( q )
00977             D( q ) = TEMP1
00978             CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
00979             IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
00980          END IF
00981  5991 CONTINUE
00982 *
00983       RETURN
00984 *     ..
00985 *     .. END OF SGSVJ0
00986 *     ..
00987       END