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