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