LAPACK: SRC/sgesvj.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
00002      +                   LDV, 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
00022       CHARACTER*1        JOBA, JOBU, JOBV
00023 *     ..
00024 *     .. Array Arguments ..
00025       REAL               A( LDA, * ), SVA( N ), V( LDV, * ),
00026      +                   WORK( LWORK )
00027 *     ..
00028 *
00029 *  Purpose
00030 *  =======
00031 *
00032 *  SGESVJ computes the singular value decomposition (SVD) of a real
00033 *  M-by-N matrix A, where M >= N. The SVD of A is written as
00034 *                                     [++]   [xx]   [x0]   [xx]
00035 *               A = U * SIGMA * V^t,  [++] = [xx] * [ox] * [xx]
00036 *                                     [++]   [xx]
00037 *  where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
00038 *  matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
00039 *  of SIGMA are the singular values of A. The columns of U and V are the
00040 *  left and the right singular vectors of A, respectively.
00041 *
00042 *  Further Details
00043 *  ~~~~~~~~~~~~~~~
00044 *  The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
00045 *  rotations. The rotations are implemented as fast scaled rotations of
00046 *  Anda and Park [1]. In the case of underflow of the Jacobi angle, a
00047 *  modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
00048 *  column interchanges of de Rijk [2]. The relative accuracy of the computed
00049 *  singular values and the accuracy of the computed singular vectors (in
00050 *  angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
00051 *  The condition number that determines the accuracy in the full rank case
00052 *  is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
00053 *  spectral condition number. The best performance of this Jacobi SVD
00054 *  procedure is achieved if used in an  accelerated version of Drmac and
00055 *  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
00056 *  Some tunning parameters (marked with [TP]) are available for the
00057 *  implementer.
00058 *  The computational range for the nonzero singular values is the  machine
00059 *  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
00060 *  denormalized singular values can be computed with the corresponding
00061 *  gradual loss of accurate digits.
00062 *
00063 *  Contributors
00064 *  ~~~~~~~~~~~~
00065 *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
00066 *
00067 *  References
00068 *  ~~~~~~~~~~
00069 * [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
00070 *     SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
00071 * [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
00072 *     singular value decomposition on a vector computer.
00073 *     SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
00074 * [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
00075 * [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
00076 *     value computation in floating point arithmetic.
00077 *     SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
00078 * [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
00079 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
00080 *     LAPACK Working note 169.
00081 * [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
00082 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
00083 *     LAPACK Working note 170.
00084 * [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
00085 *     QSVD, (H,K)-SVD computations.
00086 *     Department of Mathematics, University of Zagreb, 2008.
00087 *
00088 *  Bugs, Examples and Comments
00089 *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
00090 *  Please report all bugs and send interesting test examples and comments to
00091 *  drmac@math.hr. Thank you.
00092 *
00093 *  Arguments
00094 *  =========
00095 *
00096 *  JOBA    (input) CHARACTER* 1
00097 *          Specifies the structure of A.
00098 *          = 'L': The input matrix A is lower triangular;
00099 *          = 'U': The input matrix A is upper triangular;
00100 *          = 'G': The input matrix A is general M-by-N matrix, M >= N.
00101 *
00102 *  JOBU    (input) CHARACTER*1
00103 *          Specifies whether to compute the left singular vectors
00104 *          (columns of U):
00105 *          = 'U': The left singular vectors corresponding to the nonzero
00106 *                 singular values are computed and returned in the leading
00107 *                 columns of A. See more details in the description of A.
00108 *                 The default numerical orthogonality threshold is set to
00109 *                 approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E').
00110 *          = 'C': Analogous to JOBU='U', except that user can control the
00111 *                 level of numerical orthogonality of the computed left
00112 *                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
00113 *                 CTOL is given on input in the array WORK.
00114 *                 No CTOL smaller than ONE is allowed. CTOL greater
00115 *                 than 1 / EPS is meaningless. The option 'C'
00116 *                 can be used if M*EPS is satisfactory orthogonality
00117 *                 of the computed left singular vectors, so CTOL=M could
00118 *                 save few sweeps of Jacobi rotations.
00119 *                 See the descriptions of A and WORK(1).
00120 *          = 'N': The matrix U is not computed. However, see the
00121 *                 description of A.
00122 *
00123 *  JOBV    (input) CHARACTER*1
00124 *          Specifies whether to compute the right singular vectors, that
00125 *          is, the matrix V:
00126 *          = 'V' : the matrix V is computed and returned in the array V
00127 *          = 'A' : the Jacobi rotations are applied to the MV-by-N
00128 *                  array V. In other words, the right singular vector
00129 *                  matrix V is not computed explicitly; instead it is
00130 *                  applied to an MV-by-N matrix initially stored in the
00131 *                  first MV rows of V.
00132 *          = 'N' : the matrix V is not computed and the array V is not
00133 *                  referenced
00134 *
00135 *  M       (input) INTEGER
00136 *          The number of rows of the input matrix A.  M >= 0.
00137 *
00138 *  N       (input) INTEGER
00139 *          The number of columns of the input matrix A.
00140 *          M >= N >= 0.
00141 *
00142 *  A       (input/output) REAL array, dimension (LDA,N)
00143 *          On entry, the M-by-N matrix A.
00144 *          On exit,
00145 *          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C':
00146 *                 If INFO .EQ. 0 :
00147 *                 RANKA orthonormal columns of U are returned in the
00148 *                 leading RANKA columns of the array A. Here RANKA <= N
00149 *                 is the number of computed singular values of A that are
00150 *                 above the underflow threshold SLAMCH('S'). The singular
00151 *                 vectors corresponding to underflowed or zero singular
00152 *                 values are not computed. The value of RANKA is returned
00153 *                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
00154 *                 descriptions of SVA and WORK. The computed columns of U
00155 *                 are mutually numerically orthogonal up to approximately
00156 *                 TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
00157 *                 see the description of JOBU.
00158 *                 If INFO .GT. 0,
00159 *                 the procedure SGESVJ did not converge in the given number
00160 *                 of iterations (sweeps). In that case, the computed
00161 *                 columns of U may not be orthogonal up to TOL. The output
00162 *                 U (stored in A), SIGMA (given by the computed singular
00163 *                 values in SVA(1:N)) and V is still a decomposition of the
00164 *                 input matrix A in the sense that the residual
00165 *                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
00166 *          If JOBU .EQ. 'N':
00167 *                 If INFO .EQ. 0 :
00168 *                 Note that the left singular vectors are 'for free' in the
00169 *                 one-sided Jacobi SVD algorithm. However, if only the
00170 *                 singular values are needed, the level of numerical
00171 *                 orthogonality of U is not an issue and iterations are
00172 *                 stopped when the columns of the iterated matrix are
00173 *                 numerically orthogonal up to approximately M*EPS. Thus,
00174 *                 on exit, A contains the columns of U scaled with the
00175 *                 corresponding singular values.
00176 *                 If INFO .GT. 0 :
00177 *                 the procedure SGESVJ did not converge in the given number
00178 *                 of iterations (sweeps).
00179 *
00180 *  LDA     (input) INTEGER
00181 *          The leading dimension of the array A.  LDA >= max(1,M).
00182 *
00183 *  SVA     (workspace/output) REAL array, dimension (N)
00184 *          On exit,
00185 *          If INFO .EQ. 0 :
00186 *          depending on the value SCALE = WORK(1), we have:
00187 *                 If SCALE .EQ. ONE:
00188 *                 SVA(1:N) contains the computed singular values of A.
00189 *                 During the computation SVA contains the Euclidean column
00190 *                 norms of the iterated matrices in the array A.
00191 *                 If SCALE .NE. ONE:
00192 *                 The singular values of A are SCALE*SVA(1:N), and this
00193 *                 factored representation is due to the fact that some of the
00194 *                 singular values of A might underflow or overflow.
00195 *
00196 *          If INFO .GT. 0 :
00197 *          the procedure SGESVJ did not converge in the given number of
00198 *          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
00199 *
00200 *  MV      (input) INTEGER
00201 *          If JOBV .EQ. 'A', then the product of Jacobi rotations in SGESVJ
00202 *          is applied to the first MV rows of V. See the description of JOBV.
00203 *
00204 *  V       (input/output) REAL array, dimension (LDV,N)
00205 *          If JOBV = 'V', then V contains on exit the N-by-N matrix of
00206 *                         the right singular vectors;
00207 *          If JOBV = 'A', then V contains the product of the computed right
00208 *                         singular vector matrix and the initial matrix in
00209 *                         the array V.
00210 *          If JOBV = 'N', then V is not referenced.
00211 *
00212 *  LDV     (input) INTEGER
00213 *          The leading dimension of the array V, LDV .GE. 1.
00214 *          If JOBV .EQ. 'V', then LDV .GE. max(1,N).
00215 *          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
00216 *
00217 *  WORK    (input/workspace/output) REAL array, dimension max(4,M+N).
00218 *          On entry,
00219 *          If JOBU .EQ. 'C' :
00220 *          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
00221 *                    The process stops if all columns of A are mutually
00222 *                    orthogonal up to CTOL*EPS, EPS=SLAMCH('E').
00223 *                    It is required that CTOL >= ONE, i.e. it is not
00224 *                    allowed to force the routine to obtain orthogonality
00225 *                    below EPSILON.
00226 *          On exit,
00227 *          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
00228 *                    are the computed singular vcalues of A.
00229 *                    (See description of SVA().)
00230 *          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
00231 *                    singular values.
00232 *          WORK(3) = NINT(WORK(3)) is the number of the computed singular
00233 *                    values that are larger than the underflow threshold.
00234 *          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
00235 *                    rotations needed for numerical convergence.
00236 *          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
00237 *                    This is useful information in cases when SGESVJ did
00238 *                    not converge, as it can be used to estimate whether
00239 *                    the output is stil useful and for post festum analysis.
00240 *          WORK(6) = the largest absolute value over all sines of the
00241 *                    Jacobi rotation angles in the last sweep. It can be
00242 *                    useful for a post festum analysis.
00243 *
00244 *  LWORK   length of WORK, WORK >= MAX(6,M+N)
00245 *
00246 *  INFO    (output) INTEGER
00247 *          = 0 : successful exit.
00248 *          < 0 : if INFO = -i, then the i-th argument had an illegal value
00249 *          > 0 : SGESVJ did not converge in the maximal allowed number (30)
00250 *                of sweeps. The output may still be useful. See the
00251 *                description of WORK.
00252 *  =====================================================================
00253 *
00254 *     .. Local Parameters ..
00255       REAL               ZERO, HALF, ONE, TWO
00256       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
00257      +                   TWO = 2.0E0 )
00258       INTEGER            NSWEEP
00259       PARAMETER          ( NSWEEP = 30 )
00260 *     ..
00261 *     .. Local Scalars ..
00262       REAL               AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
00263      +                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
00264      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
00265      +                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
00266      +                   THSIGN, TOL
00267       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
00268      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
00269      +                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
00270      +                   SWBAND
00271       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
00272      +                   RSVEC, UCTOL, UPPER
00273 *     ..
00274 *     .. Local Arrays ..
00275       REAL               FASTR( 5 )
00276 *     ..
00277 *     .. Intrinsic Functions ..
00278       INTRINSIC          ABS, AMAX1, AMIN1, FLOAT, MIN0, SIGN, SQRT
00279 *     ..
00280 *     .. External Functions ..
00281 *     from BLAS
00282       REAL               SDOT, SNRM2
00283       EXTERNAL           SDOT, SNRM2
00284       INTEGER            ISAMAX
00285       EXTERNAL           ISAMAX
00286 *     from LAPACK
00287       REAL               SLAMCH
00288       EXTERNAL           SLAMCH
00289       LOGICAL            LSAME
00290       EXTERNAL           LSAME
00291 *     ..
00292 *     .. External Subroutines ..
00293 *     from BLAS
00294       EXTERNAL           SAXPY, SCOPY, SROTM, SSCAL, SSWAP
00295 *     from LAPACK
00296       EXTERNAL           SLASCL, SLASET, SLASSQ, XERBLA
00297 *
00298       EXTERNAL           SGSVJ0, SGSVJ1
00299 *     ..
00300 *     .. Executable Statements ..
00301 *
00302 *     Test the input arguments
00303 *
00304       LSVEC = LSAME( JOBU, 'U' )
00305       UCTOL = LSAME( JOBU, 'C' )
00306       RSVEC = LSAME( JOBV, 'V' )
00307       APPLV = LSAME( JOBV, 'A' )
00308       UPPER = LSAME( JOBA, 'U' )
00309       LOWER = LSAME( JOBA, 'L' )
00310 *
00311       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
00312          INFO = -1
00313       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
00314          INFO = -2
00315       ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00316          INFO = -3
00317       ELSE IF( M.LT.0 ) THEN
00318          INFO = -4
00319       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
00320          INFO = -5
00321       ELSE IF( LDA.LT.M ) THEN
00322          INFO = -7
00323       ELSE IF( MV.LT.0 ) THEN
00324          INFO = -9
00325       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
00326      +         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
00327          INFO = -11
00328       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
00329          INFO = -12
00330       ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN
00331          INFO = -13
00332       ELSE
00333          INFO = 0
00334       END IF
00335 *
00336 *     #:(
00337       IF( INFO.NE.0 ) THEN
00338          CALL XERBLA( 'SGESVJ', -INFO )
00339          RETURN
00340       END IF
00341 *
00342 * #:) Quick return for void matrix
00343 *
00344       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
00345 *
00346 *     Set numerical parameters
00347 *     The stopping criterion for Jacobi rotations is
00348 *
00349 *     max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
00350 *
00351 *     where EPS is the round-off and CTOL is defined as follows:
00352 *
00353       IF( UCTOL ) THEN
00354 *        ... user controlled
00355          CTOL = WORK( 1 )
00356       ELSE
00357 *        ... default
00358          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
00359             CTOL = SQRT( FLOAT( M ) )
00360          ELSE
00361             CTOL = FLOAT( M )
00362          END IF
00363       END IF
00364 *     ... and the machine dependent parameters are
00365 *[!]  (Make sure that SLAMCH() works properly on the target machine.)
00366 *
00367       EPSLN = SLAMCH( 'Epsilon' )
00368       ROOTEPS = SQRT( EPSLN )
00369       SFMIN = SLAMCH( 'SafeMinimum' )
00370       ROOTSFMIN = SQRT( SFMIN )
00371       SMALL = SFMIN / EPSLN
00372       BIG = SLAMCH( 'Overflow' )
00373       ROOTBIG = ONE / ROOTSFMIN
00374       LARGE = BIG / SQRT( FLOAT( M*N ) )
00375       BIGTHETA = ONE / ROOTEPS
00376 *
00377       TOL = CTOL*EPSLN
00378       ROOTTOL = SQRT( TOL )
00379 *
00380       IF( FLOAT( M )*EPSLN.GE.ONE ) THEN
00381          INFO = -5
00382          CALL XERBLA( 'SGESVJ', -INFO )
00383          RETURN
00384       END IF
00385 *
00386 *     Initialize the right singular vector matrix.
00387 *
00388       IF( RSVEC ) THEN
00389          MVL = N
00390          CALL SLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
00391       ELSE IF( APPLV ) THEN
00392          MVL = MV
00393       END IF
00394       RSVEC = RSVEC .OR. APPLV
00395 *
00396 *     Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
00397 *(!)  If necessary, scale A to protect the largest singular value
00398 *     from overflow. It is possible that saving the largest singular
00399 *     value destroys the information about the small ones.
00400 *     This initial scaling is almost minimal in the sense that the
00401 *     goal is to make sure that no column norm overflows, and that
00402 *     SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
00403 *     in A are detected, the procedure returns with INFO=-6.
00404 *
00405       SKL = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
00406       NOSCALE = .TRUE.
00407       GOSCALE = .TRUE.
00408 *
00409       IF( LOWER ) THEN
00410 *        the input matrix is M-by-N lower triangular (trapezoidal)
00411          DO 1874 p = 1, N
00412             AAPP = ZERO
00413             AAQQ = ONE
00414             CALL SLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
00415             IF( AAPP.GT.BIG ) THEN
00416                INFO = -6
00417                CALL XERBLA( 'SGESVJ', -INFO )
00418                RETURN
00419             END IF
00420             AAQQ = SQRT( AAQQ )
00421             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
00422                SVA( p ) = AAPP*AAQQ
00423             ELSE
00424                NOSCALE = .FALSE.
00425                SVA( p ) = AAPP*( AAQQ*SKL )
00426                IF( GOSCALE ) THEN
00427                   GOSCALE = .FALSE.
00428                   DO 1873 q = 1, p - 1
00429                      SVA( q ) = SVA( q )*SKL
00430  1873             CONTINUE
00431                END IF
00432             END IF
00433  1874    CONTINUE
00434       ELSE IF( UPPER ) THEN
00435 *        the input matrix is M-by-N upper triangular (trapezoidal)
00436          DO 2874 p = 1, N
00437             AAPP = ZERO
00438             AAQQ = ONE
00439             CALL SLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
00440             IF( AAPP.GT.BIG ) THEN
00441                INFO = -6
00442                CALL XERBLA( 'SGESVJ', -INFO )
00443                RETURN
00444             END IF
00445             AAQQ = SQRT( AAQQ )
00446             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
00447                SVA( p ) = AAPP*AAQQ
00448             ELSE
00449                NOSCALE = .FALSE.
00450                SVA( p ) = AAPP*( AAQQ*SKL )
00451                IF( GOSCALE ) THEN
00452                   GOSCALE = .FALSE.
00453                   DO 2873 q = 1, p - 1
00454                      SVA( q ) = SVA( q )*SKL
00455  2873             CONTINUE
00456                END IF
00457             END IF
00458  2874    CONTINUE
00459       ELSE
00460 *        the input matrix is M-by-N general dense
00461          DO 3874 p = 1, N
00462             AAPP = ZERO
00463             AAQQ = ONE
00464             CALL SLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
00465             IF( AAPP.GT.BIG ) THEN
00466                INFO = -6
00467                CALL XERBLA( 'SGESVJ', -INFO )
00468                RETURN
00469             END IF
00470             AAQQ = SQRT( AAQQ )
00471             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
00472                SVA( p ) = AAPP*AAQQ
00473             ELSE
00474                NOSCALE = .FALSE.
00475                SVA( p ) = AAPP*( AAQQ*SKL )
00476                IF( GOSCALE ) THEN
00477                   GOSCALE = .FALSE.
00478                   DO 3873 q = 1, p - 1
00479                      SVA( q ) = SVA( q )*SKL
00480  3873             CONTINUE
00481                END IF
00482             END IF
00483  3874    CONTINUE
00484       END IF
00485 *
00486       IF( NOSCALE )SKL = ONE
00487 *
00488 *     Move the smaller part of the spectrum from the underflow threshold
00489 *(!)  Start by determining the position of the nonzero entries of the
00490 *     array SVA() relative to ( SFMIN, BIG ).
00491 *
00492       AAPP = ZERO
00493       AAQQ = BIG
00494       DO 4781 p = 1, N
00495          IF( SVA( p ).NE.ZERO )AAQQ = AMIN1( AAQQ, SVA( p ) )
00496          AAPP = AMAX1( AAPP, SVA( p ) )
00497  4781 CONTINUE
00498 *
00499 * #:) Quick return for zero matrix
00500 *
00501       IF( AAPP.EQ.ZERO ) THEN
00502          IF( LSVEC )CALL SLASET( 'G', M, N, ZERO, ONE, A, LDA )
00503          WORK( 1 ) = ONE
00504          WORK( 2 ) = ZERO
00505          WORK( 3 ) = ZERO
00506          WORK( 4 ) = ZERO
00507          WORK( 5 ) = ZERO
00508          WORK( 6 ) = ZERO
00509          RETURN
00510       END IF
00511 *
00512 * #:) Quick return for one-column matrix
00513 *
00514       IF( N.EQ.1 ) THEN
00515          IF( LSVEC )CALL SLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
00516      +                           A( 1, 1 ), LDA, IERR )
00517          WORK( 1 ) = ONE / SKL
00518          IF( SVA( 1 ).GE.SFMIN ) THEN
00519             WORK( 2 ) = ONE
00520          ELSE
00521             WORK( 2 ) = ZERO
00522          END IF
00523          WORK( 3 ) = ZERO
00524          WORK( 4 ) = ZERO
00525          WORK( 5 ) = ZERO
00526          WORK( 6 ) = ZERO
00527          RETURN
00528       END IF
00529 *
00530 *     Protect small singular values from underflow, and try to
00531 *     avoid underflows/overflows in computing Jacobi rotations.
00532 *
00533       SN = SQRT( SFMIN / EPSLN )
00534       TEMP1 = SQRT( BIG / FLOAT( N ) )
00535       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
00536      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
00537          TEMP1 = AMIN1( BIG, TEMP1 / AAPP )
00538 *         AAQQ  = AAQQ*TEMP1
00539 *         AAPP  = AAPP*TEMP1
00540       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
00541          TEMP1 = AMIN1( SN / AAQQ, BIG / ( AAPP*SQRT( FLOAT( N ) ) ) )
00542 *         AAQQ  = AAQQ*TEMP1
00543 *         AAPP  = AAPP*TEMP1
00544       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
00545          TEMP1 = AMAX1( SN / AAQQ, TEMP1 / AAPP )
00546 *         AAQQ  = AAQQ*TEMP1
00547 *         AAPP  = AAPP*TEMP1
00548       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
00549          TEMP1 = AMIN1( SN / AAQQ, BIG / ( SQRT( FLOAT( N ) )*AAPP ) )
00550 *         AAQQ  = AAQQ*TEMP1
00551 *         AAPP  = AAPP*TEMP1
00552       ELSE
00553          TEMP1 = ONE
00554       END IF
00555 *
00556 *     Scale, if necessary
00557 *
00558       IF( TEMP1.NE.ONE ) THEN
00559          CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
00560       END IF
00561       SKL = TEMP1*SKL
00562       IF( SKL.NE.ONE ) THEN
00563          CALL SLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
00564          SKL = ONE / SKL
00565       END IF
00566 *
00567 *     Row-cyclic Jacobi SVD algorithm with column pivoting
00568 *
00569       EMPTSW = ( N*( N-1 ) ) / 2
00570       NOTROT = 0
00571       FASTR( 1 ) = ZERO
00572 *
00573 *     A is represented in factored form A = A * diag(WORK), where diag(WORK)
00574 *     is initialized to identity. WORK is updated during fast scaled
00575 *     rotations.
00576 *
00577       DO 1868 q = 1, N
00578          WORK( q ) = ONE
00579  1868 CONTINUE
00580 *
00581 *
00582       SWBAND = 3
00583 *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
00584 *     if SGESVJ is used as a computational routine in the preconditioned
00585 *     Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
00586 *     works on pivots inside a band-like region around the diagonal.
00587 *     The boundaries are determined dynamically, based on the number of
00588 *     pivots above a threshold.
00589 *
00590       KBL = MIN0( 8, N )
00591 *[TP] KBL is a tuning parameter that defines the tile size in the
00592 *     tiling of the p-q loops of pivot pairs. In general, an optimal
00593 *     value of KBL depends on the matrix dimensions and on the
00594 *     parameters of the computer's memory.
00595 *
00596       NBL = N / KBL
00597       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
00598 *
00599       BLSKIP = KBL**2
00600 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
00601 *
00602       ROWSKIP = MIN0( 5, KBL )
00603 *[TP] ROWSKIP is a tuning parameter.
00604 *
00605       LKAHEAD = 1
00606 *[TP] LKAHEAD is a tuning parameter.
00607 *
00608 *     Quasi block transformations, using the lower (upper) triangular
00609 *     structure of the input matrix. The quasi-block-cycling usually
00610 *     invokes cubic convergence. Big part of this cycle is done inside
00611 *     canonical subspaces of dimensions less than M.
00612 *
00613       IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
00614 *[TP] The number of partition levels and the actual partition are
00615 *     tuning parameters.
00616          N4 = N / 4
00617          N2 = N / 2
00618          N34 = 3*N4
00619          IF( APPLV ) THEN
00620             q = 0
00621          ELSE
00622             q = 1
00623          END IF
00624 *
00625          IF( LOWER ) THEN
00626 *
00627 *     This works very well on lower triangular matrices, in particular
00628 *     in the framework of the preconditioned Jacobi SVD (xGEJSV).
00629 *     The idea is simple:
00630 *     [+ 0 0 0]   Note that Jacobi transformations of [0 0]
00631 *     [+ + 0 0]                                       [0 0]
00632 *     [+ + x 0]   actually work on [x 0]              [x 0]
00633 *     [+ + x x]                    [x x].             [x x]
00634 *
00635             CALL SGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
00636      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,
00637      +                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
00638      +                   2, WORK( N+1 ), LWORK-N, IERR )
00639 *
00640             CALL SGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
00641      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
00642      +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
00643      +                   WORK( N+1 ), LWORK-N, IERR )
00644 *
00645             CALL SGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
00646      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
00647      +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
00648      +                   WORK( N+1 ), LWORK-N, IERR )
00649 *
00650             CALL SGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
00651      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,
00652      +                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
00653      +                   WORK( N+1 ), LWORK-N, IERR )
00654 *
00655             CALL SGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
00656      +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
00657      +                   IERR )
00658 *
00659             CALL SGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
00660      +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
00661      +                   LWORK-N, IERR )
00662 *
00663 *
00664          ELSE IF( UPPER ) THEN
00665 *
00666 *
00667             CALL SGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
00668      +                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
00669      +                   IERR )
00670 *
00671             CALL SGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
00672      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
00673      +                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
00674      +                   IERR )
00675 *
00676             CALL SGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
00677      +                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
00678      +                   LWORK-N, IERR )
00679 *
00680             CALL SGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
00681      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,
00682      +                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
00683      +                   WORK( N+1 ), LWORK-N, IERR )
00684 
00685          END IF
00686 *
00687       END IF
00688 *
00689 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
00690 *
00691       DO 1993 i = 1, NSWEEP
00692 *     .. go go go ...
00693 *
00694          MXAAPQ = ZERO
00695          MXSINJ = ZERO
00696          ISWROT = 0
00697 *
00698          NOTROT = 0
00699          PSKIPPED = 0
00700 *
00701 *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
00702 *     1 <= p < q <= N. This is the first step toward a blocked implementation
00703 *     of the rotations. New implementation, based on block transformations,
00704 *     is under development.
00705 *
00706          DO 2000 ibr = 1, NBL
00707 *
00708             igl = ( ibr-1 )*KBL + 1
00709 *
00710             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
00711 *
00712                igl = igl + ir1*KBL
00713 *
00714                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
00715 *
00716 *     .. de Rijk's pivoting
00717 *
00718                   q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
00719                   IF( p.NE.q ) THEN
00720                      CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
00721                      IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1,
00722      +                                      V( 1, q ), 1 )
00723                      TEMP1 = SVA( p )
00724                      SVA( p ) = SVA( q )
00725                      SVA( q ) = TEMP1
00726                      TEMP1 = WORK( p )
00727                      WORK( p ) = WORK( q )
00728                      WORK( q ) = TEMP1
00729                   END IF
00730 *
00731                   IF( ir1.EQ.0 ) THEN
00732 *
00733 *        Column norms are periodically updated by explicit
00734 *        norm computation.
00735 *        Caveat:
00736 *        Unfortunately, some BLAS implementations compute SNRM2(M,A(1,p),1)
00737 *        as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
00738 *        overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
00739 *        underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
00740 *        Hence, SNRM2 cannot be trusted, not even in the case when
00741 *        the true norm is far from the under(over)flow boundaries.
00742 *        If properly implemented SNRM2 is available, the IF-THEN-ELSE
00743 *        below should read "AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)".
00744 *
00745                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
00746      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
00747                         SVA( p ) = SNRM2( M, A( 1, p ), 1 )*WORK( p )
00748                      ELSE
00749                         TEMP1 = ZERO
00750                         AAPP = ONE
00751                         CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
00752                         SVA( p ) = TEMP1*SQRT( AAPP )*WORK( p )
00753                      END IF
00754                      AAPP = SVA( p )
00755                   ELSE
00756                      AAPP = SVA( p )
00757                   END IF
00758 *
00759                   IF( AAPP.GT.ZERO ) THEN
00760 *
00761                      PSKIPPED = 0
00762 *
00763                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
00764 *
00765                         AAQQ = SVA( q )
00766 *
00767                         IF( AAQQ.GT.ZERO ) THEN
00768 *
00769                            AAPP0 = AAPP
00770                            IF( AAQQ.GE.ONE ) THEN
00771                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
00772                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
00773                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00774      +                                  q ), 1 )*WORK( p )*WORK( q ) /
00775      +                                  AAQQ ) / AAPP
00776                               ELSE
00777                                  CALL SCOPY( M, A( 1, p ), 1,
00778      +                                       WORK( N+1 ), 1 )
00779                                  CALL SLASCL( 'G', 0, 0, AAPP,
00780      +                                        WORK( p ), M, 1,
00781      +                                        WORK( N+1 ), LDA, IERR )
00782                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
00783      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
00784                               END IF
00785                            ELSE
00786                               ROTOK = AAPP.LE.( AAQQ / SMALL )
00787                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
00788                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00789      +                                  q ), 1 )*WORK( p )*WORK( q ) /
00790      +                                  AAQQ ) / AAPP
00791                               ELSE
00792                                  CALL SCOPY( M, A( 1, q ), 1,
00793      +                                       WORK( N+1 ), 1 )
00794                                  CALL SLASCL( 'G', 0, 0, AAQQ,
00795      +                                        WORK( q ), M, 1,
00796      +                                        WORK( N+1 ), LDA, IERR )
00797                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
00798      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP
00799                               END IF
00800                            END IF
00801 *
00802                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
00803 *
00804 *        TO rotate or NOT to rotate, THAT is the question ...
00805 *
00806                            IF( ABS( AAPQ ).GT.TOL ) THEN
00807 *
00808 *           .. rotate
00809 *[RTD]      ROTATED = ROTATED + ONE
00810 *
00811                               IF( ir1.EQ.0 ) THEN
00812                                  NOTROT = 0
00813                                  PSKIPPED = 0
00814                                  ISWROT = ISWROT + 1
00815                               END IF
00816 *
00817                               IF( ROTOK ) THEN
00818 *
00819                                  AQOAP = AAQQ / AAPP
00820                                  APOAQ = AAPP / AAQQ
00821                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
00822 *
00823                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
00824 *
00825                                     T = HALF / THETA
00826                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
00827                                     FASTR( 4 ) = -T*WORK( q ) /
00828      +                                           WORK( p )
00829                                     CALL SROTM( M, A( 1, p ), 1,
00830      +                                          A( 1, q ), 1, FASTR )
00831                                     IF( RSVEC )CALL SROTM( MVL,
00832      +                                              V( 1, p ), 1,
00833      +                                              V( 1, q ), 1,
00834      +                                              FASTR )
00835                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00836      +                                         ONE+T*APOAQ*AAPQ ) )
00837                                     AAPP = AAPP*SQRT( AMAX1( ZERO, 
00838      +                                         ONE-T*AQOAP*AAPQ ) )
00839                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
00840 *
00841                                  ELSE
00842 *
00843 *                 .. choose correct signum for THETA and rotate
00844 *
00845                                     THSIGN = -SIGN( ONE, AAPQ )
00846                                     T = ONE / ( THETA+THSIGN*
00847      +                                  SQRT( ONE+THETA*THETA ) )
00848                                     CS = SQRT( ONE / ( ONE+T*T ) )
00849                                     SN = T*CS
00850 *
00851                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
00852                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00853      +                                         ONE+T*APOAQ*AAPQ ) )
00854                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
00855      +                                     ONE-T*AQOAP*AAPQ ) )
00856 *
00857                                     APOAQ = WORK( p ) / WORK( q )
00858                                     AQOAP = WORK( q ) / WORK( p )
00859                                     IF( WORK( p ).GE.ONE ) THEN
00860                                        IF( WORK( q ).GE.ONE ) THEN
00861                                           FASTR( 3 ) = T*APOAQ
00862                                           FASTR( 4 ) = -T*AQOAP
00863                                           WORK( p ) = WORK( p )*CS
00864                                           WORK( q ) = WORK( q )*CS
00865                                           CALL SROTM( M, A( 1, p ), 1,
00866      +                                                A( 1, q ), 1,
00867      +                                                FASTR )
00868                                           IF( RSVEC )CALL SROTM( MVL,
00869      +                                        V( 1, p ), 1, V( 1, q ),
00870      +                                        1, FASTR )
00871                                        ELSE
00872                                           CALL SAXPY( M, -T*AQOAP,
00873      +                                                A( 1, q ), 1,
00874      +                                                A( 1, p ), 1 )
00875                                           CALL SAXPY( M, CS*SN*APOAQ,
00876      +                                                A( 1, p ), 1,
00877      +                                                A( 1, q ), 1 )
00878                                           WORK( p ) = WORK( p )*CS
00879                                           WORK( q ) = WORK( q ) / CS
00880                                           IF( RSVEC ) THEN
00881                                              CALL SAXPY( MVL, -T*AQOAP,
00882      +                                                   V( 1, q ), 1,
00883      +                                                   V( 1, p ), 1 )
00884                                              CALL SAXPY( MVL,
00885      +                                                   CS*SN*APOAQ,
00886      +                                                   V( 1, p ), 1,
00887      +                                                   V( 1, q ), 1 )
00888                                           END IF
00889                                        END IF
00890                                     ELSE
00891                                        IF( WORK( q ).GE.ONE ) THEN
00892                                           CALL SAXPY( M, T*APOAQ,
00893      +                                                A( 1, p ), 1,
00894      +                                                A( 1, q ), 1 )
00895                                           CALL SAXPY( M, -CS*SN*AQOAP,
00896      +                                                A( 1, q ), 1,
00897      +                                                A( 1, p ), 1 )
00898                                           WORK( p ) = WORK( p ) / CS
00899                                           WORK( q ) = WORK( q )*CS
00900                                           IF( RSVEC ) THEN
00901                                              CALL SAXPY( MVL, T*APOAQ,
00902      +                                                   V( 1, p ), 1,
00903      +                                                   V( 1, q ), 1 )
00904                                              CALL SAXPY( MVL,
00905      +                                                   -CS*SN*AQOAP,
00906      +                                                   V( 1, q ), 1,
00907      +                                                   V( 1, p ), 1 )
00908                                           END IF
00909                                        ELSE
00910                                           IF( WORK( p ).GE.WORK( q ) )
00911      +                                        THEN
00912                                              CALL SAXPY( M, -T*AQOAP,
00913      +                                                   A( 1, q ), 1,
00914      +                                                   A( 1, p ), 1 )
00915                                              CALL SAXPY( M, CS*SN*APOAQ,
00916      +                                                   A( 1, p ), 1,
00917      +                                                   A( 1, q ), 1 )
00918                                              WORK( p ) = WORK( p )*CS
00919                                              WORK( q ) = WORK( q ) / CS
00920                                              IF( RSVEC ) THEN
00921                                                 CALL SAXPY( MVL,
00922      +                                               -T*AQOAP,
00923      +                                               V( 1, q ), 1,
00924      +                                               V( 1, p ), 1 )
00925                                                 CALL SAXPY( MVL,
00926      +                                               CS*SN*APOAQ,
00927      +                                               V( 1, p ), 1,
00928      +                                               V( 1, q ), 1 )
00929                                              END IF
00930                                           ELSE
00931                                              CALL SAXPY( M, T*APOAQ,
00932      +                                                   A( 1, p ), 1,
00933      +                                                   A( 1, q ), 1 )
00934                                              CALL SAXPY( M,
00935      +                                                   -CS*SN*AQOAP,
00936      +                                                   A( 1, q ), 1,
00937      +                                                   A( 1, p ), 1 )
00938                                              WORK( p ) = WORK( p ) / CS
00939                                              WORK( q ) = WORK( q )*CS
00940                                              IF( RSVEC ) THEN
00941                                                 CALL SAXPY( MVL,
00942      +                                               T*APOAQ, V( 1, p ),
00943      +                                               1, V( 1, q ), 1 )
00944                                                 CALL SAXPY( MVL,
00945      +                                               -CS*SN*AQOAP,
00946      +                                               V( 1, q ), 1,
00947      +                                               V( 1, p ), 1 )
00948                                              END IF
00949                                           END IF
00950                                        END IF
00951                                     END IF
00952                                  END IF
00953 *
00954                               ELSE
00955 *              .. have to use modified Gram-Schmidt like transformation
00956                                  CALL SCOPY( M, A( 1, p ), 1,
00957      +                                       WORK( N+1 ), 1 )
00958                                  CALL SLASCL( 'G', 0, 0, AAPP, ONE, M,
00959      +                                        1, WORK( N+1 ), LDA,
00960      +                                        IERR )
00961                                  CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M,
00962      +                                        1, A( 1, q ), LDA, IERR )
00963                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )
00964                                  CALL SAXPY( M, TEMP1, WORK( N+1 ), 1,
00965      +                                       A( 1, q ), 1 )
00966                                  CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M,
00967      +                                        1, A( 1, q ), LDA, IERR )
00968                                  SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00969      +                                      ONE-AAPQ*AAPQ ) )
00970                                  MXSINJ = AMAX1( MXSINJ, SFMIN )
00971                               END IF
00972 *           END IF ROTOK THEN ... ELSE
00973 *
00974 *           In the case of cancellation in updating SVA(q), SVA(p)
00975 *           recompute SVA(q), SVA(p).
00976 *
00977                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
00978      +                            THEN
00979                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
00980      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
00981                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
00982      +                                         WORK( q )
00983                                  ELSE
00984                                     T = ZERO
00985                                     AAQQ = ONE
00986                                     CALL SLASSQ( M, A( 1, q ), 1, T,
00987      +                                           AAQQ )
00988                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
00989                                  END IF
00990                               END IF
00991                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
00992                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
00993      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
00994                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
00995      +                                     WORK( p )
00996                                  ELSE
00997                                     T = ZERO
00998                                     AAPP = ONE
00999                                     CALL SLASSQ( M, A( 1, p ), 1, T,
01000      +                                           AAPP )
01001                                     AAPP = T*SQRT( AAPP )*WORK( p )
01002                                  END IF
01003                                  SVA( p ) = AAPP
01004                               END IF
01005 *
01006                            ELSE
01007 *        A(:,p) and A(:,q) already numerically orthogonal
01008                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
01009 *[RTD]      SKIPPED  = SKIPPED  + 1
01010                               PSKIPPED = PSKIPPED + 1
01011                            END IF
01012                         ELSE
01013 *        A(:,q) is zero column
01014                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
01015                            PSKIPPED = PSKIPPED + 1
01016                         END IF
01017 *
01018                         IF( ( i.LE.SWBAND ) .AND.
01019      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
01020                            IF( ir1.EQ.0 )AAPP = -AAPP
01021                            NOTROT = 0
01022                            GO TO 2103
01023                         END IF
01024 *
01025  2002                CONTINUE
01026 *     END q-LOOP
01027 *
01028  2103                CONTINUE
01029 *     bailed out of q-loop
01030 *
01031                      SVA( p ) = AAPP
01032 *
01033                   ELSE
01034                      SVA( p ) = AAPP
01035                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
01036      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
01037                   END IF
01038 *
01039  2001          CONTINUE
01040 *     end of the p-loop
01041 *     end of doing the block ( ibr, ibr )
01042  1002       CONTINUE
01043 *     end of ir1-loop
01044 *
01045 * ... go to the off diagonal blocks
01046 *
01047             igl = ( ibr-1 )*KBL + 1
01048 *
01049             DO 2010 jbc = ibr + 1, NBL
01050 *
01051                jgl = ( jbc-1 )*KBL + 1
01052 *
01053 *        doing the block at ( ibr, jbc )
01054 *
01055                IJBLSK = 0
01056                DO 2100 p = igl, MIN0( igl+KBL-1, N )
01057 *
01058                   AAPP = SVA( p )
01059                   IF( AAPP.GT.ZERO ) THEN
01060 *
01061                      PSKIPPED = 0
01062 *
01063                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
01064 *
01065                         AAQQ = SVA( q )
01066                         IF( AAQQ.GT.ZERO ) THEN
01067                            AAPP0 = AAPP
01068 *
01069 *     .. M x 2 Jacobi SVD ..
01070 *
01071 *        Safe Gram matrix computation
01072 *
01073                            IF( AAQQ.GE.ONE ) THEN
01074                               IF( AAPP.GE.AAQQ ) THEN
01075                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
01076                               ELSE
01077                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
01078                               END IF
01079                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
01080                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
01081      +                                  q ), 1 )*WORK( p )*WORK( q ) /
01082      +                                  AAQQ ) / AAPP
01083                               ELSE
01084                                  CALL SCOPY( M, A( 1, p ), 1,
01085      +                                       WORK( N+1 ), 1 )
01086                                  CALL SLASCL( 'G', 0, 0, AAPP,
01087      +                                        WORK( p ), M, 1,
01088      +                                        WORK( N+1 ), LDA, IERR )
01089                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
01090      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
01091                               END IF
01092                            ELSE
01093                               IF( AAPP.GE.AAQQ ) THEN
01094                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
01095                               ELSE
01096                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
01097                               END IF
01098                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
01099                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
01100      +                                  q ), 1 )*WORK( p )*WORK( q ) /
01101      +                                  AAQQ ) / AAPP
01102                               ELSE
01103                                  CALL SCOPY( M, A( 1, q ), 1,
01104      +                                       WORK( N+1 ), 1 )
01105                                  CALL SLASCL( 'G', 0, 0, AAQQ,
01106      +                                        WORK( q ), M, 1,
01107      +                                        WORK( N+1 ), LDA, IERR )
01108                                  AAPQ = SDOT( M, WORK( N+1 ), 1,
01109      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP
01110                               END IF
01111                            END IF
01112 *
01113                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
01114 *
01115 *        TO rotate or NOT to rotate, THAT is the question ...
01116 *
01117                            IF( ABS( AAPQ ).GT.TOL ) THEN
01118                               NOTROT = 0
01119 *[RTD]      ROTATED  = ROTATED + 1
01120                               PSKIPPED = 0
01121                               ISWROT = ISWROT + 1
01122 *
01123                               IF( ROTOK ) THEN
01124 *
01125                                  AQOAP = AAQQ / AAPP
01126                                  APOAQ = AAPP / AAQQ
01127                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
01128                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
01129 *
01130                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
01131                                     T = HALF / THETA
01132                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
01133                                     FASTR( 4 ) = -T*WORK( q ) /
01134      +                                           WORK( p )
01135                                     CALL SROTM( M, A( 1, p ), 1,
01136      +                                          A( 1, q ), 1, FASTR )
01137                                     IF( RSVEC )CALL SROTM( MVL,
01138      +                                              V( 1, p ), 1,
01139      +                                              V( 1, q ), 1,
01140      +                                              FASTR )
01141                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
01142      +                                         ONE+T*APOAQ*AAPQ ) )
01143                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
01144      +                                     ONE-T*AQOAP*AAPQ ) )
01145                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
01146                                  ELSE
01147 *
01148 *                 .. choose correct signum for THETA and rotate
01149 *
01150                                     THSIGN = -SIGN( ONE, AAPQ )
01151                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
01152                                     T = ONE / ( THETA+THSIGN*
01153      +                                  SQRT( ONE+THETA*THETA ) )
01154                                     CS = SQRT( ONE / ( ONE+T*T ) )
01155                                     SN = T*CS
01156                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
01157                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
01158      +                                         ONE+T*APOAQ*AAPQ ) )
01159                                     AAPP = AAPP*SQRT( AMAX1( ZERO,  
01160      +                                         ONE-T*AQOAP*AAPQ ) )
01161 *
01162                                     APOAQ = WORK( p ) / WORK( q )
01163                                     AQOAP = WORK( q ) / WORK( p )
01164                                     IF( WORK( p ).GE.ONE ) THEN
01165 *
01166                                        IF( WORK( q ).GE.ONE ) THEN
01167                                           FASTR( 3 ) = T*APOAQ
01168                                           FASTR( 4 ) = -T*AQOAP
01169                                           WORK( p ) = WORK( p )*CS
01170                                           WORK( q ) = WORK( q )*CS
01171                                           CALL SROTM( M, A( 1, p ), 1,
01172      +                                                A( 1, q ), 1,
01173      +                                                FASTR )
01174                                           IF( RSVEC )CALL SROTM( MVL,
01175      +                                        V( 1, p ), 1, V( 1, q ),
01176      +                                        1, FASTR )
01177                                        ELSE
01178                                           CALL SAXPY( M, -T*AQOAP,
01179      +                                                A( 1, q ), 1,
01180      +                                                A( 1, p ), 1 )
01181                                           CALL SAXPY( M, CS*SN*APOAQ,
01182      +                                                A( 1, p ), 1,
01183      +                                                A( 1, q ), 1 )
01184                                           IF( RSVEC ) THEN
01185                                              CALL SAXPY( MVL, -T*AQOAP,
01186      +                                                   V( 1, q ), 1,
01187      +                                                   V( 1, p ), 1 )
01188                                              CALL SAXPY( MVL,
01189      +                                                   CS*SN*APOAQ,
01190      +                                                   V( 1, p ), 1,
01191      +                                                   V( 1, q ), 1 )
01192                                           END IF
01193                                           WORK( p ) = WORK( p )*CS
01194                                           WORK( q ) = WORK( q ) / CS
01195                                        END IF
01196                                     ELSE
01197                                        IF( WORK( q ).GE.ONE ) THEN
01198                                           CALL SAXPY( M, T*APOAQ,
01199      +                                                A( 1, p ), 1,
01200      +                                                A( 1, q ), 1 )
01201                                           CALL SAXPY( M, -CS*SN*AQOAP,
01202      +                                                A( 1, q ), 1,
01203      +                                                A( 1, p ), 1 )
01204                                           IF( RSVEC ) THEN
01205                                              CALL SAXPY( MVL, T*APOAQ,
01206      +                                                   V( 1, p ), 1,
01207      +                                                   V( 1, q ), 1 )
01208                                              CALL SAXPY( MVL,
01209      +                                                   -CS*SN*AQOAP,
01210      +                                                   V( 1, q ), 1,
01211      +                                                   V( 1, p ), 1 )
01212                                           END IF
01213                                           WORK( p ) = WORK( p ) / CS
01214                                           WORK( q ) = WORK( q )*CS
01215                                        ELSE
01216                                           IF( WORK( p ).GE.WORK( q ) )
01217      +                                        THEN
01218                                              CALL SAXPY( M, -T*AQOAP,
01219      +                                                   A( 1, q ), 1,
01220      +                                                   A( 1, p ), 1 )
01221                                              CALL SAXPY( M, CS*SN*APOAQ,
01222      +                                                   A( 1, p ), 1,
01223      +                                                   A( 1, q ), 1 )
01224                                              WORK( p ) = WORK( p )*CS
01225                                              WORK( q ) = WORK( q ) / CS
01226                                              IF( RSVEC ) THEN
01227                                                 CALL SAXPY( MVL,
01228      +                                               -T*AQOAP,
01229      +                                               V( 1, q ), 1,
01230      +                                               V( 1, p ), 1 )
01231                                                 CALL SAXPY( MVL,
01232      +                                               CS*SN*APOAQ,
01233      +                                               V( 1, p ), 1,
01234      +                                               V( 1, q ), 1 )
01235                                              END IF
01236                                           ELSE
01237                                              CALL SAXPY( M, T*APOAQ,
01238      +                                                   A( 1, p ), 1,
01239      +                                                   A( 1, q ), 1 )
01240                                              CALL SAXPY( M,
01241      +                                                   -CS*SN*AQOAP,
01242      +                                                   A( 1, q ), 1,
01243      +                                                   A( 1, p ), 1 )
01244                                              WORK( p ) = WORK( p ) / CS
01245                                              WORK( q ) = WORK( q )*CS
01246                                              IF( RSVEC ) THEN
01247                                                 CALL SAXPY( MVL,
01248      +                                               T*APOAQ, V( 1, p ),
01249      +                                               1, V( 1, q ), 1 )
01250                                                 CALL SAXPY( MVL,
01251      +                                               -CS*SN*AQOAP,
01252      +                                               V( 1, q ), 1,
01253      +                                               V( 1, p ), 1 )
01254                                              END IF
01255                                           END IF
01256                                        END IF
01257                                     END IF
01258                                  END IF
01259 *
01260                               ELSE
01261                                  IF( AAPP.GT.AAQQ ) THEN
01262                                     CALL SCOPY( M, A( 1, p ), 1,
01263      +                                          WORK( N+1 ), 1 )
01264                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
01265      +                                           M, 1, WORK( N+1 ), LDA,
01266      +                                           IERR )
01267                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
01268      +                                           M, 1, A( 1, q ), LDA,
01269      +                                           IERR )
01270                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )
01271                                     CALL SAXPY( M, TEMP1, WORK( N+1 ),
01272      +                                          1, A( 1, q ), 1 )
01273                                     CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
01274      +                                           M, 1, A( 1, q ), LDA,
01275      +                                           IERR )
01276                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
01277      +                                         ONE-AAPQ*AAPQ ) )
01278                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
01279                                  ELSE
01280                                     CALL SCOPY( M, A( 1, q ), 1,
01281      +                                          WORK( N+1 ), 1 )
01282                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
01283      +                                           M, 1, WORK( N+1 ), LDA,
01284      +                                           IERR )
01285                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
01286      +                                           M, 1, A( 1, p ), LDA,
01287      +                                           IERR )
01288                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )
01289                                     CALL SAXPY( M, TEMP1, WORK( N+1 ),
01290      +                                          1, A( 1, p ), 1 )
01291                                     CALL SLASCL( 'G', 0, 0, ONE, AAPP,
01292      +                                           M, 1, A( 1, p ), LDA,
01293      +                                           IERR )
01294                                     SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
01295      +                                         ONE-AAPQ*AAPQ ) )
01296                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
01297                                  END IF
01298                               END IF
01299 *           END IF ROTOK THEN ... ELSE
01300 *
01301 *           In the case of cancellation in updating SVA(q)
01302 *           .. recompute SVA(q)
01303                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
01304      +                            THEN
01305                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
01306      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
01307                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
01308      +                                         WORK( q )
01309                                  ELSE
01310                                     T = ZERO
01311                                     AAQQ = ONE
01312                                     CALL SLASSQ( M, A( 1, q ), 1, T,
01313      +                                           AAQQ )
01314                                     SVA( q ) = T*SQRT( AAQQ )*WORK( q )
01315                                  END IF
01316                               END IF
01317                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
01318                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
01319      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
01320                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
01321      +                                     WORK( p )
01322                                  ELSE
01323                                     T = ZERO
01324                                     AAPP = ONE
01325                                     CALL SLASSQ( M, A( 1, p ), 1, T,
01326      +                                           AAPP )
01327                                     AAPP = T*SQRT( AAPP )*WORK( p )
01328                                  END IF
01329                                  SVA( p ) = AAPP
01330                               END IF
01331 *              end of OK rotation
01332                            ELSE
01333                               NOTROT = NOTROT + 1
01334 *[RTD]      SKIPPED  = SKIPPED  + 1
01335                               PSKIPPED = PSKIPPED + 1
01336                               IJBLSK = IJBLSK + 1
01337                            END IF
01338                         ELSE
01339                            NOTROT = NOTROT + 1
01340                            PSKIPPED = PSKIPPED + 1
01341                            IJBLSK = IJBLSK + 1
01342                         END IF
01343 *
01344                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
01345      +                      THEN
01346                            SVA( p ) = AAPP
01347                            NOTROT = 0
01348                            GO TO 2011
01349                         END IF
01350                         IF( ( i.LE.SWBAND ) .AND.
01351      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
01352                            AAPP = -AAPP
01353                            NOTROT = 0
01354                            GO TO 2203
01355                         END IF
01356 *
01357  2200                CONTINUE
01358 *        end of the q-loop
01359  2203                CONTINUE
01360 *
01361                      SVA( p ) = AAPP
01362 *
01363                   ELSE
01364 *
01365                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
01366      +                   MIN0( jgl+KBL-1, N ) - jgl + 1
01367                      IF( AAPP.LT.ZERO )NOTROT = 0
01368 *
01369                   END IF
01370 *
01371  2100          CONTINUE
01372 *     end of the p-loop
01373  2010       CONTINUE
01374 *     end of the jbc-loop
01375  2011       CONTINUE
01376 *2011 bailed out of the jbc-loop
01377             DO 2012 p = igl, MIN0( igl+KBL-1, N )
01378                SVA( p ) = ABS( SVA( p ) )
01379  2012       CONTINUE
01380 ***
01381  2000    CONTINUE
01382 *2000 :: end of the ibr-loop
01383 *
01384 *     .. update SVA(N)
01385          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
01386      +       THEN
01387             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*WORK( N )
01388          ELSE
01389             T = ZERO
01390             AAPP = ONE
01391             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
01392             SVA( N ) = T*SQRT( AAPP )*WORK( N )
01393          END IF
01394 *
01395 *     Additional steering devices
01396 *
01397          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
01398      +       ( ISWROT.LE.N ) ) )SWBAND = i
01399 *
01400          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
01401      +       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
01402             GO TO 1994
01403          END IF
01404 *
01405          IF( NOTROT.GE.EMPTSW )GO TO 1994
01406 *
01407  1993 CONTINUE
01408 *     end i=1:NSWEEP loop
01409 *
01410 * #:( Reaching this point means that the procedure has not converged.
01411       INFO = NSWEEP - 1
01412       GO TO 1995
01413 *
01414  1994 CONTINUE
01415 * #:) Reaching this point means numerical convergence after the i-th
01416 *     sweep.
01417 *
01418       INFO = 0
01419 * #:) INFO = 0 confirms successful iterations.
01420  1995 CONTINUE
01421 *
01422 *     Sort the singular values and find how many are above
01423 *     the underflow threshold.
01424 *
01425       N2 = 0
01426       N4 = 0
01427       DO 5991 p = 1, N - 1
01428          q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
01429          IF( p.NE.q ) THEN
01430             TEMP1 = SVA( p )
01431             SVA( p ) = SVA( q )
01432             SVA( q ) = TEMP1
01433             TEMP1 = WORK( p )
01434             WORK( p ) = WORK( q )
01435             WORK( q ) = TEMP1
01436             CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
01437             IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
01438          END IF
01439          IF( SVA( p ).NE.ZERO ) THEN
01440             N4 = N4 + 1
01441             IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
01442          END IF
01443  5991 CONTINUE
01444       IF( SVA( N ).NE.ZERO ) THEN
01445          N4 = N4 + 1
01446          IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
01447       END IF
01448 *
01449 *     Normalize the left singular vectors.
01450 *
01451       IF( LSVEC .OR. UCTOL ) THEN
01452          DO 1998 p = 1, N2
01453             CALL SSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
01454  1998    CONTINUE
01455       END IF
01456 *
01457 *     Scale the product of Jacobi rotations (assemble the fast rotations).
01458 *
01459       IF( RSVEC ) THEN
01460          IF( APPLV ) THEN
01461             DO 2398 p = 1, N
01462                CALL SSCAL( MVL, WORK( p ), V( 1, p ), 1 )
01463  2398       CONTINUE
01464          ELSE
01465             DO 2399 p = 1, N
01466                TEMP1 = ONE / SNRM2( MVL, V( 1, p ), 1 )
01467                CALL SSCAL( MVL, TEMP1, V( 1, p ), 1 )
01468  2399       CONTINUE
01469          END IF
01470       END IF
01471 *
01472 *     Undo scaling, if necessary (and possible).
01473       IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
01474      +    SKL ) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
01475      +    ( SFMIN / SKL ) ) ) ) THEN
01476          DO 2400 p = 1, N
01477             SVA( p ) = SKL*SVA( p )
01478  2400    CONTINUE
01479          SKL = ONE
01480       END IF
01481 *
01482       WORK( 1 ) = SKL
01483 *     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
01484 *     then some of the singular values may overflow or underflow and
01485 *     the spectrum is given in this factored representation.
01486 *
01487       WORK( 2 ) = FLOAT( N4 )
01488 *     N4 is the number of computed nonzero singular values of A.
01489 *
01490       WORK( 3 ) = FLOAT( N2 )
01491 *     N2 is the number of singular values of A greater than SFMIN.
01492 *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
01493 *     that may carry some information.
01494 *
01495       WORK( 4 ) = FLOAT( i )
01496 *     i is the index of the last sweep before declaring convergence.
01497 *
01498       WORK( 5 ) = MXAAPQ
01499 *     MXAAPQ is the largest absolute value of scaled pivots in the
01500 *     last sweep
01501 *
01502       WORK( 6 ) = MXSINJ
01503 *     MXSINJ is the largest absolute value of the sines of Jacobi angles
01504 *     in the last sweep
01505 *
01506       RETURN
01507 *     ..
01508 *     .. END OF SGESVJ
01509 *     ..
01510       END