LAPACK: SRC/dgesvj.f Source File

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