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