![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLASD3 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLASD3 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd3.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd3.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd3.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, 00022 * LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, 00023 * INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR, 00027 * $ SQRE 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER CTOT( * ), IDXC( * ) 00031 * DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ), 00032 * $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ), 00033 * $ Z( * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> DLASD3 finds all the square roots of the roots of the secular 00043 *> equation, as defined by the values in D and Z. It makes the 00044 *> appropriate calls to DLASD4 and then updates the singular 00045 *> vectors by matrix multiplication. 00046 *> 00047 *> This code makes very mild assumptions about floating point 00048 *> arithmetic. It will work on machines with a guard digit in 00049 *> add/subtract, or on those binary machines without guard digits 00050 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. 00051 *> It could conceivably fail on hexadecimal or decimal machines 00052 *> without guard digits, but we know of none. 00053 *> 00054 *> DLASD3 is called from DLASD1. 00055 *> \endverbatim 00056 * 00057 * Arguments: 00058 * ========== 00059 * 00060 *> \param[in] NL 00061 *> \verbatim 00062 *> NL is INTEGER 00063 *> The row dimension of the upper block. NL >= 1. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] NR 00067 *> \verbatim 00068 *> NR is INTEGER 00069 *> The row dimension of the lower block. NR >= 1. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] SQRE 00073 *> \verbatim 00074 *> SQRE is INTEGER 00075 *> = 0: the lower block is an NR-by-NR square matrix. 00076 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00077 *> 00078 *> The bidiagonal matrix has N = NL + NR + 1 rows and 00079 *> M = N + SQRE >= N columns. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] K 00083 *> \verbatim 00084 *> K is INTEGER 00085 *> The size of the secular equation, 1 =< K = < N. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] D 00089 *> \verbatim 00090 *> D is DOUBLE PRECISION array, dimension(K) 00091 *> On exit the square roots of the roots of the secular equation, 00092 *> in ascending order. 00093 *> \endverbatim 00094 *> 00095 *> \param[out] Q 00096 *> \verbatim 00097 *> Q is DOUBLE PRECISION array, 00098 *> dimension at least (LDQ,K). 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDQ 00102 *> \verbatim 00103 *> LDQ is INTEGER 00104 *> The leading dimension of the array Q. LDQ >= K. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] DSIGMA 00108 *> \verbatim 00109 *> DSIGMA is DOUBLE PRECISION array, dimension(K) 00110 *> The first K elements of this array contain the old roots 00111 *> of the deflated updating problem. These are the poles 00112 *> of the secular equation. 00113 *> \endverbatim 00114 *> 00115 *> \param[out] U 00116 *> \verbatim 00117 *> U is DOUBLE PRECISION array, dimension (LDU, N) 00118 *> The last N - K columns of this matrix contain the deflated 00119 *> left singular vectors. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] LDU 00123 *> \verbatim 00124 *> LDU is INTEGER 00125 *> The leading dimension of the array U. LDU >= N. 00126 *> \endverbatim 00127 *> 00128 *> \param[in,out] U2 00129 *> \verbatim 00130 *> U2 is DOUBLE PRECISION array, dimension (LDU2, N) 00131 *> The first K columns of this matrix contain the non-deflated 00132 *> left singular vectors for the split problem. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] LDU2 00136 *> \verbatim 00137 *> LDU2 is INTEGER 00138 *> The leading dimension of the array U2. LDU2 >= N. 00139 *> \endverbatim 00140 *> 00141 *> \param[out] VT 00142 *> \verbatim 00143 *> VT is DOUBLE PRECISION array, dimension (LDVT, M) 00144 *> The last M - K columns of VT**T contain the deflated 00145 *> right singular vectors. 00146 *> \endverbatim 00147 *> 00148 *> \param[in] LDVT 00149 *> \verbatim 00150 *> LDVT is INTEGER 00151 *> The leading dimension of the array VT. LDVT >= N. 00152 *> \endverbatim 00153 *> 00154 *> \param[in,out] VT2 00155 *> \verbatim 00156 *> VT2 is DOUBLE PRECISION array, dimension (LDVT2, N) 00157 *> The first K columns of VT2**T contain the non-deflated 00158 *> right singular vectors for the split problem. 00159 *> \endverbatim 00160 *> 00161 *> \param[in] LDVT2 00162 *> \verbatim 00163 *> LDVT2 is INTEGER 00164 *> The leading dimension of the array VT2. LDVT2 >= N. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] IDXC 00168 *> \verbatim 00169 *> IDXC is INTEGER array, dimension ( N ) 00170 *> The permutation used to arrange the columns of U (and rows of 00171 *> VT) into three groups: the first group contains non-zero 00172 *> entries only at and above (or before) NL +1; the second 00173 *> contains non-zero entries only at and below (or after) NL+2; 00174 *> and the third is dense. The first column of U and the row of 00175 *> VT are treated separately, however. 00176 *> 00177 *> The rows of the singular vectors found by DLASD4 00178 *> must be likewise permuted before the matrix multiplies can 00179 *> take place. 00180 *> \endverbatim 00181 *> 00182 *> \param[in] CTOT 00183 *> \verbatim 00184 *> CTOT is INTEGER array, dimension ( 4 ) 00185 *> A count of the total number of the various types of columns 00186 *> in U (or rows in VT), as described in IDXC. The fourth column 00187 *> type is any column which has been deflated. 00188 *> \endverbatim 00189 *> 00190 *> \param[in] Z 00191 *> \verbatim 00192 *> Z is DOUBLE PRECISION array, dimension (K) 00193 *> The first K elements of this array contain the components 00194 *> of the deflation-adjusted updating row vector. 00195 *> \endverbatim 00196 *> 00197 *> \param[out] INFO 00198 *> \verbatim 00199 *> INFO is INTEGER 00200 *> = 0: successful exit. 00201 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00202 *> > 0: if INFO = 1, a singular value did not converge 00203 *> \endverbatim 00204 * 00205 * Authors: 00206 * ======== 00207 * 00208 *> \author Univ. of Tennessee 00209 *> \author Univ. of California Berkeley 00210 *> \author Univ. of Colorado Denver 00211 *> \author NAG Ltd. 00212 * 00213 *> \date November 2011 00214 * 00215 *> \ingroup auxOTHERauxiliary 00216 * 00217 *> \par Contributors: 00218 * ================== 00219 *> 00220 *> Ming Gu and Huan Ren, Computer Science Division, University of 00221 *> California at Berkeley, USA 00222 *> 00223 * ===================================================================== 00224 SUBROUTINE DLASD3( NL, NR, SQRE, K, D, Q, LDQ, DSIGMA, U, LDU, U2, 00225 $ LDU2, VT, LDVT, VT2, LDVT2, IDXC, CTOT, Z, 00226 $ INFO ) 00227 * 00228 * -- LAPACK auxiliary routine (version 3.4.0) -- 00229 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00230 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00231 * November 2011 00232 * 00233 * .. Scalar Arguments .. 00234 INTEGER INFO, K, LDQ, LDU, LDU2, LDVT, LDVT2, NL, NR, 00235 $ SQRE 00236 * .. 00237 * .. Array Arguments .. 00238 INTEGER CTOT( * ), IDXC( * ) 00239 DOUBLE PRECISION D( * ), DSIGMA( * ), Q( LDQ, * ), U( LDU, * ), 00240 $ U2( LDU2, * ), VT( LDVT, * ), VT2( LDVT2, * ), 00241 $ Z( * ) 00242 * .. 00243 * 00244 * ===================================================================== 00245 * 00246 * .. Parameters .. 00247 DOUBLE PRECISION ONE, ZERO, NEGONE 00248 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, 00249 $ NEGONE = -1.0D+0 ) 00250 * .. 00251 * .. Local Scalars .. 00252 INTEGER CTEMP, I, J, JC, KTEMP, M, N, NLP1, NLP2, NRP1 00253 DOUBLE PRECISION RHO, TEMP 00254 * .. 00255 * .. External Functions .. 00256 DOUBLE PRECISION DLAMC3, DNRM2 00257 EXTERNAL DLAMC3, DNRM2 00258 * .. 00259 * .. External Subroutines .. 00260 EXTERNAL DCOPY, DGEMM, DLACPY, DLASCL, DLASD4, XERBLA 00261 * .. 00262 * .. Intrinsic Functions .. 00263 INTRINSIC ABS, SIGN, SQRT 00264 * .. 00265 * .. Executable Statements .. 00266 * 00267 * Test the input parameters. 00268 * 00269 INFO = 0 00270 * 00271 IF( NL.LT.1 ) THEN 00272 INFO = -1 00273 ELSE IF( NR.LT.1 ) THEN 00274 INFO = -2 00275 ELSE IF( ( SQRE.NE.1 ) .AND. ( SQRE.NE.0 ) ) THEN 00276 INFO = -3 00277 END IF 00278 * 00279 N = NL + NR + 1 00280 M = N + SQRE 00281 NLP1 = NL + 1 00282 NLP2 = NL + 2 00283 * 00284 IF( ( K.LT.1 ) .OR. ( K.GT.N ) ) THEN 00285 INFO = -4 00286 ELSE IF( LDQ.LT.K ) THEN 00287 INFO = -7 00288 ELSE IF( LDU.LT.N ) THEN 00289 INFO = -10 00290 ELSE IF( LDU2.LT.N ) THEN 00291 INFO = -12 00292 ELSE IF( LDVT.LT.M ) THEN 00293 INFO = -14 00294 ELSE IF( LDVT2.LT.M ) THEN 00295 INFO = -16 00296 END IF 00297 IF( INFO.NE.0 ) THEN 00298 CALL XERBLA( 'DLASD3', -INFO ) 00299 RETURN 00300 END IF 00301 * 00302 * Quick return if possible 00303 * 00304 IF( K.EQ.1 ) THEN 00305 D( 1 ) = ABS( Z( 1 ) ) 00306 CALL DCOPY( M, VT2( 1, 1 ), LDVT2, VT( 1, 1 ), LDVT ) 00307 IF( Z( 1 ).GT.ZERO ) THEN 00308 CALL DCOPY( N, U2( 1, 1 ), 1, U( 1, 1 ), 1 ) 00309 ELSE 00310 DO 10 I = 1, N 00311 U( I, 1 ) = -U2( I, 1 ) 00312 10 CONTINUE 00313 END IF 00314 RETURN 00315 END IF 00316 * 00317 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can 00318 * be computed with high relative accuracy (barring over/underflow). 00319 * This is a problem on machines without a guard digit in 00320 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). 00321 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), 00322 * which on any of these machines zeros out the bottommost 00323 * bit of DSIGMA(I) if it is 1; this makes the subsequent 00324 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation 00325 * occurs. On binary machines with a guard digit (almost all 00326 * machines) it does not change DSIGMA(I) at all. On hexadecimal 00327 * and decimal machines with a guard digit, it slightly 00328 * changes the bottommost bits of DSIGMA(I). It does not account 00329 * for hexadecimal or decimal machines without guard digits 00330 * (we know of none). We use a subroutine call to compute 00331 * 2*DSIGMA(I) to prevent optimizing compilers from eliminating 00332 * this code. 00333 * 00334 DO 20 I = 1, K 00335 DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) 00336 20 CONTINUE 00337 * 00338 * Keep a copy of Z. 00339 * 00340 CALL DCOPY( K, Z, 1, Q, 1 ) 00341 * 00342 * Normalize Z. 00343 * 00344 RHO = DNRM2( K, Z, 1 ) 00345 CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO ) 00346 RHO = RHO*RHO 00347 * 00348 * Find the new singular values. 00349 * 00350 DO 30 J = 1, K 00351 CALL DLASD4( K, J, DSIGMA, Z, U( 1, J ), RHO, D( J ), 00352 $ VT( 1, J ), INFO ) 00353 * 00354 * If the zero finder fails, the computation is terminated. 00355 * 00356 IF( INFO.NE.0 ) THEN 00357 RETURN 00358 END IF 00359 30 CONTINUE 00360 * 00361 * Compute updated Z. 00362 * 00363 DO 60 I = 1, K 00364 Z( I ) = U( I, K )*VT( I, K ) 00365 DO 40 J = 1, I - 1 00366 Z( I ) = Z( I )*( U( I, J )*VT( I, J ) / 00367 $ ( DSIGMA( I )-DSIGMA( J ) ) / 00368 $ ( DSIGMA( I )+DSIGMA( J ) ) ) 00369 40 CONTINUE 00370 DO 50 J = I, K - 1 00371 Z( I ) = Z( I )*( U( I, J )*VT( I, J ) / 00372 $ ( DSIGMA( I )-DSIGMA( J+1 ) ) / 00373 $ ( DSIGMA( I )+DSIGMA( J+1 ) ) ) 00374 50 CONTINUE 00375 Z( I ) = SIGN( SQRT( ABS( Z( I ) ) ), Q( I, 1 ) ) 00376 60 CONTINUE 00377 * 00378 * Compute left singular vectors of the modified diagonal matrix, 00379 * and store related information for the right singular vectors. 00380 * 00381 DO 90 I = 1, K 00382 VT( 1, I ) = Z( 1 ) / U( 1, I ) / VT( 1, I ) 00383 U( 1, I ) = NEGONE 00384 DO 70 J = 2, K 00385 VT( J, I ) = Z( J ) / U( J, I ) / VT( J, I ) 00386 U( J, I ) = DSIGMA( J )*VT( J, I ) 00387 70 CONTINUE 00388 TEMP = DNRM2( K, U( 1, I ), 1 ) 00389 Q( 1, I ) = U( 1, I ) / TEMP 00390 DO 80 J = 2, K 00391 JC = IDXC( J ) 00392 Q( J, I ) = U( JC, I ) / TEMP 00393 80 CONTINUE 00394 90 CONTINUE 00395 * 00396 * Update the left singular vector matrix. 00397 * 00398 IF( K.EQ.2 ) THEN 00399 CALL DGEMM( 'N', 'N', N, K, K, ONE, U2, LDU2, Q, LDQ, ZERO, U, 00400 $ LDU ) 00401 GO TO 100 00402 END IF 00403 IF( CTOT( 1 ).GT.0 ) THEN 00404 CALL DGEMM( 'N', 'N', NL, K, CTOT( 1 ), ONE, U2( 1, 2 ), LDU2, 00405 $ Q( 2, 1 ), LDQ, ZERO, U( 1, 1 ), LDU ) 00406 IF( CTOT( 3 ).GT.0 ) THEN 00407 KTEMP = 2 + CTOT( 1 ) + CTOT( 2 ) 00408 CALL DGEMM( 'N', 'N', NL, K, CTOT( 3 ), ONE, U2( 1, KTEMP ), 00409 $ LDU2, Q( KTEMP, 1 ), LDQ, ONE, U( 1, 1 ), LDU ) 00410 END IF 00411 ELSE IF( CTOT( 3 ).GT.0 ) THEN 00412 KTEMP = 2 + CTOT( 1 ) + CTOT( 2 ) 00413 CALL DGEMM( 'N', 'N', NL, K, CTOT( 3 ), ONE, U2( 1, KTEMP ), 00414 $ LDU2, Q( KTEMP, 1 ), LDQ, ZERO, U( 1, 1 ), LDU ) 00415 ELSE 00416 CALL DLACPY( 'F', NL, K, U2, LDU2, U, LDU ) 00417 END IF 00418 CALL DCOPY( K, Q( 1, 1 ), LDQ, U( NLP1, 1 ), LDU ) 00419 KTEMP = 2 + CTOT( 1 ) 00420 CTEMP = CTOT( 2 ) + CTOT( 3 ) 00421 CALL DGEMM( 'N', 'N', NR, K, CTEMP, ONE, U2( NLP2, KTEMP ), LDU2, 00422 $ Q( KTEMP, 1 ), LDQ, ZERO, U( NLP2, 1 ), LDU ) 00423 * 00424 * Generate the right singular vectors. 00425 * 00426 100 CONTINUE 00427 DO 120 I = 1, K 00428 TEMP = DNRM2( K, VT( 1, I ), 1 ) 00429 Q( I, 1 ) = VT( 1, I ) / TEMP 00430 DO 110 J = 2, K 00431 JC = IDXC( J ) 00432 Q( I, J ) = VT( JC, I ) / TEMP 00433 110 CONTINUE 00434 120 CONTINUE 00435 * 00436 * Update the right singular vector matrix. 00437 * 00438 IF( K.EQ.2 ) THEN 00439 CALL DGEMM( 'N', 'N', K, M, K, ONE, Q, LDQ, VT2, LDVT2, ZERO, 00440 $ VT, LDVT ) 00441 RETURN 00442 END IF 00443 KTEMP = 1 + CTOT( 1 ) 00444 CALL DGEMM( 'N', 'N', K, NLP1, KTEMP, ONE, Q( 1, 1 ), LDQ, 00445 $ VT2( 1, 1 ), LDVT2, ZERO, VT( 1, 1 ), LDVT ) 00446 KTEMP = 2 + CTOT( 1 ) + CTOT( 2 ) 00447 IF( KTEMP.LE.LDVT2 ) 00448 $ CALL DGEMM( 'N', 'N', K, NLP1, CTOT( 3 ), ONE, Q( 1, KTEMP ), 00449 $ LDQ, VT2( KTEMP, 1 ), LDVT2, ONE, VT( 1, 1 ), 00450 $ LDVT ) 00451 * 00452 KTEMP = CTOT( 1 ) + 1 00453 NRP1 = NR + SQRE 00454 IF( KTEMP.GT.1 ) THEN 00455 DO 130 I = 1, K 00456 Q( I, KTEMP ) = Q( I, 1 ) 00457 130 CONTINUE 00458 DO 140 I = NLP2, M 00459 VT2( KTEMP, I ) = VT2( 1, I ) 00460 140 CONTINUE 00461 END IF 00462 CTEMP = 1 + CTOT( 2 ) + CTOT( 3 ) 00463 CALL DGEMM( 'N', 'N', K, NRP1, CTEMP, ONE, Q( 1, KTEMP ), LDQ, 00464 $ VT2( KTEMP, NLP2 ), LDVT2, ZERO, VT( 1, NLP2 ), LDVT ) 00465 * 00466 RETURN 00467 * 00468 * End of DLASD3 00469 * 00470 END