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