00001 SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, 00002 $ CTOT, W, S, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, K, LDQ, N, N1 00011 DOUBLE PRECISION RHO 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER CTOT( * ), INDX( * ) 00015 DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), 00016 $ S( * ), W( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DLAED3 finds the roots of the secular equation, as defined by the 00023 * values in D, W, and RHO, between 1 and K. It makes the 00024 * appropriate calls to DLAED4 and then updates the eigenvectors by 00025 * multiplying the matrix of eigenvectors of the pair of eigensystems 00026 * being combined by the matrix of eigenvectors of the K-by-K system 00027 * which is solved here. 00028 * 00029 * This code makes very mild assumptions about floating point 00030 * arithmetic. It will work on machines with a guard digit in 00031 * add/subtract, or on those binary machines without guard digits 00032 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 00033 * It could conceivably fail on hexadecimal or decimal machines 00034 * without guard digits, but we know of none. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * K (input) INTEGER 00040 * The number of terms in the rational function to be solved by 00041 * DLAED4. K >= 0. 00042 * 00043 * N (input) INTEGER 00044 * The number of rows and columns in the Q matrix. 00045 * N >= K (deflation may result in N>K). 00046 * 00047 * N1 (input) INTEGER 00048 * The location of the last eigenvalue in the leading submatrix. 00049 * min(1,N) <= N1 <= N/2. 00050 * 00051 * D (output) DOUBLE PRECISION array, dimension (N) 00052 * D(I) contains the updated eigenvalues for 00053 * 1 <= I <= K. 00054 * 00055 * Q (output) DOUBLE PRECISION array, dimension (LDQ,N) 00056 * Initially the first K columns are used as workspace. 00057 * On output the columns 1 to K contain 00058 * the updated eigenvectors. 00059 * 00060 * LDQ (input) INTEGER 00061 * The leading dimension of the array Q. LDQ >= max(1,N). 00062 * 00063 * RHO (input) DOUBLE PRECISION 00064 * The value of the parameter in the rank one update equation. 00065 * RHO >= 0 required. 00066 * 00067 * DLAMDA (input/output) DOUBLE PRECISION array, dimension (K) 00068 * The first K elements of this array contain the old roots 00069 * of the deflated updating problem. These are the poles 00070 * of the secular equation. May be changed on output by 00071 * having lowest order bit set to zero on Cray X-MP, Cray Y-MP, 00072 * Cray-2, or Cray C-90, as described above. 00073 * 00074 * Q2 (input) DOUBLE PRECISION array, dimension (LDQ2, N) 00075 * The first K columns of this matrix contain the non-deflated 00076 * eigenvectors for the split problem. 00077 * 00078 * INDX (input) INTEGER array, dimension (N) 00079 * The permutation used to arrange the columns of the deflated 00080 * Q matrix into three groups (see DLAED2). 00081 * The rows of the eigenvectors found by DLAED4 must be likewise 00082 * permuted before the matrix multiply can take place. 00083 * 00084 * CTOT (input) INTEGER array, dimension (4) 00085 * A count of the total number of the various types of columns 00086 * in Q, as described in INDX. The fourth column type is any 00087 * column which has been deflated. 00088 * 00089 * W (input/output) DOUBLE PRECISION array, dimension (K) 00090 * The first K elements of this array contain the components 00091 * of the deflation-adjusted updating vector. Destroyed on 00092 * output. 00093 * 00094 * S (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K 00095 * Will contain the eigenvectors of the repaired matrix which 00096 * will be multiplied by the previously accumulated eigenvectors 00097 * to update the system. 00098 * 00099 * LDS (input) INTEGER 00100 * The leading dimension of S. LDS >= max(1,K). 00101 * 00102 * INFO (output) INTEGER 00103 * = 0: successful exit. 00104 * < 0: if INFO = -i, the i-th argument had an illegal value. 00105 * > 0: if INFO = 1, an eigenvalue did not converge 00106 * 00107 * Further Details 00108 * =============== 00109 * 00110 * Based on contributions by 00111 * Jeff Rutter, Computer Science Division, University of California 00112 * at Berkeley, USA 00113 * Modified by Francoise Tisseur, University of Tennessee. 00114 * 00115 * ===================================================================== 00116 * 00117 * .. Parameters .. 00118 DOUBLE PRECISION ONE, ZERO 00119 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) 00120 * .. 00121 * .. Local Scalars .. 00122 INTEGER I, II, IQ2, J, N12, N2, N23 00123 DOUBLE PRECISION TEMP 00124 * .. 00125 * .. External Functions .. 00126 DOUBLE PRECISION DLAMC3, DNRM2 00127 EXTERNAL DLAMC3, DNRM2 00128 * .. 00129 * .. External Subroutines .. 00130 EXTERNAL DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA 00131 * .. 00132 * .. Intrinsic Functions .. 00133 INTRINSIC MAX, SIGN, SQRT 00134 * .. 00135 * .. Executable Statements .. 00136 * 00137 * Test the input parameters. 00138 * 00139 INFO = 0 00140 * 00141 IF( K.LT.0 ) THEN 00142 INFO = -1 00143 ELSE IF( N.LT.K ) THEN 00144 INFO = -2 00145 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 00146 INFO = -6 00147 END IF 00148 IF( INFO.NE.0 ) THEN 00149 CALL XERBLA( 'DLAED3', -INFO ) 00150 RETURN 00151 END IF 00152 * 00153 * Quick return if possible 00154 * 00155 IF( K.EQ.0 ) 00156 $ RETURN 00157 * 00158 * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can 00159 * be computed with high relative accuracy (barring over/underflow). 00160 * This is a problem on machines without a guard digit in 00161 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). 00162 * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), 00163 * which on any of these machines zeros out the bottommost 00164 * bit of DLAMDA(I) if it is 1; this makes the subsequent 00165 * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation 00166 * occurs. On binary machines with a guard digit (almost all 00167 * machines) it does not change DLAMDA(I) at all. On hexadecimal 00168 * and decimal machines with a guard digit, it slightly 00169 * changes the bottommost bits of DLAMDA(I). It does not account 00170 * for hexadecimal or decimal machines without guard digits 00171 * (we know of none). We use a subroutine call to compute 00172 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating 00173 * this code. 00174 * 00175 DO 10 I = 1, K 00176 DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) 00177 10 CONTINUE 00178 * 00179 DO 20 J = 1, K 00180 CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) 00181 * 00182 * If the zero finder fails, the computation is terminated. 00183 * 00184 IF( INFO.NE.0 ) 00185 $ GO TO 120 00186 20 CONTINUE 00187 * 00188 IF( K.EQ.1 ) 00189 $ GO TO 110 00190 IF( K.EQ.2 ) THEN 00191 DO 30 J = 1, K 00192 W( 1 ) = Q( 1, J ) 00193 W( 2 ) = Q( 2, J ) 00194 II = INDX( 1 ) 00195 Q( 1, J ) = W( II ) 00196 II = INDX( 2 ) 00197 Q( 2, J ) = W( II ) 00198 30 CONTINUE 00199 GO TO 110 00200 END IF 00201 * 00202 * Compute updated W. 00203 * 00204 CALL DCOPY( K, W, 1, S, 1 ) 00205 * 00206 * Initialize W(I) = Q(I,I) 00207 * 00208 CALL DCOPY( K, Q, LDQ+1, W, 1 ) 00209 DO 60 J = 1, K 00210 DO 40 I = 1, J - 1 00211 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 00212 40 CONTINUE 00213 DO 50 I = J + 1, K 00214 W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) 00215 50 CONTINUE 00216 60 CONTINUE 00217 DO 70 I = 1, K 00218 W( I ) = SIGN( SQRT( -W( I ) ), S( I ) ) 00219 70 CONTINUE 00220 * 00221 * Compute eigenvectors of the modified rank-1 modification. 00222 * 00223 DO 100 J = 1, K 00224 DO 80 I = 1, K 00225 S( I ) = W( I ) / Q( I, J ) 00226 80 CONTINUE 00227 TEMP = DNRM2( K, S, 1 ) 00228 DO 90 I = 1, K 00229 II = INDX( I ) 00230 Q( I, J ) = S( II ) / TEMP 00231 90 CONTINUE 00232 100 CONTINUE 00233 * 00234 * Compute the updated eigenvectors. 00235 * 00236 110 CONTINUE 00237 * 00238 N2 = N - N1 00239 N12 = CTOT( 1 ) + CTOT( 2 ) 00240 N23 = CTOT( 2 ) + CTOT( 3 ) 00241 * 00242 CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 ) 00243 IQ2 = N1*N12 + 1 00244 IF( N23.NE.0 ) THEN 00245 CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23, 00246 $ ZERO, Q( N1+1, 1 ), LDQ ) 00247 ELSE 00248 CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ ) 00249 END IF 00250 * 00251 CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 ) 00252 IF( N12.NE.0 ) THEN 00253 CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q, 00254 $ LDQ ) 00255 ELSE 00256 CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ ) 00257 END IF 00258 * 00259 * 00260 120 CONTINUE 00261 RETURN 00262 * 00263 * End of DLAED3 00264 * 00265 END
1.7.2