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