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