|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, 00002 & M, N, A, LDA, SVA, U, LDU, V, LDV, 00003 & WORK, LWORK, IWORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.3.0) -- 00006 * 00007 * -- Contributed by Zlatko Drmac of the University of Zagreb and -- 00008 * -- Kresimir Veselic of the Fernuniversitaet Hagen -- 00009 * November 2010 00010 * 00011 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00012 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00013 * 00014 * This routine is also part of SIGMA (version 1.23, October 23. 2008.) 00015 * SIGMA is a library of algorithms for highly accurate algorithms for 00016 * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the 00017 * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0. 00018 * 00019 * .. Scalar Arguments .. 00020 IMPLICIT NONE 00021 INTEGER INFO, LDA, LDU, LDV, LWORK, M, N 00022 * .. 00023 * .. Array Arguments .. 00024 DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ), 00025 & WORK( LWORK ) 00026 INTEGER IWORK( * ) 00027 CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV 00028 * .. 00029 * 00030 * Purpose 00031 * ======= 00032 * 00033 * DGEJSV computes the singular value decomposition (SVD) of a real M-by-N 00034 * matrix [A], where M >= N. The SVD of [A] is written as 00035 * 00036 * [A] = [U] * [SIGMA] * [V]^t, 00037 * 00038 * where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N 00039 * diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and 00040 * [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are 00041 * the singular values of [A]. The columns of [U] and [V] are the left and 00042 * the right singular vectors of [A], respectively. The matrices [U] and [V] 00043 * are computed and stored in the arrays U and V, respectively. The diagonal 00044 * of [SIGMA] is computed and stored in the array SVA. 00045 * 00046 * Arguments 00047 * ========= 00048 * 00049 * JOBA (input) CHARACTER*1 00050 * Specifies the level of accuracy: 00051 * = 'C': This option works well (high relative accuracy) if A = B * D, 00052 * with well-conditioned B and arbitrary diagonal matrix D. 00053 * The accuracy cannot be spoiled by COLUMN scaling. The 00054 * accuracy of the computed output depends on the condition of 00055 * B, and the procedure aims at the best theoretical accuracy. 00056 * The relative error max_{i=1:N}|d sigma_i| / sigma_i is 00057 * bounded by f(M,N)*epsilon* cond(B), independent of D. 00058 * The input matrix is preprocessed with the QRF with column 00059 * pivoting. This initial preprocessing and preconditioning by 00060 * a rank revealing QR factorization is common for all values of 00061 * JOBA. Additional actions are specified as follows: 00062 * = 'E': Computation as with 'C' with an additional estimate of the 00063 * condition number of B. It provides a realistic error bound. 00064 * = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings 00065 * D1, D2, and well-conditioned matrix C, this option gives 00066 * higher accuracy than the 'C' option. If the structure of the 00067 * input matrix is not known, and relative accuracy is 00068 * desirable, then this option is advisable. The input matrix A 00069 * is preprocessed with QR factorization with FULL (row and 00070 * column) pivoting. 00071 * = 'G' Computation as with 'F' with an additional estimate of the 00072 * condition number of B, where A=D*B. If A has heavily weighted 00073 * rows, then using this condition number gives too pessimistic 00074 * error bound. 00075 * = 'A': Small singular values are the noise and the matrix is treated 00076 * as numerically rank defficient. The error in the computed 00077 * singular values is bounded by f(m,n)*epsilon*||A||. 00078 * The computed SVD A = U * S * V^t restores A up to 00079 * f(m,n)*epsilon*||A||. 00080 * This gives the procedure the licence to discard (set to zero) 00081 * all singular values below N*epsilon*||A||. 00082 * = 'R': Similar as in 'A'. Rank revealing property of the initial 00083 * QR factorization is used do reveal (using triangular factor) 00084 * a gap sigma_{r+1} < epsilon * sigma_r in which case the 00085 * numerical RANK is declared to be r. The SVD is computed with 00086 * absolute error bounds, but more accurately than with 'A'. 00087 * 00088 * JOBU (input) CHARACTER*1 00089 * Specifies whether to compute the columns of U: 00090 * = 'U': N columns of U are returned in the array U. 00091 * = 'F': full set of M left sing. vectors is returned in the array U. 00092 * = 'W': U may be used as workspace of length M*N. See the description 00093 * of U. 00094 * = 'N': U is not computed. 00095 * 00096 * JOBV (input) CHARACTER*1 00097 * Specifies whether to compute the matrix V: 00098 * = 'V': N columns of V are returned in the array V; Jacobi rotations 00099 * are not explicitly accumulated. 00100 * = 'J': N columns of V are returned in the array V, but they are 00101 * computed as the product of Jacobi rotations. This option is 00102 * allowed only if JOBU .NE. 'N', i.e. in computing the full SVD. 00103 * = 'W': V may be used as workspace of length N*N. See the description 00104 * of V. 00105 * = 'N': V is not computed. 00106 * 00107 * JOBR (input) CHARACTER*1 00108 * Specifies the RANGE for the singular values. Issues the licence to 00109 * set to zero small positive singular values if they are outside 00110 * specified range. If A .NE. 0 is scaled so that the largest singular 00111 * value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues 00112 * the licence to kill columns of A whose norm in c*A is less than 00113 * DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN, 00114 * where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E'). 00115 * = 'N': Do not kill small columns of c*A. This option assumes that 00116 * BLAS and QR factorizations and triangular solvers are 00117 * implemented to work in that range. If the condition of A 00118 * is greater than BIG, use DGESVJ. 00119 * = 'R': RESTRICTED range for sigma(c*A) is [DSQRT(SFMIN), DSQRT(BIG)] 00120 * (roughly, as described above). This option is recommended. 00121 * ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 00122 * For computing the singular values in the FULL range [SFMIN,BIG] 00123 * use DGESVJ. 00124 * 00125 * JOBT (input) CHARACTER*1 00126 * If the matrix is square then the procedure may determine to use 00127 * transposed A if A^t seems to be better with respect to convergence. 00128 * If the matrix is not square, JOBT is ignored. This is subject to 00129 * changes in the future. 00130 * The decision is based on two values of entropy over the adjoint 00131 * orbit of A^t * A. See the descriptions of WORK(6) and WORK(7). 00132 * = 'T': transpose if entropy test indicates possibly faster 00133 * convergence of Jacobi process if A^t is taken as input. If A is 00134 * replaced with A^t, then the row pivoting is included automatically. 00135 * = 'N': do not speculate. 00136 * This option can be used to compute only the singular values, or the 00137 * full SVD (U, SIGMA and V). For only one set of singular vectors 00138 * (U or V), the caller should provide both U and V, as one of the 00139 * matrices is used as workspace if the matrix A is transposed. 00140 * The implementer can easily remove this constraint and make the 00141 * code more complicated. See the descriptions of U and V. 00142 * 00143 * JOBP (input) CHARACTER*1 00144 * Issues the licence to introduce structured perturbations to drown 00145 * denormalized numbers. This licence should be active if the 00146 * denormals are poorly implemented, causing slow computation, 00147 * especially in cases of fast convergence (!). For details see [1,2]. 00148 * For the sake of simplicity, this perturbations are included only 00149 * when the full SVD or only the singular values are requested. The 00150 * implementer/user can easily add the perturbation for the cases of 00151 * computing one set of singular vectors. 00152 * = 'P': introduce perturbation 00153 * = 'N': do not perturb 00154 * 00155 * M (input) INTEGER 00156 * The number of rows of the input matrix A. M >= 0. 00157 * 00158 * N (input) INTEGER 00159 * The number of columns of the input matrix A. M >= N >= 0. 00160 * 00161 * A (input/workspace) DOUBLE PRECISION array, dimension (LDA,N) 00162 * On entry, the M-by-N matrix A. 00163 * 00164 * LDA (input) INTEGER 00165 * The leading dimension of the array A. LDA >= max(1,M). 00166 * 00167 * SVA (workspace/output) DOUBLE PRECISION array, dimension (N) 00168 * On exit, 00169 * - For WORK(1)/WORK(2) = ONE: The singular values of A. During the 00170 * computation SVA contains Euclidean column norms of the 00171 * iterated matrices in the array A. 00172 * - For WORK(1) .NE. WORK(2): The singular values of A are 00173 * (WORK(1)/WORK(2)) * SVA(1:N). This factored form is used if 00174 * sigma_max(A) overflows or if small singular values have been 00175 * saved from underflow by scaling the input matrix A. 00176 * - If JOBR='R' then some of the singular values may be returned 00177 * as exact zeros obtained by "set to zero" because they are 00178 * below the numerical rank threshold or are denormalized numbers. 00179 * 00180 * U (workspace/output) DOUBLE PRECISION array, dimension ( LDU, N ) 00181 * If JOBU = 'U', then U contains on exit the M-by-N matrix of 00182 * the left singular vectors. 00183 * If JOBU = 'F', then U contains on exit the M-by-M matrix of 00184 * the left singular vectors, including an ONB 00185 * of the orthogonal complement of the Range(A). 00186 * If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N), 00187 * then U is used as workspace if the procedure 00188 * replaces A with A^t. In that case, [V] is computed 00189 * in U as left singular vectors of A^t and then 00190 * copied back to the V array. This 'W' option is just 00191 * a reminder to the caller that in this case U is 00192 * reserved as workspace of length N*N. 00193 * If JOBU = 'N' U is not referenced. 00194 * 00195 * LDU (input) INTEGER 00196 * The leading dimension of the array U, LDU >= 1. 00197 * IF JOBU = 'U' or 'F' or 'W', then LDU >= M. 00198 * 00199 * V (workspace/output) DOUBLE PRECISION array, dimension ( LDV, N ) 00200 * If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of 00201 * the right singular vectors; 00202 * If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N), 00203 * then V is used as workspace if the pprocedure 00204 * replaces A with A^t. In that case, [U] is computed 00205 * in V as right singular vectors of A^t and then 00206 * copied back to the U array. This 'W' option is just 00207 * a reminder to the caller that in this case V is 00208 * reserved as workspace of length N*N. 00209 * If JOBV = 'N' V is not referenced. 00210 * 00211 * LDV (input) INTEGER 00212 * The leading dimension of the array V, LDV >= 1. 00213 * If JOBV = 'V' or 'J' or 'W', then LDV >= N. 00214 * 00215 * WORK (workspace/output) DOUBLE PRECISION array, dimension at least LWORK. 00216 * On exit, 00217 * WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such 00218 * that SCALE*SVA(1:N) are the computed singular values 00219 * of A. (See the description of SVA().) 00220 * WORK(2) = See the description of WORK(1). 00221 * WORK(3) = SCONDA is an estimate for the condition number of 00222 * column equilibrated A. (If JOBA .EQ. 'E' or 'G') 00223 * SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). 00224 * It is computed using DPOCON. It holds 00225 * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA 00226 * where R is the triangular factor from the QRF of A. 00227 * However, if R is truncated and the numerical rank is 00228 * determined to be strictly smaller than N, SCONDA is 00229 * returned as -1, thus indicating that the smallest 00230 * singular values might be lost. 00231 * 00232 * If full SVD is needed, the following two condition numbers are 00233 * useful for the analysis of the algorithm. They are provied for 00234 * a developer/implementer who is familiar with the details of 00235 * the method. 00236 * 00237 * WORK(4) = an estimate of the scaled condition number of the 00238 * triangular factor in the first QR factorization. 00239 * WORK(5) = an estimate of the scaled condition number of the 00240 * triangular factor in the second QR factorization. 00241 * The following two parameters are computed if JOBT .EQ. 'T'. 00242 * They are provided for a developer/implementer who is familiar 00243 * with the details of the method. 00244 * 00245 * WORK(6) = the entropy of A^t*A :: this is the Shannon entropy 00246 * of diag(A^t*A) / Trace(A^t*A) taken as point in the 00247 * probability simplex. 00248 * WORK(7) = the entropy of A*A^t. 00249 * 00250 * LWORK (input) INTEGER 00251 * Length of WORK to confirm proper allocation of work space. 00252 * LWORK depends on the job: 00253 * 00254 * If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and 00255 * -> .. no scaled condition estimate required ( JOBE.EQ.'N'): 00256 * LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement. 00257 * For optimal performance (blocked code) the optimal value 00258 * is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal 00259 * block size for xGEQP3/xGEQRF. 00260 * -> .. an estimate of the scaled condition number of A is 00261 * required (JOBA='E', 'G'). In this case, LWORK is the maximum 00262 * of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7). 00263 * 00264 * If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), 00265 * -> the minimal requirement is LWORK >= max(2*N+M,7). 00266 * -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), 00267 * where NB is the optimal block size. 00268 * 00269 * If SIGMA and the left singular vectors are needed 00270 * -> the minimal requirement is LWORK >= max(2*N+M,7). 00271 * -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), 00272 * where NB is the optimal block size. 00273 * 00274 * If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and 00275 * -> .. the singular vectors are computed without explicit 00276 * accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N 00277 * -> .. in the iterative part, the Jacobi rotations are 00278 * explicitly accumulated (option, see the description of JOBV), 00279 * then the minimal requirement is LWORK >= max(M+3*N+N*N,7). 00280 * For better performance, if NB is the optimal block size, 00281 * LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7). 00282 * 00283 * IWORK (workspace/output) INTEGER array, dimension M+3*N. 00284 * On exit, 00285 * IWORK(1) = the numerical rank determined after the initial 00286 * QR factorization with pivoting. See the descriptions 00287 * of JOBA and JOBR. 00288 * IWORK(2) = the number of the computed nonzero singular values 00289 * IWORK(3) = if nonzero, a warning message: 00290 * If IWORK(3).EQ.1 then some of the column norms of A 00291 * were denormalized floats. The requested high accuracy 00292 * is not warranted by the data. 00293 * 00294 * INFO (output) INTEGER 00295 * < 0 : if INFO = -i, then the i-th argument had an illegal value. 00296 * = 0 : successfull exit; 00297 * > 0 : DGEJSV did not converge in the maximal allowed number 00298 * of sweeps. The computed values may be inaccurate. 00299 * 00300 * Further Details 00301 * =============== 00302 * 00303 * DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3, 00304 * SGEQRF, and SGELQF as preprocessors and preconditioners. Optionally, an 00305 * additional row pivoting can be used as a preprocessor, which in some 00306 * cases results in much higher accuracy. An example is matrix A with the 00307 * structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned 00308 * diagonal matrices and C is well-conditioned matrix. In that case, complete 00309 * pivoting in the first QR factorizations provides accuracy dependent on the 00310 * condition number of C, and independent of D1, D2. Such higher accuracy is 00311 * not completely understood theoretically, but it works well in practice. 00312 * Further, if A can be written as A = B*D, with well-conditioned B and some 00313 * diagonal D, then the high accuracy is guaranteed, both theoretically and 00314 * in software, independent of D. For more details see [1], [2]. 00315 * The computational range for the singular values can be the full range 00316 * ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS 00317 * & LAPACK routines called by DGEJSV are implemented to work in that range. 00318 * If that is not the case, then the restriction for safe computation with 00319 * the singular values in the range of normalized IEEE numbers is that the 00320 * spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not 00321 * overflow. This code (DGEJSV) is best used in this restricted range, 00322 * meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are 00323 * returned as zeros. See JOBR for details on this. 00324 * Further, this implementation is somewhat slower than the one described 00325 * in [1,2] due to replacement of some non-LAPACK components, and because 00326 * the choice of some tuning parameters in the iterative part (DGESVJ) is 00327 * left to the implementer on a particular machine. 00328 * The rank revealing QR factorization (in this code: SGEQP3) should be 00329 * implemented as in [3]. We have a new version of SGEQP3 under development 00330 * that is more robust than the current one in LAPACK, with a cleaner cut in 00331 * rank defficient cases. It will be available in the SIGMA library [4]. 00332 * If M is much larger than N, it is obvious that the inital QRF with 00333 * column pivoting can be preprocessed by the QRF without pivoting. That 00334 * well known trick is not used in DGEJSV because in some cases heavy row 00335 * weighting can be treated with complete pivoting. The overhead in cases 00336 * M much larger than N is then only due to pivoting, but the benefits in 00337 * terms of accuracy have prevailed. The implementer/user can incorporate 00338 * this extra QRF step easily. The implementer can also improve data movement 00339 * (matrix transpose, matrix copy, matrix transposed copy) - this 00340 * implementation of DGEJSV uses only the simplest, naive data movement. 00341 * 00342 * Contributors 00343 * 00344 * Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) 00345 * 00346 * References 00347 * 00348 * [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. 00349 * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. 00350 * LAPACK Working note 169. 00351 * [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. 00352 * SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. 00353 * LAPACK Working note 170. 00354 * [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR 00355 * factorization software - a case study. 00356 * ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. 00357 * LAPACK Working note 176. 00358 * [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, 00359 * QSVD, (H,K)-SVD computations. 00360 * Department of Mathematics, University of Zagreb, 2008. 00361 * 00362 * Bugs, examples and comments 00363 * 00364 * Please report all bugs and send interesting examples and/or comments to 00365 * drmac@math.hr. Thank you. 00366 * 00367 * ========================================================================== 00368 * 00369 * .. Local Parameters .. 00370 DOUBLE PRECISION ZERO, ONE 00371 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00372 * .. 00373 * .. Local Scalars .. 00374 DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK, 00375 & CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM, 00376 & SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC 00377 INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING 00378 LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC, 00379 & L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, 00380 & NOSCAL, ROWPIV, RSVEC, TRANSP 00381 * .. 00382 * .. Intrinsic Functions .. 00383 INTRINSIC DABS, DLOG, DMAX1, DMIN1, DBLE, 00384 & MAX0, MIN0, IDNINT, DSIGN, DSQRT 00385 * .. 00386 * .. External Functions .. 00387 DOUBLE PRECISION DLAMCH, DNRM2 00388 INTEGER IDAMAX 00389 LOGICAL LSAME 00390 EXTERNAL IDAMAX, LSAME, DLAMCH, DNRM2 00391 * .. 00392 * .. External Subroutines .. 00393 EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL, 00394 & DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ, 00395 & DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA 00396 * 00397 EXTERNAL DGESVJ 00398 * .. 00399 * 00400 * Test the input arguments 00401 * 00402 LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' ) 00403 JRACC = LSAME( JOBV, 'J' ) 00404 RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC 00405 ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' ) 00406 L2RANK = LSAME( JOBA, 'R' ) 00407 L2ABER = LSAME( JOBA, 'A' ) 00408 ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' ) 00409 L2TRAN = LSAME( JOBT, 'T' ) 00410 L2KILL = LSAME( JOBR, 'R' ) 00411 DEFR = LSAME( JOBR, 'N' ) 00412 L2PERT = LSAME( JOBP, 'P' ) 00413 * 00414 IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR. 00415 & ERREST .OR. LSAME( JOBA, 'C' ) )) THEN 00416 INFO = - 1 00417 ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR. 00418 & LSAME( JOBU, 'W' )) ) THEN 00419 INFO = - 2 00420 ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR. 00421 & LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN 00422 INFO = - 3 00423 ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN 00424 INFO = - 4 00425 ELSE IF ( .NOT. ( L2TRAN .OR. LSAME( JOBT, 'N' ) ) ) THEN 00426 INFO = - 5 00427 ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN 00428 INFO = - 6 00429 ELSE IF ( M .LT. 0 ) THEN 00430 INFO = - 7 00431 ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN 00432 INFO = - 8 00433 ELSE IF ( LDA .LT. M ) THEN 00434 INFO = - 10 00435 ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN 00436 INFO = - 13 00437 ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN 00438 INFO = - 14 00439 ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. 00440 & (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR. 00441 & (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND. 00442 & (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR. 00443 & (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. 00444 & (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. 00445 & (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N)) 00446 & .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N))) 00447 & THEN 00448 INFO = - 17 00449 ELSE 00450 * #:) 00451 INFO = 0 00452 END IF 00453 * 00454 IF ( INFO .NE. 0 ) THEN 00455 * #:( 00456 CALL XERBLA( 'DGEJSV', - INFO ) 00457 END IF 00458 * 00459 * Quick return for void matrix (Y3K safe) 00460 * #:) 00461 IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN 00462 * 00463 * Determine whether the matrix U should be M x N or M x M 00464 * 00465 IF ( LSVEC ) THEN 00466 N1 = N 00467 IF ( LSAME( JOBU, 'F' ) ) N1 = M 00468 END IF 00469 * 00470 * Set numerical parameters 00471 * 00472 *! NOTE: Make sure DLAMCH() does not fail on the target architecture. 00473 * 00474 00475 EPSLN = DLAMCH('Epsilon') 00476 SFMIN = DLAMCH('SafeMinimum') 00477 SMALL = SFMIN / EPSLN 00478 BIG = DLAMCH('O') 00479 * BIG = ONE / SFMIN 00480 * 00481 * Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N 00482 * 00483 *(!) If necessary, scale SVA() to protect the largest norm from 00484 * overflow. It is possible that this scaling pushes the smallest 00485 * column norm left from the underflow threshold (extreme case). 00486 * 00487 SCALEM = ONE / DSQRT(DBLE(M)*DBLE(N)) 00488 NOSCAL = .TRUE. 00489 GOSCAL = .TRUE. 00490 DO 1874 p = 1, N 00491 AAPP = ZERO 00492 AAQQ = ONE 00493 CALL DLASSQ( M, A(1,p), 1, AAPP, AAQQ ) 00494 IF ( AAPP .GT. BIG ) THEN 00495 INFO = - 9 00496 CALL XERBLA( 'DGEJSV', -INFO ) 00497 RETURN 00498 END IF 00499 AAQQ = DSQRT(AAQQ) 00500 IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN 00501 SVA(p) = AAPP * AAQQ 00502 ELSE 00503 NOSCAL = .FALSE. 00504 SVA(p) = AAPP * ( AAQQ * SCALEM ) 00505 IF ( GOSCAL ) THEN 00506 GOSCAL = .FALSE. 00507 CALL DSCAL( p-1, SCALEM, SVA, 1 ) 00508 END IF 00509 END IF 00510 1874 CONTINUE 00511 * 00512 IF ( NOSCAL ) SCALEM = ONE 00513 * 00514 AAPP = ZERO 00515 AAQQ = BIG 00516 DO 4781 p = 1, N 00517 AAPP = DMAX1( AAPP, SVA(p) ) 00518 IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) ) 00519 4781 CONTINUE 00520 * 00521 * Quick return for zero M x N matrix 00522 * #:) 00523 IF ( AAPP .EQ. ZERO ) THEN 00524 IF ( LSVEC ) CALL DLASET( 'G', M, N1, ZERO, ONE, U, LDU ) 00525 IF ( RSVEC ) CALL DLASET( 'G', N, N, ZERO, ONE, V, LDV ) 00526 WORK(1) = ONE 00527 WORK(2) = ONE 00528 IF ( ERREST ) WORK(3) = ONE 00529 IF ( LSVEC .AND. RSVEC ) THEN 00530 WORK(4) = ONE 00531 WORK(5) = ONE 00532 END IF 00533 IF ( L2TRAN ) THEN 00534 WORK(6) = ZERO 00535 WORK(7) = ZERO 00536 END IF 00537 IWORK(1) = 0 00538 IWORK(2) = 0 00539 RETURN 00540 END IF 00541 * 00542 * Issue warning if denormalized column norms detected. Override the 00543 * high relative accuracy request. Issue licence to kill columns 00544 * (set them to zero) whose norm is less than sigma_max / BIG (roughly). 00545 * #:( 00546 WARNING = 0 00547 IF ( AAQQ .LE. SFMIN ) THEN 00548 L2RANK = .TRUE. 00549 L2KILL = .TRUE. 00550 WARNING = 1 00551 END IF 00552 * 00553 * Quick return for one-column matrix 00554 * #:) 00555 IF ( N .EQ. 1 ) THEN 00556 * 00557 IF ( LSVEC ) THEN 00558 CALL DLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR ) 00559 CALL DLACPY( 'A', M, 1, A, LDA, U, LDU ) 00560 * computing all M left singular vectors of the M x 1 matrix 00561 IF ( N1 .NE. N ) THEN 00562 CALL DGEQRF( M, N, U,LDU, WORK, WORK(N+1),LWORK-N,IERR ) 00563 CALL DORGQR( M,N1,1, U,LDU,WORK,WORK(N+1),LWORK-N,IERR ) 00564 CALL DCOPY( M, A(1,1), 1, U(1,1), 1 ) 00565 END IF 00566 END IF 00567 IF ( RSVEC ) THEN 00568 V(1,1) = ONE 00569 END IF 00570 IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN 00571 SVA(1) = SVA(1) / SCALEM 00572 SCALEM = ONE 00573 END IF 00574 WORK(1) = ONE / SCALEM 00575 WORK(2) = ONE 00576 IF ( SVA(1) .NE. ZERO ) THEN 00577 IWORK(1) = 1 00578 IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN 00579 IWORK(2) = 1 00580 ELSE 00581 IWORK(2) = 0 00582 END IF 00583 ELSE 00584 IWORK(1) = 0 00585 IWORK(2) = 0 00586 END IF 00587 IF ( ERREST ) WORK(3) = ONE 00588 IF ( LSVEC .AND. RSVEC ) THEN 00589 WORK(4) = ONE 00590 WORK(5) = ONE 00591 END IF 00592 IF ( L2TRAN ) THEN 00593 WORK(6) = ZERO 00594 WORK(7) = ZERO 00595 END IF 00596 RETURN 00597 * 00598 END IF 00599 * 00600 TRANSP = .FALSE. 00601 L2TRAN = L2TRAN .AND. ( M .EQ. N ) 00602 * 00603 AATMAX = -ONE 00604 AATMIN = BIG 00605 IF ( ROWPIV .OR. L2TRAN ) THEN 00606 * 00607 * Compute the row norms, needed to determine row pivoting sequence 00608 * (in the case of heavily row weighted A, row pivoting is strongly 00609 * advised) and to collect information needed to compare the 00610 * structures of A * A^t and A^t * A (in the case L2TRAN.EQ..TRUE.). 00611 * 00612 IF ( L2TRAN ) THEN 00613 DO 1950 p = 1, M 00614 XSC = ZERO 00615 TEMP1 = ONE 00616 CALL DLASSQ( N, A(p,1), LDA, XSC, TEMP1 ) 00617 * DLASSQ gets both the ell_2 and the ell_infinity norm 00618 * in one pass through the vector 00619 WORK(M+N+p) = XSC * SCALEM 00620 WORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1)) 00621 AATMAX = DMAX1( AATMAX, WORK(N+p) ) 00622 IF (WORK(N+p) .NE. ZERO) AATMIN = DMIN1(AATMIN,WORK(N+p)) 00623 1950 CONTINUE 00624 ELSE 00625 DO 1904 p = 1, M 00626 WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) ) 00627 AATMAX = DMAX1( AATMAX, WORK(M+N+p) ) 00628 AATMIN = DMIN1( AATMIN, WORK(M+N+p) ) 00629 1904 CONTINUE 00630 END IF 00631 * 00632 END IF 00633 * 00634 * For square matrix A try to determine whether A^t would be better 00635 * input for the preconditioned Jacobi SVD, with faster convergence. 00636 * The decision is based on an O(N) function of the vector of column 00637 * and row norms of A, based on the Shannon entropy. This should give 00638 * the right choice in most cases when the difference actually matters. 00639 * It may fail and pick the slower converging side. 00640 * 00641 ENTRA = ZERO 00642 ENTRAT = ZERO 00643 IF ( L2TRAN ) THEN 00644 * 00645 XSC = ZERO 00646 TEMP1 = ONE 00647 CALL DLASSQ( N, SVA, 1, XSC, TEMP1 ) 00648 TEMP1 = ONE / TEMP1 00649 * 00650 ENTRA = ZERO 00651 DO 1113 p = 1, N 00652 BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1 00653 IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1) 00654 1113 CONTINUE 00655 ENTRA = - ENTRA / DLOG(DBLE(N)) 00656 * 00657 * Now, SVA().^2/Trace(A^t * A) is a point in the probability simplex. 00658 * It is derived from the diagonal of A^t * A. Do the same with the 00659 * diagonal of A * A^t, compute the entropy of the corresponding 00660 * probability distribution. Note that A * A^t and A^t * A have the 00661 * same trace. 00662 * 00663 ENTRAT = ZERO 00664 DO 1114 p = N+1, N+M 00665 BIG1 = ( ( WORK(p) / XSC )**2 ) * TEMP1 00666 IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1) 00667 1114 CONTINUE 00668 ENTRAT = - ENTRAT / DLOG(DBLE(M)) 00669 * 00670 * Analyze the entropies and decide A or A^t. Smaller entropy 00671 * usually means better input for the algorithm. 00672 * 00673 TRANSP = ( ENTRAT .LT. ENTRA ) 00674 * 00675 * If A^t is better than A, transpose A. 00676 * 00677 IF ( TRANSP ) THEN 00678 * In an optimal implementation, this trivial transpose 00679 * should be replaced with faster transpose. 00680 DO 1115 p = 1, N - 1 00681 DO 1116 q = p + 1, N 00682 TEMP1 = A(q,p) 00683 A(q,p) = A(p,q) 00684 A(p,q) = TEMP1 00685 1116 CONTINUE 00686 1115 CONTINUE 00687 DO 1117 p = 1, N 00688 WORK(M+N+p) = SVA(p) 00689 SVA(p) = WORK(N+p) 00690 1117 CONTINUE 00691 TEMP1 = AAPP 00692 AAPP = AATMAX 00693 AATMAX = TEMP1 00694 TEMP1 = AAQQ 00695 AAQQ = AATMIN 00696 AATMIN = TEMP1 00697 KILL = LSVEC 00698 LSVEC = RSVEC 00699 RSVEC = KILL 00700 IF ( LSVEC ) N1 = N 00701 * 00702 ROWPIV = .TRUE. 00703 END IF 00704 * 00705 END IF 00706 * END IF L2TRAN 00707 * 00708 * Scale the matrix so that its maximal singular value remains less 00709 * than DSQRT(BIG) -- the matrix is scaled so that its maximal column 00710 * has Euclidean norm equal to DSQRT(BIG/N). The only reason to keep 00711 * DSQRT(BIG) instead of BIG is the fact that DGEJSV uses LAPACK and 00712 * BLAS routines that, in some implementations, are not capable of 00713 * working in the full interval [SFMIN,BIG] and that they may provoke 00714 * overflows in the intermediate results. If the singular values spread 00715 * from SFMIN to BIG, then DGESVJ will compute them. So, in that case, 00716 * one should use DGESVJ instead of DGEJSV. 00717 * 00718 BIG1 = DSQRT( BIG ) 00719 TEMP1 = DSQRT( BIG / DBLE(N) ) 00720 * 00721 CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR ) 00722 IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN 00723 AAQQ = ( AAQQ / AAPP ) * TEMP1 00724 ELSE 00725 AAQQ = ( AAQQ * TEMP1 ) / AAPP 00726 END IF 00727 TEMP1 = TEMP1 * SCALEM 00728 CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR ) 00729 * 00730 * To undo scaling at the end of this procedure, multiply the 00731 * computed singular values with USCAL2 / USCAL1. 00732 * 00733 USCAL1 = TEMP1 00734 USCAL2 = AAPP 00735 * 00736 IF ( L2KILL ) THEN 00737 * L2KILL enforces computation of nonzero singular values in 00738 * the restricted range of condition number of the initial A, 00739 * sigma_max(A) / sigma_min(A) approx. DSQRT(BIG)/DSQRT(SFMIN). 00740 XSC = DSQRT( SFMIN ) 00741 ELSE 00742 XSC = SMALL 00743 * 00744 * Now, if the condition number of A is too big, 00745 * sigma_max(A) / sigma_min(A) .GT. DSQRT(BIG/N) * EPSLN / SFMIN, 00746 * as a precaution measure, the full SVD is computed using DGESVJ 00747 * with accumulated Jacobi rotations. This provides numerically 00748 * more robust computation, at the cost of slightly increased run 00749 * time. Depending on the concrete implementation of BLAS and LAPACK 00750 * (i.e. how they behave in presence of extreme ill-conditioning) the 00751 * implementor may decide to remove this switch. 00752 IF ( ( AAQQ.LT.DSQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN 00753 JRACC = .TRUE. 00754 END IF 00755 * 00756 END IF 00757 IF ( AAQQ .LT. XSC ) THEN 00758 DO 700 p = 1, N 00759 IF ( SVA(p) .LT. XSC ) THEN 00760 CALL DLASET( 'A', M, 1, ZERO, ZERO, A(1,p), LDA ) 00761 SVA(p) = ZERO 00762 END IF 00763 700 CONTINUE 00764 END IF 00765 * 00766 * Preconditioning using QR factorization with pivoting 00767 * 00768 IF ( ROWPIV ) THEN 00769 * Optional row permutation (Bjoerck row pivoting): 00770 * A result by Cox and Higham shows that the Bjoerck's 00771 * row pivoting combined with standard column pivoting 00772 * has similar effect as Powell-Reid complete pivoting. 00773 * The ell-infinity norms of A are made nonincreasing. 00774 DO 1952 p = 1, M - 1 00775 q = IDAMAX( M-p+1, WORK(M+N+p), 1 ) + p - 1 00776 IWORK(2*N+p) = q 00777 IF ( p .NE. q ) THEN 00778 TEMP1 = WORK(M+N+p) 00779 WORK(M+N+p) = WORK(M+N+q) 00780 WORK(M+N+q) = TEMP1 00781 END IF 00782 1952 CONTINUE 00783 CALL DLASWP( N, A, LDA, 1, M-1, IWORK(2*N+1), 1 ) 00784 END IF 00785 * 00786 * End of the preparation phase (scaling, optional sorting and 00787 * transposing, optional flushing of small columns). 00788 * 00789 * Preconditioning 00790 * 00791 * If the full SVD is needed, the right singular vectors are computed 00792 * from a matrix equation, and for that we need theoretical analysis 00793 * of the Businger-Golub pivoting. So we use DGEQP3 as the first RR QRF. 00794 * In all other cases the first RR QRF can be chosen by other criteria 00795 * (eg speed by replacing global with restricted window pivoting, such 00796 * as in SGEQPX from TOMS # 782). Good results will be obtained using 00797 * SGEQPX with properly (!) chosen numerical parameters. 00798 * Any improvement of DGEQP3 improves overal performance of DGEJSV. 00799 * 00800 * A * P1 = Q1 * [ R1^t 0]^t: 00801 DO 1963 p = 1, N 00802 * .. all columns are free columns 00803 IWORK(p) = 0 00804 1963 CONTINUE 00805 CALL DGEQP3( M,N,A,LDA, IWORK,WORK, WORK(N+1),LWORK-N, IERR ) 00806 * 00807 * The upper triangular matrix R1 from the first QRF is inspected for 00808 * rank deficiency and possibilities for deflation, or possible 00809 * ill-conditioning. Depending on the user specified flag L2RANK, 00810 * the procedure explores possibilities to reduce the numerical 00811 * rank by inspecting the computed upper triangular factor. If 00812 * L2RANK or L2ABER are up, then DGEJSV will compute the SVD of 00813 * A + dA, where ||dA|| <= f(M,N)*EPSLN. 00814 * 00815 NR = 1 00816 IF ( L2ABER ) THEN 00817 * Standard absolute error bound suffices. All sigma_i with 00818 * sigma_i < N*EPSLN*||A|| are flushed to zero. This is an 00819 * agressive enforcement of lower numerical rank by introducing a 00820 * backward error of the order of N*EPSLN*||A||. 00821 TEMP1 = DSQRT(DBLE(N))*EPSLN 00822 DO 3001 p = 2, N 00823 IF ( DABS(A(p,p)) .GE. (TEMP1*DABS(A(1,1))) ) THEN 00824 NR = NR + 1 00825 ELSE 00826 GO TO 3002 00827 END IF 00828 3001 CONTINUE 00829 3002 CONTINUE 00830 ELSE IF ( L2RANK ) THEN 00831 * .. similarly as above, only slightly more gentle (less agressive). 00832 * Sudden drop on the diagonal of R1 is used as the criterion for 00833 * close-to-rank-defficient. 00834 TEMP1 = DSQRT(SFMIN) 00835 DO 3401 p = 2, N 00836 IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR. 00837 & ( DABS(A(p,p)) .LT. SMALL ) .OR. 00838 & ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402 00839 NR = NR + 1 00840 3401 CONTINUE 00841 3402 CONTINUE 00842 * 00843 ELSE 00844 * The goal is high relative accuracy. However, if the matrix 00845 * has high scaled condition number the relative accuracy is in 00846 * general not feasible. Later on, a condition number estimator 00847 * will be deployed to estimate the scaled condition number. 00848 * Here we just remove the underflowed part of the triangular 00849 * factor. This prevents the situation in which the code is 00850 * working hard to get the accuracy not warranted by the data. 00851 TEMP1 = DSQRT(SFMIN) 00852 DO 3301 p = 2, N 00853 IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR. 00854 & ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302 00855 NR = NR + 1 00856 3301 CONTINUE 00857 3302 CONTINUE 00858 * 00859 END IF 00860 * 00861 ALMORT = .FALSE. 00862 IF ( NR .EQ. N ) THEN 00863 MAXPRJ = ONE 00864 DO 3051 p = 2, N 00865 TEMP1 = DABS(A(p,p)) / SVA(IWORK(p)) 00866 MAXPRJ = DMIN1( MAXPRJ, TEMP1 ) 00867 3051 CONTINUE 00868 IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. 00869 END IF 00870 * 00871 * 00872 SCONDA = - ONE 00873 CONDR1 = - ONE 00874 CONDR2 = - ONE 00875 * 00876 IF ( ERREST ) THEN 00877 IF ( N .EQ. NR ) THEN 00878 IF ( RSVEC ) THEN 00879 * .. V is available as workspace 00880 CALL DLACPY( 'U', N, N, A, LDA, V, LDV ) 00881 DO 3053 p = 1, N 00882 TEMP1 = SVA(IWORK(p)) 00883 CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 ) 00884 3053 CONTINUE 00885 CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1, 00886 & WORK(N+1), IWORK(2*N+M+1), IERR ) 00887 ELSE IF ( LSVEC ) THEN 00888 * .. U is available as workspace 00889 CALL DLACPY( 'U', N, N, A, LDA, U, LDU ) 00890 DO 3054 p = 1, N 00891 TEMP1 = SVA(IWORK(p)) 00892 CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 ) 00893 3054 CONTINUE 00894 CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1, 00895 & WORK(N+1), IWORK(2*N+M+1), IERR ) 00896 ELSE 00897 CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N ) 00898 DO 3052 p = 1, N 00899 TEMP1 = SVA(IWORK(p)) 00900 CALL DSCAL( p, ONE/TEMP1, WORK(N+(p-1)*N+1), 1 ) 00901 3052 CONTINUE 00902 * .. the columns of R are scaled to have unit Euclidean lengths. 00903 CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1, 00904 & WORK(N+N*N+1), IWORK(2*N+M+1), IERR ) 00905 END IF 00906 SCONDA = ONE / DSQRT(TEMP1) 00907 * SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). 00908 * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA 00909 ELSE 00910 SCONDA = - ONE 00911 END IF 00912 END IF 00913 * 00914 L2PERT = L2PERT .AND. ( DABS( A(1,1)/A(NR,NR) ) .GT. DSQRT(BIG1) ) 00915 * If there is no violent scaling, artificial perturbation is not needed. 00916 * 00917 * Phase 3: 00918 * 00919 00920 IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN 00921 * 00922 * Singular Values only 00923 * 00924 * .. transpose A(1:NR,1:N) 00925 DO 1946 p = 1, MIN0( N-1, NR ) 00926 CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 ) 00927 1946 CONTINUE 00928 * 00929 * The following two DO-loops introduce small relative perturbation 00930 * into the strict upper triangle of the lower triangular matrix. 00931 * Small entries below the main diagonal are also changed. 00932 * This modification is useful if the computing environment does not 00933 * provide/allow FLUSH TO ZERO underflow, for it prevents many 00934 * annoying denormalized numbers in case of strongly scaled matrices. 00935 * The perturbation is structured so that it does not introduce any 00936 * new perturbation of the singular values, and it does not destroy 00937 * the job done by the preconditioner. 00938 * The licence for this perturbation is in the variable L2PERT, which 00939 * should be .FALSE. if FLUSH TO ZERO underflow is active. 00940 * 00941 IF ( .NOT. ALMORT ) THEN 00942 * 00943 IF ( L2PERT ) THEN 00944 * XSC = DSQRT(SMALL) 00945 XSC = EPSLN / DBLE(N) 00946 DO 4947 q = 1, NR 00947 TEMP1 = XSC*DABS(A(q,q)) 00948 DO 4949 p = 1, N 00949 IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) 00950 & .OR. ( p .LT. q ) ) 00951 & A(p,q) = DSIGN( TEMP1, A(p,q) ) 00952 4949 CONTINUE 00953 4947 CONTINUE 00954 ELSE 00955 CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, A(1,2),LDA ) 00956 END IF 00957 * 00958 * .. second preconditioning using the QR factorization 00959 * 00960 CALL DGEQRF( N,NR, A,LDA, WORK, WORK(N+1),LWORK-N, IERR ) 00961 * 00962 * .. and transpose upper to lower triangular 00963 DO 1948 p = 1, NR - 1 00964 CALL DCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 ) 00965 1948 CONTINUE 00966 * 00967 END IF 00968 * 00969 * Row-cyclic Jacobi SVD algorithm with column pivoting 00970 * 00971 * .. again some perturbation (a "background noise") is added 00972 * to drown denormals 00973 IF ( L2PERT ) THEN 00974 * XSC = DSQRT(SMALL) 00975 XSC = EPSLN / DBLE(N) 00976 DO 1947 q = 1, NR 00977 TEMP1 = XSC*DABS(A(q,q)) 00978 DO 1949 p = 1, NR 00979 IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) 00980 & .OR. ( p .LT. q ) ) 00981 & A(p,q) = DSIGN( TEMP1, A(p,q) ) 00982 1949 CONTINUE 00983 1947 CONTINUE 00984 ELSE 00985 CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, A(1,2), LDA ) 00986 END IF 00987 * 00988 * .. and one-sided Jacobi rotations are started on a lower 00989 * triangular matrix (plus perturbation which is ignored in 00990 * the part which destroys triangular form (confusing?!)) 00991 * 00992 CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA, 00993 & N, V, LDV, WORK, LWORK, INFO ) 00994 * 00995 SCALEM = WORK(1) 00996 NUMRANK = IDNINT(WORK(2)) 00997 * 00998 * 00999 ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN 01000 * 01001 * -> Singular Values and Right Singular Vectors <- 01002 * 01003 IF ( ALMORT ) THEN 01004 * 01005 * .. in this case NR equals N 01006 DO 1998 p = 1, NR 01007 CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) 01008 1998 CONTINUE 01009 CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) 01010 * 01011 CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA, 01012 & WORK, LWORK, INFO ) 01013 SCALEM = WORK(1) 01014 NUMRANK = IDNINT(WORK(2)) 01015 01016 ELSE 01017 * 01018 * .. two more QR factorizations ( one QRF is not enough, two require 01019 * accumulated product of Jacobi rotations, three are perfect ) 01020 * 01021 CALL DLASET( 'Lower', NR-1, NR-1, ZERO, ZERO, A(2,1), LDA ) 01022 CALL DGELQF( NR, N, A, LDA, WORK, WORK(N+1), LWORK-N, IERR) 01023 CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV ) 01024 CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) 01025 CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1), 01026 & LWORK-2*N, IERR ) 01027 DO 8998 p = 1, NR 01028 CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) 01029 8998 CONTINUE 01030 CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) 01031 * 01032 CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U, 01033 & LDU, WORK(N+1), LWORK, INFO ) 01034 SCALEM = WORK(N+1) 01035 NUMRANK = IDNINT(WORK(N+2)) 01036 IF ( NR .LT. N ) THEN 01037 CALL DLASET( 'A',N-NR, NR, ZERO,ZERO, V(NR+1,1), LDV ) 01038 CALL DLASET( 'A',NR, N-NR, ZERO,ZERO, V(1,NR+1), LDV ) 01039 CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE, V(NR+1,NR+1), LDV ) 01040 END IF 01041 * 01042 CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK, 01043 & V, LDV, WORK(N+1), LWORK-N, IERR ) 01044 * 01045 END IF 01046 * 01047 DO 8991 p = 1, N 01048 CALL DCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA ) 01049 8991 CONTINUE 01050 CALL DLACPY( 'All', N, N, A, LDA, V, LDV ) 01051 * 01052 IF ( TRANSP ) THEN 01053 CALL DLACPY( 'All', N, N, V, LDV, U, LDU ) 01054 END IF 01055 * 01056 ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN 01057 * 01058 * .. Singular Values and Left Singular Vectors .. 01059 * 01060 * .. second preconditioning step to avoid need to accumulate 01061 * Jacobi rotations in the Jacobi iterations. 01062 DO 1965 p = 1, NR 01063 CALL DCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 ) 01064 1965 CONTINUE 01065 CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) 01066 * 01067 CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1), 01068 & LWORK-2*N, IERR ) 01069 * 01070 DO 1967 p = 1, NR - 1 01071 CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 ) 01072 1967 CONTINUE 01073 CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) 01074 * 01075 CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A, 01076 & LDA, WORK(N+1), LWORK-N, INFO ) 01077 SCALEM = WORK(N+1) 01078 NUMRANK = IDNINT(WORK(N+2)) 01079 * 01080 IF ( NR .LT. M ) THEN 01081 CALL DLASET( 'A', M-NR, NR,ZERO, ZERO, U(NR+1,1), LDU ) 01082 IF ( NR .LT. N1 ) THEN 01083 CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1), LDU ) 01084 CALL DLASET( 'A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1), LDU ) 01085 END IF 01086 END IF 01087 * 01088 CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, 01089 & LDU, WORK(N+1), LWORK-N, IERR ) 01090 * 01091 IF ( ROWPIV ) 01092 & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) 01093 * 01094 DO 1974 p = 1, N1 01095 XSC = ONE / DNRM2( M, U(1,p), 1 ) 01096 CALL DSCAL( M, XSC, U(1,p), 1 ) 01097 1974 CONTINUE 01098 * 01099 IF ( TRANSP ) THEN 01100 CALL DLACPY( 'All', N, N, U, LDU, V, LDV ) 01101 END IF 01102 * 01103 ELSE 01104 * 01105 * .. Full SVD .. 01106 * 01107 IF ( .NOT. JRACC ) THEN 01108 * 01109 IF ( .NOT. ALMORT ) THEN 01110 * 01111 * Second Preconditioning Step (QRF [with pivoting]) 01112 * Note that the composition of TRANSPOSE, QRF and TRANSPOSE is 01113 * equivalent to an LQF CALL. Since in many libraries the QRF 01114 * seems to be better optimized than the LQF, we do explicit 01115 * transpose and use the QRF. This is subject to changes in an 01116 * optimized implementation of DGEJSV. 01117 * 01118 DO 1968 p = 1, NR 01119 CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) 01120 1968 CONTINUE 01121 * 01122 * .. the following two loops perturb small entries to avoid 01123 * denormals in the second QR factorization, where they are 01124 * as good as zeros. This is done to avoid painfully slow 01125 * computation with denormals. The relative size of the perturbation 01126 * is a parameter that can be changed by the implementer. 01127 * This perturbation device will be obsolete on machines with 01128 * properly implemented arithmetic. 01129 * To switch it off, set L2PERT=.FALSE. To remove it from the 01130 * code, remove the action under L2PERT=.TRUE., leave the ELSE part. 01131 * The following two loops should be blocked and fused with the 01132 * transposed copy above. 01133 * 01134 IF ( L2PERT ) THEN 01135 XSC = DSQRT(SMALL) 01136 DO 2969 q = 1, NR 01137 TEMP1 = XSC*DABS( V(q,q) ) 01138 DO 2968 p = 1, N 01139 IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) 01140 & .OR. ( p .LT. q ) ) 01141 & V(p,q) = DSIGN( TEMP1, V(p,q) ) 01142 IF ( p. LT. q ) V(p,q) = - V(p,q) 01143 2968 CONTINUE 01144 2969 CONTINUE 01145 ELSE 01146 CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) 01147 END IF 01148 * 01149 * Estimate the row scaled condition number of R1 01150 * (If R1 is rectangular, N > NR, then the condition number 01151 * of the leading NR x NR submatrix is estimated.) 01152 * 01153 CALL DLACPY( 'L', NR, NR, V, LDV, WORK(2*N+1), NR ) 01154 DO 3950 p = 1, NR 01155 TEMP1 = DNRM2(NR-p+1,WORK(2*N+(p-1)*NR+p),1) 01156 CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1) 01157 3950 CONTINUE 01158 CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1, 01159 & WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR) 01160 CONDR1 = ONE / DSQRT(TEMP1) 01161 * .. here need a second oppinion on the condition number 01162 * .. then assume worst case scenario 01163 * R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) 01164 * more conservative <=> CONDR1 .LT. DSQRT(DBLE(N)) 01165 * 01166 COND_OK = DSQRT(DBLE(NR)) 01167 *[TP] COND_OK is a tuning parameter. 01168 01169 IF ( CONDR1 .LT. COND_OK ) THEN 01170 * .. the second QRF without pivoting. Note: in an optimized 01171 * implementation, this QRF should be implemented as the QRF 01172 * of a lower triangular matrix. 01173 * R1^t = Q2 * R2 01174 CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), 01175 & LWORK-2*N, IERR ) 01176 * 01177 IF ( L2PERT ) THEN 01178 XSC = DSQRT(SMALL)/EPSLN 01179 DO 3959 p = 2, NR 01180 DO 3958 q = 1, p - 1 01181 TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) 01182 IF ( DABS(V(q,p)) .LE. TEMP1 ) 01183 & V(q,p) = DSIGN( TEMP1, V(q,p) ) 01184 3958 CONTINUE 01185 3959 CONTINUE 01186 END IF 01187 * 01188 IF ( NR .NE. N ) 01189 * .. save ... 01190 & CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N ) 01191 * 01192 * .. this transposed copy should be better than naive 01193 DO 1969 p = 1, NR - 1 01194 CALL DCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 ) 01195 1969 CONTINUE 01196 * 01197 CONDR2 = CONDR1 01198 * 01199 ELSE 01200 * 01201 * .. ill-conditioned case: second QRF with pivoting 01202 * Note that windowed pivoting would be equaly good 01203 * numerically, and more run-time efficient. So, in 01204 * an optimal implementation, the next call to DGEQP3 01205 * should be replaced with eg. CALL SGEQPX (ACM TOMS #782) 01206 * with properly (carefully) chosen parameters. 01207 * 01208 * R1^t * P2 = Q2 * R2 01209 DO 3003 p = 1, NR 01210 IWORK(N+p) = 0 01211 3003 CONTINUE 01212 CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1), 01213 & WORK(2*N+1), LWORK-2*N, IERR ) 01214 ** CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), 01215 ** & LWORK-2*N, IERR ) 01216 IF ( L2PERT ) THEN 01217 XSC = DSQRT(SMALL) 01218 DO 3969 p = 2, NR 01219 DO 3968 q = 1, p - 1 01220 TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) 01221 IF ( DABS(V(q,p)) .LE. TEMP1 ) 01222 & V(q,p) = DSIGN( TEMP1, V(q,p) ) 01223 3968 CONTINUE 01224 3969 CONTINUE 01225 END IF 01226 * 01227 CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N ) 01228 * 01229 IF ( L2PERT ) THEN 01230 XSC = DSQRT(SMALL) 01231 DO 8970 p = 2, NR 01232 DO 8971 q = 1, p - 1 01233 TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) 01234 V(p,q) = - DSIGN( TEMP1, V(q,p) ) 01235 8971 CONTINUE 01236 8970 CONTINUE 01237 ELSE 01238 CALL DLASET( 'L',NR-1,NR-1,ZERO,ZERO,V(2,1),LDV ) 01239 END IF 01240 * Now, compute R2 = L3 * Q3, the LQ factorization. 01241 CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1), 01242 & WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR ) 01243 * .. and estimate the condition number 01244 CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR ) 01245 DO 4950 p = 1, NR 01246 TEMP1 = DNRM2( p, WORK(2*N+N*NR+NR+p), NR ) 01247 CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR ) 01248 4950 CONTINUE 01249 CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1, 01250 & WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR ) 01251 CONDR2 = ONE / DSQRT(TEMP1) 01252 * 01253 IF ( CONDR2 .GE. COND_OK ) THEN 01254 * .. save the Householder vectors used for Q3 01255 * (this overwrittes the copy of R2, as it will not be 01256 * needed in this branch, but it does not overwritte the 01257 * Huseholder vectors of Q2.). 01258 CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N ) 01259 * .. and the rest of the information on Q3 is in 01260 * WORK(2*N+N*NR+1:2*N+N*NR+N) 01261 END IF 01262 * 01263 END IF 01264 * 01265 IF ( L2PERT ) THEN 01266 XSC = DSQRT(SMALL) 01267 DO 4968 q = 2, NR 01268 TEMP1 = XSC * V(q,q) 01269 DO 4969 p = 1, q - 1 01270 * V(p,q) = - DSIGN( TEMP1, V(q,p) ) 01271 V(p,q) = - DSIGN( TEMP1, V(p,q) ) 01272 4969 CONTINUE 01273 4968 CONTINUE 01274 ELSE 01275 CALL DLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV ) 01276 END IF 01277 * 01278 * Second preconditioning finished; continue with Jacobi SVD 01279 * The input matrix is lower trinagular. 01280 * 01281 * Recover the right singular vectors as solution of a well 01282 * conditioned triangular matrix equation. 01283 * 01284 IF ( CONDR1 .LT. COND_OK ) THEN 01285 * 01286 CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, 01287 & LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO ) 01288 SCALEM = WORK(2*N+N*NR+NR+1) 01289 NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) 01290 DO 3970 p = 1, NR 01291 CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 ) 01292 CALL DSCAL( NR, SVA(p), V(1,p), 1 ) 01293 3970 CONTINUE 01294 01295 * .. pick the right matrix equation and solve it 01296 * 01297 IF ( NR. EQ. N ) THEN 01298 * :)) .. best case, R1 is inverted. The solution of this matrix 01299 * equation is Q2*V2 = the product of the Jacobi rotations 01300 * used in DGESVJ, premultiplied with the orthogonal matrix 01301 * from the second QR factorization. 01302 CALL DTRSM( 'L','U','N','N', NR,NR,ONE, A,LDA, V,LDV ) 01303 ELSE 01304 * .. R1 is well conditioned, but non-square. Transpose(R2) 01305 * is inverted to get the product of the Jacobi rotations 01306 * used in DGESVJ. The Q-factor from the second QR 01307 * factorization is then built in explicitly. 01308 CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1), 01309 & N,V,LDV) 01310 IF ( NR .LT. N ) THEN 01311 CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV) 01312 CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV) 01313 CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV) 01314 END IF 01315 CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), 01316 & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) 01317 END IF 01318 * 01319 ELSE IF ( CONDR2 .LT. COND_OK ) THEN 01320 * 01321 * :) .. the input matrix A is very likely a relative of 01322 * the Kahan matrix :) 01323 * The matrix R2 is inverted. The solution of the matrix equation 01324 * is Q3^T*V3 = the product of the Jacobi rotations (appplied to 01325 * the lower triangular L3 from the LQ factorization of 01326 * R2=L3*Q3), pre-multiplied with the transposed Q3. 01327 CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U, 01328 & LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) 01329 SCALEM = WORK(2*N+N*NR+NR+1) 01330 NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) 01331 DO 3870 p = 1, NR 01332 CALL DCOPY( NR, V(1,p), 1, U(1,p), 1 ) 01333 CALL DSCAL( NR, SVA(p), U(1,p), 1 ) 01334 3870 CONTINUE 01335 CALL DTRSM('L','U','N','N',NR,NR,ONE,WORK(2*N+1),N,U,LDU) 01336 * .. apply the permutation from the second QR factorization 01337 DO 873 q = 1, NR 01338 DO 872 p = 1, NR 01339 WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) 01340 872 CONTINUE 01341 DO 874 p = 1, NR 01342 U(p,q) = WORK(2*N+N*NR+NR+p) 01343 874 CONTINUE 01344 873 CONTINUE 01345 IF ( NR .LT. N ) THEN 01346 CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) 01347 CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) 01348 CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) 01349 END IF 01350 CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), 01351 & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) 01352 ELSE 01353 * Last line of defense. 01354 * #:( This is a rather pathological case: no scaled condition 01355 * improvement after two pivoted QR factorizations. Other 01356 * possibility is that the rank revealing QR factorization 01357 * or the condition estimator has failed, or the COND_OK 01358 * is set very close to ONE (which is unnecessary). Normally, 01359 * this branch should never be executed, but in rare cases of 01360 * failure of the RRQR or condition estimator, the last line of 01361 * defense ensures that DGEJSV completes the task. 01362 * Compute the full SVD of L3 using DGESVJ with explicit 01363 * accumulation of Jacobi rotations. 01364 CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U, 01365 & LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) 01366 SCALEM = WORK(2*N+N*NR+NR+1) 01367 NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) 01368 IF ( NR .LT. N ) THEN 01369 CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) 01370 CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) 01371 CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) 01372 END IF 01373 CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), 01374 & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) 01375 * 01376 CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N, 01377 & WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1), 01378 & LWORK-2*N-N*NR-NR, IERR ) 01379 DO 773 q = 1, NR 01380 DO 772 p = 1, NR 01381 WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) 01382 772 CONTINUE 01383 DO 774 p = 1, NR 01384 U(p,q) = WORK(2*N+N*NR+NR+p) 01385 774 CONTINUE 01386 773 CONTINUE 01387 * 01388 END IF 01389 * 01390 * Permute the rows of V using the (column) permutation from the 01391 * first QRF. Also, scale the columns to make them unit in 01392 * Euclidean norm. This applies to all cases. 01393 * 01394 TEMP1 = DSQRT(DBLE(N)) * EPSLN 01395 DO 1972 q = 1, N 01396 DO 972 p = 1, N 01397 WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) 01398 972 CONTINUE 01399 DO 973 p = 1, N 01400 V(p,q) = WORK(2*N+N*NR+NR+p) 01401 973 CONTINUE 01402 XSC = ONE / DNRM2( N, V(1,q), 1 ) 01403 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) 01404 & CALL DSCAL( N, XSC, V(1,q), 1 ) 01405 1972 CONTINUE 01406 * At this moment, V contains the right singular vectors of A. 01407 * Next, assemble the left singular vector matrix U (M x N). 01408 IF ( NR .LT. M ) THEN 01409 CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU ) 01410 IF ( NR .LT. N1 ) THEN 01411 CALL DLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU) 01412 CALL DLASET('A',M-NR,N1-NR,ZERO,ONE,U(NR+1,NR+1),LDU) 01413 END IF 01414 END IF 01415 * 01416 * The Q matrix from the first QRF is built into the left singular 01417 * matrix U. This applies to all cases. 01418 * 01419 CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U, 01420 & LDU, WORK(N+1), LWORK-N, IERR ) 01421 01422 * The columns of U are normalized. The cost is O(M*N) flops. 01423 TEMP1 = DSQRT(DBLE(M)) * EPSLN 01424 DO 1973 p = 1, NR 01425 XSC = ONE / DNRM2( M, U(1,p), 1 ) 01426 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) 01427 & CALL DSCAL( M, XSC, U(1,p), 1 ) 01428 1973 CONTINUE 01429 * 01430 * If the initial QRF is computed with row pivoting, the left 01431 * singular vectors must be adjusted. 01432 * 01433 IF ( ROWPIV ) 01434 & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) 01435 * 01436 ELSE 01437 * 01438 * .. the initial matrix A has almost orthogonal columns and 01439 * the second QRF is not needed 01440 * 01441 CALL DLACPY( 'Upper', N, N, A, LDA, WORK(N+1), N ) 01442 IF ( L2PERT ) THEN 01443 XSC = DSQRT(SMALL) 01444 DO 5970 p = 2, N 01445 TEMP1 = XSC * WORK( N + (p-1)*N + p ) 01446 DO 5971 q = 1, p - 1 01447 WORK(N+(q-1)*N+p)=-DSIGN(TEMP1,WORK(N+(p-1)*N+q)) 01448 5971 CONTINUE 01449 5970 CONTINUE 01450 ELSE 01451 CALL DLASET( 'Lower',N-1,N-1,ZERO,ZERO,WORK(N+2),N ) 01452 END IF 01453 * 01454 CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA, 01455 & N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO ) 01456 * 01457 SCALEM = WORK(N+N*N+1) 01458 NUMRANK = IDNINT(WORK(N+N*N+2)) 01459 DO 6970 p = 1, N 01460 CALL DCOPY( N, WORK(N+(p-1)*N+1), 1, U(1,p), 1 ) 01461 CALL DSCAL( N, SVA(p), WORK(N+(p-1)*N+1), 1 ) 01462 6970 CONTINUE 01463 * 01464 CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N, 01465 & ONE, A, LDA, WORK(N+1), N ) 01466 DO 6972 p = 1, N 01467 CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV ) 01468 6972 CONTINUE 01469 TEMP1 = DSQRT(DBLE(N))*EPSLN 01470 DO 6971 p = 1, N 01471 XSC = ONE / DNRM2( N, V(1,p), 1 ) 01472 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) 01473 & CALL DSCAL( N, XSC, V(1,p), 1 ) 01474 6971 CONTINUE 01475 * 01476 * Assemble the left singular vector matrix U (M x N). 01477 * 01478 IF ( N .LT. M ) THEN 01479 CALL DLASET( 'A', M-N, N, ZERO, ZERO, U(N+1,1), LDU ) 01480 IF ( N .LT. N1 ) THEN 01481 CALL DLASET( 'A',N, N1-N, ZERO, ZERO, U(1,N+1),LDU ) 01482 CALL DLASET( 'A',M-N,N1-N, ZERO, ONE,U(N+1,N+1),LDU ) 01483 END IF 01484 END IF 01485 CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, 01486 & LDU, WORK(N+1), LWORK-N, IERR ) 01487 TEMP1 = DSQRT(DBLE(M))*EPSLN 01488 DO 6973 p = 1, N1 01489 XSC = ONE / DNRM2( M, U(1,p), 1 ) 01490 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) 01491 & CALL DSCAL( M, XSC, U(1,p), 1 ) 01492 6973 CONTINUE 01493 * 01494 IF ( ROWPIV ) 01495 & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) 01496 * 01497 END IF 01498 * 01499 * end of the >> almost orthogonal case << in the full SVD 01500 * 01501 ELSE 01502 * 01503 * This branch deploys a preconditioned Jacobi SVD with explicitly 01504 * accumulated rotations. It is included as optional, mainly for 01505 * experimental purposes. It does perfom well, and can also be used. 01506 * In this implementation, this branch will be automatically activated 01507 * if the condition number sigma_max(A) / sigma_min(A) is predicted 01508 * to be greater than the overflow threshold. This is because the 01509 * a posteriori computation of the singular vectors assumes robust 01510 * implementation of BLAS and some LAPACK procedures, capable of working 01511 * in presence of extreme values. Since that is not always the case, ... 01512 * 01513 DO 7968 p = 1, NR 01514 CALL DCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) 01515 7968 CONTINUE 01516 * 01517 IF ( L2PERT ) THEN 01518 XSC = DSQRT(SMALL/EPSLN) 01519 DO 5969 q = 1, NR 01520 TEMP1 = XSC*DABS( V(q,q) ) 01521 DO 5968 p = 1, N 01522 IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) 01523 & .OR. ( p .LT. q ) ) 01524 & V(p,q) = DSIGN( TEMP1, V(p,q) ) 01525 IF ( p. LT. q ) V(p,q) = - V(p,q) 01526 5968 CONTINUE 01527 5969 CONTINUE 01528 ELSE 01529 CALL DLASET( 'U', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) 01530 END IF 01531 01532 CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), 01533 & LWORK-2*N, IERR ) 01534 CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N ) 01535 * 01536 DO 7969 p = 1, NR 01537 CALL DCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 ) 01538 7969 CONTINUE 01539 01540 IF ( L2PERT ) THEN 01541 XSC = DSQRT(SMALL/EPSLN) 01542 DO 9970 q = 2, NR 01543 DO 9971 p = 1, q - 1 01544 TEMP1 = XSC * DMIN1(DABS(U(p,p)),DABS(U(q,q))) 01545 U(p,q) = - DSIGN( TEMP1, U(q,p) ) 01546 9971 CONTINUE 01547 9970 CONTINUE 01548 ELSE 01549 CALL DLASET('U', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) 01550 END IF 01551 01552 CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA, 01553 & N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO ) 01554 SCALEM = WORK(2*N+N*NR+1) 01555 NUMRANK = IDNINT(WORK(2*N+N*NR+2)) 01556 01557 IF ( NR .LT. N ) THEN 01558 CALL DLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV ) 01559 CALL DLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV ) 01560 CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) 01561 END IF 01562 01563 CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), 01564 & V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) 01565 * 01566 * Permute the rows of V using the (column) permutation from the 01567 * first QRF. Also, scale the columns to make them unit in 01568 * Euclidean norm. This applies to all cases. 01569 * 01570 TEMP1 = DSQRT(DBLE(N)) * EPSLN 01571 DO 7972 q = 1, N 01572 DO 8972 p = 1, N 01573 WORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) 01574 8972 CONTINUE 01575 DO 8973 p = 1, N 01576 V(p,q) = WORK(2*N+N*NR+NR+p) 01577 8973 CONTINUE 01578 XSC = ONE / DNRM2( N, V(1,q), 1 ) 01579 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) 01580 & CALL DSCAL( N, XSC, V(1,q), 1 ) 01581 7972 CONTINUE 01582 * 01583 * At this moment, V contains the right singular vectors of A. 01584 * Next, assemble the left singular vector matrix U (M x N). 01585 * 01586 IF ( NR .LT. M ) THEN 01587 CALL DLASET( 'A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU ) 01588 IF ( NR .LT. N1 ) THEN 01589 CALL DLASET( 'A',NR, N1-NR, ZERO, ZERO, U(1,NR+1),LDU ) 01590 CALL DLASET( 'A',M-NR,N1-NR, ZERO, ONE,U(NR+1,NR+1),LDU ) 01591 END IF 01592 END IF 01593 * 01594 CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, 01595 & LDU, WORK(N+1), LWORK-N, IERR ) 01596 * 01597 IF ( ROWPIV ) 01598 & CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) 01599 * 01600 * 01601 END IF 01602 IF ( TRANSP ) THEN 01603 * .. swap U and V because the procedure worked on A^t 01604 DO 6974 p = 1, N 01605 CALL DSWAP( N, U(1,p), 1, V(1,p), 1 ) 01606 6974 CONTINUE 01607 END IF 01608 * 01609 END IF 01610 * end of the full SVD 01611 * 01612 * Undo scaling, if necessary (and possible) 01613 * 01614 IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN 01615 CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR ) 01616 USCAL1 = ONE 01617 USCAL2 = ONE 01618 END IF 01619 * 01620 IF ( NR .LT. N ) THEN 01621 DO 3004 p = NR+1, N 01622 SVA(p) = ZERO 01623 3004 CONTINUE 01624 END IF 01625 * 01626 WORK(1) = USCAL2 * SCALEM 01627 WORK(2) = USCAL1 01628 IF ( ERREST ) WORK(3) = SCONDA 01629 IF ( LSVEC .AND. RSVEC ) THEN 01630 WORK(4) = CONDR1 01631 WORK(5) = CONDR2 01632 END IF 01633 IF ( L2TRAN ) THEN 01634 WORK(6) = ENTRA 01635 WORK(7) = ENTRAT 01636 END IF 01637 * 01638 IWORK(1) = NR 01639 IWORK(2) = NUMRANK 01640 IWORK(3) = WARNING 01641 * 01642 RETURN 01643 * .. 01644 * .. END OF DGEJSV 01645 * .. 01646 END 01647 *
1.7.3 
