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