![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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