LAPACK 3.3.1
Linear Algebra PACKage

dgesvj.f

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.1)                                  --
00005 *
00006 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
00007 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
00008 *  -- April 2011                                                      --
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. 1/DLAMCH('E') > 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 = -4
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)/AAPQ
00828 *
00829                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
00830 *
00831                                     T = HALF / THETA
00832                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
00833                                     FASTR( 4 ) = -T*WORK( q ) /
00834      $                                           WORK( p )
00835                                     CALL DROTM( M, A( 1, p ), 1,
00836      $                                          A( 1, q ), 1, FASTR )
00837                                     IF( RSVEC )CALL DROTM( MVL,
00838      $                                              V( 1, p ), 1,
00839      $                                              V( 1, q ), 1,
00840      $                                              FASTR )
00841                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
00842      $                                         ONE+T*APOAQ*AAPQ ) )
00843                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
00844      $                                     ONE-T*AQOAP*AAPQ ) )
00845                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
00846 *
00847                                  ELSE
00848 *
00849 *                 .. choose correct signum for THETA and rotate
00850 *
00851                                     THSIGN = -DSIGN( ONE, AAPQ )
00852                                     T = ONE / ( THETA+THSIGN*
00853      $                                  DSQRT( ONE+THETA*THETA ) )
00854                                     CS = DSQRT( ONE / ( ONE+T*T ) )
00855                                     SN = T*CS
00856 *
00857                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
00858                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
00859      $                                         ONE+T*APOAQ*AAPQ ) )
00860                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
00861      $                                     ONE-T*AQOAP*AAPQ ) )
00862 *
00863                                     APOAQ = WORK( p ) / WORK( q )
00864                                     AQOAP = WORK( q ) / WORK( p )
00865                                     IF( WORK( p ).GE.ONE ) THEN
00866                                        IF( WORK( q ).GE.ONE ) THEN
00867                                           FASTR( 3 ) = T*APOAQ
00868                                           FASTR( 4 ) = -T*AQOAP
00869                                           WORK( p ) = WORK( p )*CS
00870                                           WORK( q ) = WORK( q )*CS
00871                                           CALL DROTM( M, A( 1, p ), 1,
00872      $                                                A( 1, q ), 1,
00873      $                                                FASTR )
00874                                           IF( RSVEC )CALL DROTM( MVL,
00875      $                                        V( 1, p ), 1, V( 1, q ),
00876      $                                        1, FASTR )
00877                                        ELSE
00878                                           CALL DAXPY( M, -T*AQOAP,
00879      $                                                A( 1, q ), 1,
00880      $                                                A( 1, p ), 1 )
00881                                           CALL DAXPY( M, CS*SN*APOAQ,
00882      $                                                A( 1, p ), 1,
00883      $                                                A( 1, q ), 1 )
00884                                           WORK( p ) = WORK( p )*CS
00885                                           WORK( q ) = WORK( q ) / CS
00886                                           IF( RSVEC ) THEN
00887                                              CALL DAXPY( MVL, -T*AQOAP,
00888      $                                                   V( 1, q ), 1,
00889      $                                                   V( 1, p ), 1 )
00890                                              CALL DAXPY( MVL,
00891      $                                                   CS*SN*APOAQ,
00892      $                                                   V( 1, p ), 1,
00893      $                                                   V( 1, q ), 1 )
00894                                           END IF
00895                                        END IF
00896                                     ELSE
00897                                        IF( WORK( q ).GE.ONE ) THEN
00898                                           CALL DAXPY( M, T*APOAQ,
00899      $                                                A( 1, p ), 1,
00900      $                                                A( 1, q ), 1 )
00901                                           CALL DAXPY( M, -CS*SN*AQOAP,
00902      $                                                A( 1, q ), 1,
00903      $                                                A( 1, p ), 1 )
00904                                           WORK( p ) = WORK( p ) / CS
00905                                           WORK( q ) = WORK( q )*CS
00906                                           IF( RSVEC ) THEN
00907                                              CALL DAXPY( MVL, T*APOAQ,
00908      $                                                   V( 1, p ), 1,
00909      $                                                   V( 1, q ), 1 )
00910                                              CALL DAXPY( MVL,
00911      $                                                   -CS*SN*AQOAP,
00912      $                                                   V( 1, q ), 1,
00913      $                                                   V( 1, p ), 1 )
00914                                           END IF
00915                                        ELSE
00916                                           IF( WORK( p ).GE.WORK( q ) )
00917      $                                        THEN
00918                                              CALL DAXPY( M, -T*AQOAP,
00919      $                                                   A( 1, q ), 1,
00920      $                                                   A( 1, p ), 1 )
00921                                              CALL DAXPY( M, CS*SN*APOAQ,
00922      $                                                   A( 1, p ), 1,
00923      $                                                   A( 1, q ), 1 )
00924                                              WORK( p ) = WORK( p )*CS
00925                                              WORK( q ) = WORK( q ) / CS
00926                                              IF( RSVEC ) THEN
00927                                                 CALL DAXPY( MVL,
00928      $                                               -T*AQOAP,
00929      $                                               V( 1, q ), 1,
00930      $                                               V( 1, p ), 1 )
00931                                                 CALL DAXPY( MVL,
00932      $                                               CS*SN*APOAQ,
00933      $                                               V( 1, p ), 1,
00934      $                                               V( 1, q ), 1 )
00935                                              END IF
00936                                           ELSE
00937                                              CALL DAXPY( M, T*APOAQ,
00938      $                                                   A( 1, p ), 1,
00939      $                                                   A( 1, q ), 1 )
00940                                              CALL DAXPY( M,
00941      $                                                   -CS*SN*AQOAP,
00942      $                                                   A( 1, q ), 1,
00943      $                                                   A( 1, p ), 1 )
00944                                              WORK( p ) = WORK( p ) / CS
00945                                              WORK( q ) = WORK( q )*CS
00946                                              IF( RSVEC ) THEN
00947                                                 CALL DAXPY( MVL,
00948      $                                               T*APOAQ, V( 1, p ),
00949      $                                               1, V( 1, q ), 1 )
00950                                                 CALL DAXPY( MVL,
00951      $                                               -CS*SN*AQOAP,
00952      $                                               V( 1, q ), 1,
00953      $                                               V( 1, p ), 1 )
00954                                              END IF
00955                                           END IF
00956                                        END IF
00957                                     END IF
00958                                  END IF
00959 *
00960                               ELSE
00961 *              .. have to use modified Gram-Schmidt like transformation
00962                                  CALL DCOPY( M, A( 1, p ), 1,
00963      $                                       WORK( N+1 ), 1 )
00964                                  CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
00965      $                                        1, WORK( N+1 ), LDA,
00966      $                                        IERR )
00967                                  CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
00968      $                                        1, A( 1, q ), LDA, IERR )
00969                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )
00970                                  CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,
00971      $                                       A( 1, q ), 1 )
00972                                  CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
00973      $                                        1, A( 1, q ), LDA, IERR )
00974                                  SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
00975      $                                      ONE-AAPQ*AAPQ ) )
00976                                  MXSINJ = DMAX1( MXSINJ, SFMIN )
00977                               END IF
00978 *           END IF ROTOK THEN ... ELSE
00979 *
00980 *           In the case of cancellation in updating SVA(q), SVA(p)
00981 *           recompute SVA(q), SVA(p).
00982 *
00983                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
00984      $                            THEN
00985                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
00986      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
00987                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
00988      $                                         WORK( q )
00989                                  ELSE
00990                                     T = ZERO
00991                                     AAQQ = ONE
00992                                     CALL DLASSQ( M, A( 1, q ), 1, T,
00993      $                                           AAQQ )
00994                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
00995                                  END IF
00996                               END IF
00997                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
00998                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
00999      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
01000                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
01001      $                                     WORK( p )
01002                                  ELSE
01003                                     T = ZERO
01004                                     AAPP = ONE
01005                                     CALL DLASSQ( M, A( 1, p ), 1, T,
01006      $                                           AAPP )
01007                                     AAPP = T*DSQRT( AAPP )*WORK( p )
01008                                  END IF
01009                                  SVA( p ) = AAPP
01010                               END IF
01011 *
01012                            ELSE
01013 *        A(:,p) and A(:,q) already numerically orthogonal
01014                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
01015 *[RTD]      SKIPPED  = SKIPPED  + 1
01016                               PSKIPPED = PSKIPPED + 1
01017                            END IF
01018                         ELSE
01019 *        A(:,q) is zero column
01020                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
01021                            PSKIPPED = PSKIPPED + 1
01022                         END IF
01023 *
01024                         IF( ( i.LE.SWBAND ) .AND.
01025      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
01026                            IF( ir1.EQ.0 )AAPP = -AAPP
01027                            NOTROT = 0
01028                            GO TO 2103
01029                         END IF
01030 *
01031  2002                CONTINUE
01032 *     END q-LOOP
01033 *
01034  2103                CONTINUE
01035 *     bailed out of q-loop
01036 *
01037                      SVA( p ) = AAPP
01038 *
01039                   ELSE
01040                      SVA( p ) = AAPP
01041                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
01042      $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
01043                   END IF
01044 *
01045  2001          CONTINUE
01046 *     end of the p-loop
01047 *     end of doing the block ( ibr, ibr )
01048  1002       CONTINUE
01049 *     end of ir1-loop
01050 *
01051 * ... go to the off diagonal blocks
01052 *
01053             igl = ( ibr-1 )*KBL + 1
01054 *
01055             DO 2010 jbc = ibr + 1, NBL
01056 *
01057                jgl = ( jbc-1 )*KBL + 1
01058 *
01059 *        doing the block at ( ibr, jbc )
01060 *
01061                IJBLSK = 0
01062                DO 2100 p = igl, MIN0( igl+KBL-1, N )
01063 *
01064                   AAPP = SVA( p )
01065                   IF( AAPP.GT.ZERO ) THEN
01066 *
01067                      PSKIPPED = 0
01068 *
01069                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
01070 *
01071                         AAQQ = SVA( q )
01072                         IF( AAQQ.GT.ZERO ) THEN
01073                            AAPP0 = AAPP
01074 *
01075 *     .. M x 2 Jacobi SVD ..
01076 *
01077 *        Safe Gram matrix computation
01078 *
01079                            IF( AAQQ.GE.ONE ) THEN
01080                               IF( AAPP.GE.AAQQ ) THEN
01081                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
01082                               ELSE
01083                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
01084                               END IF
01085                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
01086                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
01087      $                                  q ), 1 )*WORK( p )*WORK( q ) /
01088      $                                  AAQQ ) / AAPP
01089                               ELSE
01090                                  CALL DCOPY( M, A( 1, p ), 1,
01091      $                                       WORK( N+1 ), 1 )
01092                                  CALL DLASCL( 'G', 0, 0, AAPP,
01093      $                                        WORK( p ), M, 1,
01094      $                                        WORK( N+1 ), LDA, IERR )
01095                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
01096      $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
01097                               END IF
01098                            ELSE
01099                               IF( AAPP.GE.AAQQ ) THEN
01100                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
01101                               ELSE
01102                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
01103                               END IF
01104                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
01105                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
01106      $                                  q ), 1 )*WORK( p )*WORK( q ) /
01107      $                                  AAQQ ) / AAPP
01108                               ELSE
01109                                  CALL DCOPY( M, A( 1, q ), 1,
01110      $                                       WORK( N+1 ), 1 )
01111                                  CALL DLASCL( 'G', 0, 0, AAQQ,
01112      $                                        WORK( q ), M, 1,
01113      $                                        WORK( N+1 ), LDA, IERR )
01114                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
01115      $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
01116                               END IF
01117                            END IF
01118 *
01119                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
01120 *
01121 *        TO rotate or NOT to rotate, THAT is the question ...
01122 *
01123                            IF( DABS( AAPQ ).GT.TOL ) THEN
01124                               NOTROT = 0
01125 *[RTD]      ROTATED  = ROTATED + 1
01126                               PSKIPPED = 0
01127                               ISWROT = ISWROT + 1
01128 *
01129                               IF( ROTOK ) THEN
01130 *
01131                                  AQOAP = AAQQ / AAPP
01132                                  APOAQ = AAPP / AAQQ
01133                                  THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
01134                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
01135 *
01136                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
01137                                     T = HALF / THETA
01138                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
01139                                     FASTR( 4 ) = -T*WORK( q ) /
01140      $                                           WORK( p )
01141                                     CALL DROTM( M, A( 1, p ), 1,
01142      $                                          A( 1, q ), 1, FASTR )
01143                                     IF( RSVEC )CALL DROTM( MVL,
01144      $                                              V( 1, p ), 1,
01145      $                                              V( 1, q ), 1,
01146      $                                              FASTR )
01147                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
01148      $                                         ONE+T*APOAQ*AAPQ ) )
01149                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
01150      $                                     ONE-T*AQOAP*AAPQ ) )
01151                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
01152                                  ELSE
01153 *
01154 *                 .. choose correct signum for THETA and rotate
01155 *
01156                                     THSIGN = -DSIGN( ONE, AAPQ )
01157                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
01158                                     T = ONE / ( THETA+THSIGN*
01159      $                                  DSQRT( ONE+THETA*THETA ) )
01160                                     CS = DSQRT( ONE / ( ONE+T*T ) )
01161                                     SN = T*CS
01162                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
01163                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
01164      $                                         ONE+T*APOAQ*AAPQ ) )
01165                                     AAPP = AAPP*DSQRT( DMAX1( ZERO, 
01166      $                                     ONE-T*AQOAP*AAPQ ) )
01167 *
01168                                     APOAQ = WORK( p ) / WORK( q )
01169                                     AQOAP = WORK( q ) / WORK( p )
01170                                     IF( WORK( p ).GE.ONE ) THEN
01171 *
01172                                        IF( WORK( q ).GE.ONE ) THEN
01173                                           FASTR( 3 ) = T*APOAQ
01174                                           FASTR( 4 ) = -T*AQOAP
01175                                           WORK( p ) = WORK( p )*CS
01176                                           WORK( q ) = WORK( q )*CS
01177                                           CALL DROTM( M, A( 1, p ), 1,
01178      $                                                A( 1, q ), 1,
01179      $                                                FASTR )
01180                                           IF( RSVEC )CALL DROTM( MVL,
01181      $                                        V( 1, p ), 1, V( 1, q ),
01182      $                                        1, FASTR )
01183                                        ELSE
01184                                           CALL DAXPY( M, -T*AQOAP,
01185      $                                                A( 1, q ), 1,
01186      $                                                A( 1, p ), 1 )
01187                                           CALL DAXPY( M, CS*SN*APOAQ,
01188      $                                                A( 1, p ), 1,
01189      $                                                A( 1, q ), 1 )
01190                                           IF( RSVEC ) THEN
01191                                              CALL DAXPY( MVL, -T*AQOAP,
01192      $                                                   V( 1, q ), 1,
01193      $                                                   V( 1, p ), 1 )
01194                                              CALL DAXPY( MVL,
01195      $                                                   CS*SN*APOAQ,
01196      $                                                   V( 1, p ), 1,
01197      $                                                   V( 1, q ), 1 )
01198                                           END IF
01199                                           WORK( p ) = WORK( p )*CS
01200                                           WORK( q ) = WORK( q ) / CS
01201                                        END IF
01202                                     ELSE
01203                                        IF( WORK( q ).GE.ONE ) THEN
01204                                           CALL DAXPY( M, T*APOAQ,
01205      $                                                A( 1, p ), 1,
01206      $                                                A( 1, q ), 1 )
01207                                           CALL DAXPY( M, -CS*SN*AQOAP,
01208      $                                                A( 1, q ), 1,
01209      $                                                A( 1, p ), 1 )
01210                                           IF( RSVEC ) THEN
01211                                              CALL DAXPY( MVL, T*APOAQ,
01212      $                                                   V( 1, p ), 1,
01213      $                                                   V( 1, q ), 1 )
01214                                              CALL DAXPY( MVL,
01215      $                                                   -CS*SN*AQOAP,
01216      $                                                   V( 1, q ), 1,
01217      $                                                   V( 1, p ), 1 )
01218                                           END IF
01219                                           WORK( p ) = WORK( p ) / CS
01220                                           WORK( q ) = WORK( q )*CS
01221                                        ELSE
01222                                           IF( WORK( p ).GE.WORK( q ) )
01223      $                                        THEN
01224                                              CALL DAXPY( M, -T*AQOAP,
01225      $                                                   A( 1, q ), 1,
01226      $                                                   A( 1, p ), 1 )
01227                                              CALL DAXPY( M, CS*SN*APOAQ,
01228      $                                                   A( 1, p ), 1,
01229      $                                                   A( 1, q ), 1 )
01230                                              WORK( p ) = WORK( p )*CS
01231                                              WORK( q ) = WORK( q ) / CS
01232                                              IF( RSVEC ) THEN
01233                                                 CALL DAXPY( MVL,
01234      $                                               -T*AQOAP,
01235      $                                               V( 1, q ), 1,
01236      $                                               V( 1, p ), 1 )
01237                                                 CALL DAXPY( MVL,
01238      $                                               CS*SN*APOAQ,
01239      $                                               V( 1, p ), 1,
01240      $                                               V( 1, q ), 1 )
01241                                              END IF
01242                                           ELSE
01243                                              CALL DAXPY( M, T*APOAQ,
01244      $                                                   A( 1, p ), 1,
01245      $                                                   A( 1, q ), 1 )
01246                                              CALL DAXPY( M,
01247      $                                                   -CS*SN*AQOAP,
01248      $                                                   A( 1, q ), 1,
01249      $                                                   A( 1, p ), 1 )
01250                                              WORK( p ) = WORK( p ) / CS
01251                                              WORK( q ) = WORK( q )*CS
01252                                              IF( RSVEC ) THEN
01253                                                 CALL DAXPY( MVL,
01254      $                                               T*APOAQ, V( 1, p ),
01255      $                                               1, V( 1, q ), 1 )
01256                                                 CALL DAXPY( MVL,
01257      $                                               -CS*SN*AQOAP,
01258      $                                               V( 1, q ), 1,
01259      $                                               V( 1, p ), 1 )
01260                                              END IF
01261                                           END IF
01262                                        END IF
01263                                     END IF
01264                                  END IF
01265 *
01266                               ELSE
01267                                  IF( AAPP.GT.AAQQ ) THEN
01268                                     CALL DCOPY( M, A( 1, p ), 1,
01269      $                                          WORK( N+1 ), 1 )
01270                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
01271      $                                           M, 1, WORK( N+1 ), LDA,
01272      $                                           IERR )
01273                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
01274      $                                           M, 1, A( 1, q ), LDA,
01275      $                                           IERR )
01276                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )
01277                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),
01278      $                                          1, A( 1, q ), 1 )
01279                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
01280      $                                           M, 1, A( 1, q ), LDA,
01281      $                                           IERR )
01282                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
01283      $                                         ONE-AAPQ*AAPQ ) )
01284                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
01285                                  ELSE
01286                                     CALL DCOPY( M, A( 1, q ), 1,
01287      $                                          WORK( N+1 ), 1 )
01288                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
01289      $                                           M, 1, WORK( N+1 ), LDA,
01290      $                                           IERR )
01291                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
01292      $                                           M, 1, A( 1, p ), LDA,
01293      $                                           IERR )
01294                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )
01295                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),
01296      $                                          1, A( 1, p ), 1 )
01297                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,
01298      $                                           M, 1, A( 1, p ), LDA,
01299      $                                           IERR )
01300                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
01301      $                                         ONE-AAPQ*AAPQ ) )
01302                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
01303                                  END IF
01304                               END IF
01305 *           END IF ROTOK THEN ... ELSE
01306 *
01307 *           In the case of cancellation in updating SVA(q)
01308 *           .. recompute SVA(q)
01309                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
01310      $                            THEN
01311                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
01312      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
01313                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
01314      $                                         WORK( q )
01315                                  ELSE
01316                                     T = ZERO
01317                                     AAQQ = ONE
01318                                     CALL DLASSQ( M, A( 1, q ), 1, T,
01319      $                                           AAQQ )
01320                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
01321                                  END IF
01322                               END IF
01323                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
01324                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
01325      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
01326                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
01327      $                                     WORK( p )
01328                                  ELSE
01329                                     T = ZERO
01330                                     AAPP = ONE
01331                                     CALL DLASSQ( M, A( 1, p ), 1, T,
01332      $                                           AAPP )
01333                                     AAPP = T*DSQRT( AAPP )*WORK( p )
01334                                  END IF
01335                                  SVA( p ) = AAPP
01336                               END IF
01337 *              end of OK rotation
01338                            ELSE
01339                               NOTROT = NOTROT + 1
01340 *[RTD]      SKIPPED  = SKIPPED  + 1
01341                               PSKIPPED = PSKIPPED + 1
01342                               IJBLSK = IJBLSK + 1
01343                            END IF
01344                         ELSE
01345                            NOTROT = NOTROT + 1
01346                            PSKIPPED = PSKIPPED + 1
01347                            IJBLSK = IJBLSK + 1
01348                         END IF
01349 *
01350                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
01351      $                      THEN
01352                            SVA( p ) = AAPP
01353                            NOTROT = 0
01354                            GO TO 2011
01355                         END IF
01356                         IF( ( i.LE.SWBAND ) .AND.
01357      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
01358                            AAPP = -AAPP
01359                            NOTROT = 0
01360                            GO TO 2203
01361                         END IF
01362 *
01363  2200                CONTINUE
01364 *        end of the q-loop
01365  2203                CONTINUE
01366 *
01367                      SVA( p ) = AAPP
01368 *
01369                   ELSE
01370 *
01371                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
01372      $                   MIN0( jgl+KBL-1, N ) - jgl + 1
01373                      IF( AAPP.LT.ZERO )NOTROT = 0
01374 *
01375                   END IF
01376 *
01377  2100          CONTINUE
01378 *     end of the p-loop
01379  2010       CONTINUE
01380 *     end of the jbc-loop
01381  2011       CONTINUE
01382 *2011 bailed out of the jbc-loop
01383             DO 2012 p = igl, MIN0( igl+KBL-1, N )
01384                SVA( p ) = DABS( SVA( p ) )
01385  2012       CONTINUE
01386 ***
01387  2000    CONTINUE
01388 *2000 :: end of the ibr-loop
01389 *
01390 *     .. update SVA(N)
01391          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
01392      $       THEN
01393             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
01394          ELSE
01395             T = ZERO
01396             AAPP = ONE
01397             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
01398             SVA( N ) = T*DSQRT( AAPP )*WORK( N )
01399          END IF
01400 *
01401 *     Additional steering devices
01402 *
01403          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
01404      $       ( ISWROT.LE.N ) ) )SWBAND = i
01405 *
01406          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
01407      $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
01408             GO TO 1994
01409          END IF
01410 *
01411          IF( NOTROT.GE.EMPTSW )GO TO 1994
01412 *
01413  1993 CONTINUE
01414 *     end i=1:NSWEEP loop
01415 *
01416 * #:( Reaching this point means that the procedure has not converged.
01417       INFO = NSWEEP - 1
01418       GO TO 1995
01419 *
01420  1994 CONTINUE
01421 * #:) Reaching this point means numerical convergence after the i-th
01422 *     sweep.
01423 *
01424       INFO = 0
01425 * #:) INFO = 0 confirms successful iterations.
01426  1995 CONTINUE
01427 *
01428 *     Sort the singular values and find how many are above
01429 *     the underflow threshold.
01430 *
01431       N2 = 0
01432       N4 = 0
01433       DO 5991 p = 1, N - 1
01434          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
01435          IF( p.NE.q ) THEN
01436             TEMP1 = SVA( p )
01437             SVA( p ) = SVA( q )
01438             SVA( q ) = TEMP1
01439             TEMP1 = WORK( p )
01440             WORK( p ) = WORK( q )
01441             WORK( q ) = TEMP1
01442             CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
01443             IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
01444          END IF
01445          IF( SVA( p ).NE.ZERO ) THEN
01446             N4 = N4 + 1
01447             IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
01448          END IF
01449  5991 CONTINUE
01450       IF( SVA( N ).NE.ZERO ) THEN
01451          N4 = N4 + 1
01452          IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
01453       END IF
01454 *
01455 *     Normalize the left singular vectors.
01456 *
01457       IF( LSVEC .OR. UCTOL ) THEN
01458          DO 1998 p = 1, N2
01459             CALL DSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
01460  1998    CONTINUE
01461       END IF
01462 *
01463 *     Scale the product of Jacobi rotations (assemble the fast rotations).
01464 *
01465       IF( RSVEC ) THEN
01466          IF( APPLV ) THEN
01467             DO 2398 p = 1, N
01468                CALL DSCAL( MVL, WORK( p ), V( 1, p ), 1 )
01469  2398       CONTINUE
01470          ELSE
01471             DO 2399 p = 1, N
01472                TEMP1 = ONE / DNRM2( MVL, V( 1, p ), 1 )
01473                CALL DSCAL( MVL, TEMP1, V( 1, p ), 1 )
01474  2399       CONTINUE
01475          END IF
01476       END IF
01477 *
01478 *     Undo scaling, if necessary (and possible).
01479       IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
01480      $    SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
01481      $    ( SFMIN / SKL) ) ) ) THEN
01482          DO 2400 p = 1, N
01483             SVA( p ) = SKL*SVA( p )
01484  2400    CONTINUE
01485          SKL= ONE
01486       END IF
01487 *
01488       WORK( 1 ) = SKL
01489 *     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
01490 *     then some of the singular values may overflow or underflow and
01491 *     the spectrum is given in this factored representation.
01492 *
01493       WORK( 2 ) = DBLE( N4 )
01494 *     N4 is the number of computed nonzero singular values of A.
01495 *
01496       WORK( 3 ) = DBLE( N2 )
01497 *     N2 is the number of singular values of A greater than SFMIN.
01498 *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
01499 *     that may carry some information.
01500 *
01501       WORK( 4 ) = DBLE( i )
01502 *     i is the index of the last sweep before declaring convergence.
01503 *
01504       WORK( 5 ) = MXAAPQ
01505 *     MXAAPQ is the largest absolute value of scaled pivots in the
01506 *     last sweep
01507 *
01508       WORK( 6 ) = MXSINJ
01509 *     MXSINJ is the largest absolute value of the sines of Jacobi angles
01510 *     in the last sweep
01511 *
01512       RETURN
01513 *     ..
01514 *     .. END OF DGESVJ
01515 *     ..
01516       END
 All Files Functions