![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAED3 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAED3 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed3.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed3.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed3.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, 00022 * CTOT, W, S, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, K, LDQ, N, N1 00026 * DOUBLE PRECISION RHO 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER CTOT( * ), INDX( * ) 00030 * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), 00031 * $ S( * ), W( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> DLAED3 finds the roots of the secular equation, as defined by the 00041 *> values in D, W, and RHO, between 1 and K. It makes the 00042 *> appropriate calls to DLAED4 and then updates the eigenvectors by 00043 *> multiplying the matrix of eigenvectors of the pair of eigensystems 00044 *> being combined by the matrix of eigenvectors of the K-by-K system 00045 *> which is solved here. 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 X-MP, Cray Y-MP, 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 *> \endverbatim 00054 * 00055 * Arguments: 00056 * ========== 00057 * 00058 *> \param[in] K 00059 *> \verbatim 00060 *> K is INTEGER 00061 *> The number of terms in the rational function to be solved by 00062 *> DLAED4. K >= 0. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] N 00066 *> \verbatim 00067 *> N is INTEGER 00068 *> The number of rows and columns in the Q matrix. 00069 *> N >= K (deflation may result in N>K). 00070 *> \endverbatim 00071 *> 00072 *> \param[in] N1 00073 *> \verbatim 00074 *> N1 is INTEGER 00075 *> The location of the last eigenvalue in the leading submatrix. 00076 *> min(1,N) <= N1 <= N/2. 00077 *> \endverbatim 00078 *> 00079 *> \param[out] D 00080 *> \verbatim 00081 *> D is DOUBLE PRECISION array, dimension (N) 00082 *> D(I) contains the updated eigenvalues for 00083 *> 1 <= I <= K. 00084 *> \endverbatim 00085 *> 00086 *> \param[out] Q 00087 *> \verbatim 00088 *> Q is DOUBLE PRECISION array, dimension (LDQ,N) 00089 *> Initially the first K columns are used as workspace. 00090 *> On output the columns 1 to K contain 00091 *> the updated eigenvectors. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] LDQ 00095 *> \verbatim 00096 *> LDQ is INTEGER 00097 *> The leading dimension of the array Q. LDQ >= max(1,N). 00098 *> \endverbatim 00099 *> 00100 *> \param[in] RHO 00101 *> \verbatim 00102 *> RHO is DOUBLE PRECISION 00103 *> The value of the parameter in the rank one update equation. 00104 *> RHO >= 0 required. 00105 *> \endverbatim 00106 *> 00107 *> \param[in,out] DLAMDA 00108 *> \verbatim 00109 *> DLAMDA 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. May be changed on output by 00113 *> having lowest order bit set to zero on Cray X-MP, Cray Y-MP, 00114 *> Cray-2, or Cray C-90, as described above. 00115 *> \endverbatim 00116 *> 00117 *> \param[in] Q2 00118 *> \verbatim 00119 *> Q2 is DOUBLE PRECISION array, dimension (LDQ2, N) 00120 *> The first K columns of this matrix contain the non-deflated 00121 *> eigenvectors for the split problem. 00122 *> \endverbatim 00123 *> 00124 *> \param[in] INDX 00125 *> \verbatim 00126 *> INDX is INTEGER array, dimension (N) 00127 *> The permutation used to arrange the columns of the deflated 00128 *> Q matrix into three groups (see DLAED2). 00129 *> The rows of the eigenvectors found by DLAED4 must be likewise 00130 *> permuted before the matrix multiply can take place. 00131 *> \endverbatim 00132 *> 00133 *> \param[in] CTOT 00134 *> \verbatim 00135 *> CTOT is INTEGER array, dimension (4) 00136 *> A count of the total number of the various types of columns 00137 *> in Q, as described in INDX. The fourth column type is any 00138 *> column which has been deflated. 00139 *> \endverbatim 00140 *> 00141 *> \param[in,out] W 00142 *> \verbatim 00143 *> W is DOUBLE PRECISION array, dimension (K) 00144 *> The first K elements of this array contain the components 00145 *> of the deflation-adjusted updating vector. Destroyed on 00146 *> output. 00147 *> \endverbatim 00148 *> 00149 *> \param[out] S 00150 *> \verbatim 00151 *> S is DOUBLE PRECISION array, dimension (N1 + 1)*K 00152 *> Will contain the eigenvectors of the repaired matrix which 00153 *> will be multiplied by the previously accumulated eigenvectors 00154 *> to update the system. 00155 *> \endverbatim 00156 *> 00157 *> \param[out] INFO 00158 *> \verbatim 00159 *> INFO is INTEGER 00160 *> = 0: successful exit. 00161 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00162 *> > 0: if INFO = 1, an eigenvalue did not converge 00163 *> \endverbatim 00164 * 00165 * Authors: 00166 * ======== 00167 * 00168 *> \author Univ. of Tennessee 00169 *> \author Univ. of California Berkeley 00170 *> \author Univ. of Colorado Denver 00171 *> \author NAG Ltd. 00172 * 00173 *> \date November 2011 00174 * 00175 *> \ingroup auxOTHERcomputational 00176 * 00177 *> \par Contributors: 00178 * ================== 00179 *> 00180 *> Jeff Rutter, Computer Science Division, University of California 00181 *> at Berkeley, USA \n 00182 *> Modified by Francoise Tisseur, University of Tennessee 00183 *> 00184 * ===================================================================== 00185 SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, 00186 $ CTOT, W, S, INFO ) 00187 * 00188 * -- LAPACK computational routine (version 3.4.0) -- 00189 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00191 * November 2011 00192 * 00193 * .. Scalar Arguments .. 00194 INTEGER INFO, K, LDQ, N, N1 00195 DOUBLE PRECISION RHO 00196 * .. 00197 * .. Array Arguments .. 00198 INTEGER CTOT( * ), INDX( * ) 00199 DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), 00200 $ S( * ), W( * ) 00201 * .. 00202 * 00203 * ===================================================================== 00204 * 00205 * .. Parameters .. 00206 DOUBLE PRECISION ONE, ZERO 00207 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) 00208 * .. 00209 * .. Local Scalars .. 00210 INTEGER I, II, IQ2, J, N12, N2, N23 00211 DOUBLE PRECISION TEMP 00212 * .. 00213 * .. External Functions .. 00214 DOUBLE PRECISION DLAMC3, DNRM2 00215 EXTERNAL DLAMC3, DNRM2 00216 * .. 00217 * .. External Subroutines .. 00218 EXTERNAL DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA 00219 * .. 00220 * .. Intrinsic Functions .. 00221 INTRINSIC MAX, SIGN, SQRT 00222 * .. 00223 * .. Executable Statements .. 00224 * 00225 * Test the input parameters. 00226 * 00227 INFO = 0 00228 * 00229 IF( K.LT.0 ) THEN 00230 INFO = -1 00231 ELSE IF( N.LT.K ) THEN 00232 INFO = -2 00233 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 00234 INFO = -6 00235 END IF 00236 IF( INFO.NE.0 ) THEN 00237 CALL XERBLA( 'DLAED3', -INFO ) 00238 RETURN 00239 END IF 00240 * 00241 * Quick return if possible 00242 * 00243 IF( K.EQ.0 ) 00244 $ RETURN 00245 * 00246 * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can 00247 * be computed with high relative accuracy (barring over/underflow). 00248 * This is a problem on machines without a guard digit in 00249 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). 00250 * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), 00251 * which on any of these machines zeros out the bottommost 00252 * bit of DLAMDA(I) if it is 1; this makes the subsequent 00253 * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation 00254 * occurs. On binary machines with a guard digit (almost all 00255 * machines) it does not change DLAMDA(I) at all. On hexadecimal 00256 * and decimal machines with a guard digit, it slightly 00257 * changes the bottommost bits of DLAMDA(I). It does not account 00258 * for hexadecimal or decimal machines without guard digits 00259 * (we know of none). We use a subroutine call to compute 00260 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating 00261 * this code. 00262 * 00263 DO 10 I = 1, K 00264 DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) 00265 10 CONTINUE 00266 * 00267 DO 20 J = 1, K 00268 CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) 00269 * 00270 * If the zero finder fails, the computation is terminated. 00271 * 00272 IF( INFO.NE.0 ) 00273 $ GO TO 120 00274 20 CONTINUE 00275 * 00276 IF( K.EQ.1 ) 00277 $ GO TO 110 00278 IF( K.EQ.2 ) THEN 00279 DO 30 J = 1, K 00280 W( 1 ) = Q( 1, J ) 00281 W( 2 ) = Q( 2, J ) 00282 II = INDX( 1 ) 00283 Q( 1, J ) = W( II ) 00284 II = INDX( 2 ) 00285 Q( 2, J ) = W( II ) 00286 30 CONTINUE 00287 GO TO 110 00288 END IF 00289 * 00290 * Compute updated W. 00291 * 00292 CALL DCOPY( K, W, 1, S, 1 ) 00293 * 00294 * Initialize W(I) = Q(I,I) 00295 * 00296 CALL DCOPY( K, Q, LDQ+1, W, 1 ) 00297 DO 60 J = 1, K 00298 DO 40 I = 1, J - 1 00299 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 00300 40 CONTINUE 00301 DO 50 I = J + 1, K 00302 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 00303 50 CONTINUE 00304 60 CONTINUE 00305 DO 70 I = 1, K 00306 W( I ) = SIGN( SQRT( -W( I ) ), S( I ) ) 00307 70 CONTINUE 00308 * 00309 * Compute eigenvectors of the modified rank-1 modification. 00310 * 00311 DO 100 J = 1, K 00312 DO 80 I = 1, K 00313 S( I ) = W( I ) / Q( I, J ) 00314 80 CONTINUE 00315 TEMP = DNRM2( K, S, 1 ) 00316 DO 90 I = 1, K 00317 II = INDX( I ) 00318 Q( I, J ) = S( II ) / TEMP 00319 90 CONTINUE 00320 100 CONTINUE 00321 * 00322 * Compute the updated eigenvectors. 00323 * 00324 110 CONTINUE 00325 * 00326 N2 = N - N1 00327 N12 = CTOT( 1 ) + CTOT( 2 ) 00328 N23 = CTOT( 2 ) + CTOT( 3 ) 00329 * 00330 CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 ) 00331 IQ2 = N1*N12 + 1 00332 IF( N23.NE.0 ) THEN 00333 CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23, 00334 $ ZERO, Q( N1+1, 1 ), LDQ ) 00335 ELSE 00336 CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ ) 00337 END IF 00338 * 00339 CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 ) 00340 IF( N12.NE.0 ) THEN 00341 CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q, 00342 $ LDQ ) 00343 ELSE 00344 CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ ) 00345 END IF 00346 * 00347 * 00348 120 CONTINUE 00349 RETURN 00350 * 00351 * End of DLAED3 00352 * 00353 END