1: *> \brief \b DLAED3 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAED3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
22: * CTOT, W, S, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INFO, K, LDQ, N, N1
26: * DOUBLE PRECISION RHO
27: * ..
28: * .. Array Arguments ..
29: * INTEGER CTOT( * ), INDX( * )
30: * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
31: * $ S( * ), W( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DLAED3 finds the roots of the secular equation, as defined by the
41: *> values in D, W, and RHO, between 1 and K. It makes the
42: *> appropriate calls to DLAED4 and then updates the eigenvectors by
43: *> multiplying the matrix of eigenvectors of the pair of eigensystems
44: *> being combined by the matrix of eigenvectors of the K-by-K system
45: *> which is solved here.
46: *>
47: *> This code makes very mild assumptions about floating point
48: *> arithmetic. It will work on machines with a guard digit in
49: *> add/subtract, or on those binary machines without guard digits
50: *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
51: *> It could conceivably fail on hexadecimal or decimal machines
52: *> without guard digits, but we know of none.
53: *> \endverbatim
54: *
55: * Arguments:
56: * ==========
57: *
58: *> \param[in] K
59: *> \verbatim
60: *> K is INTEGER
61: *> The number of terms in the rational function to be solved by
62: *> DLAED4. K >= 0.
63: *> \endverbatim
64: *>
65: *> \param[in] N
66: *> \verbatim
67: *> N is INTEGER
68: *> The number of rows and columns in the Q matrix.
69: *> N >= K (deflation may result in N>K).
70: *> \endverbatim
71: *>
72: *> \param[in] N1
73: *> \verbatim
74: *> N1 is INTEGER
75: *> The location of the last eigenvalue in the leading submatrix.
76: *> min(1,N) <= N1 <= N/2.
77: *> \endverbatim
78: *>
79: *> \param[out] D
80: *> \verbatim
81: *> D is DOUBLE PRECISION array, dimension (N)
82: *> D(I) contains the updated eigenvalues for
83: *> 1 <= I <= K.
84: *> \endverbatim
85: *>
86: *> \param[out] Q
87: *> \verbatim
88: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
89: *> Initially the first K columns are used as workspace.
90: *> On output the columns 1 to K contain
91: *> the updated eigenvectors.
92: *> \endverbatim
93: *>
94: *> \param[in] LDQ
95: *> \verbatim
96: *> LDQ is INTEGER
97: *> The leading dimension of the array Q. LDQ >= max(1,N).
98: *> \endverbatim
99: *>
100: *> \param[in] RHO
101: *> \verbatim
102: *> RHO is DOUBLE PRECISION
103: *> The value of the parameter in the rank one update equation.
104: *> RHO >= 0 required.
105: *> \endverbatim
106: *>
107: *> \param[in,out] DLAMDA
108: *> \verbatim
109: *> DLAMDA is DOUBLE PRECISION array, dimension (K)
110: *> The first K elements of this array contain the old roots
111: *> of the deflated updating problem. These are the poles
112: *> of the secular equation. May be changed on output by
113: *> having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
114: *> Cray-2, or Cray C-90, as described above.
115: *> \endverbatim
116: *>
117: *> \param[in] Q2
118: *> \verbatim
119: *> Q2 is DOUBLE PRECISION array, dimension (LDQ2*N)
120: *> The first K columns of this matrix contain the non-deflated
121: *> eigenvectors for the split problem.
122: *> \endverbatim
123: *>
124: *> \param[in] INDX
125: *> \verbatim
126: *> INDX is INTEGER array, dimension (N)
127: *> The permutation used to arrange the columns of the deflated
128: *> Q matrix into three groups (see DLAED2).
129: *> The rows of the eigenvectors found by DLAED4 must be likewise
130: *> permuted before the matrix multiply can take place.
131: *> \endverbatim
132: *>
133: *> \param[in] CTOT
134: *> \verbatim
135: *> CTOT is INTEGER array, dimension (4)
136: *> A count of the total number of the various types of columns
137: *> in Q, as described in INDX. The fourth column type is any
138: *> column which has been deflated.
139: *> \endverbatim
140: *>
141: *> \param[in,out] W
142: *> \verbatim
143: *> W is DOUBLE PRECISION array, dimension (K)
144: *> The first K elements of this array contain the components
145: *> of the deflation-adjusted updating vector. Destroyed on
146: *> output.
147: *> \endverbatim
148: *>
149: *> \param[out] S
150: *> \verbatim
151: *> S is DOUBLE PRECISION array, dimension (N1 + 1)*K
152: *> Will contain the eigenvectors of the repaired matrix which
153: *> will be multiplied by the previously accumulated eigenvectors
154: *> to update the system.
155: *> \endverbatim
156: *>
157: *> \param[out] INFO
158: *> \verbatim
159: *> INFO is INTEGER
160: *> = 0: successful exit.
161: *> < 0: if INFO = -i, the i-th argument had an illegal value.
162: *> > 0: if INFO = 1, an eigenvalue did not converge
163: *> \endverbatim
164: *
165: * Authors:
166: * ========
167: *
168: *> \author Univ. of Tennessee
169: *> \author Univ. of California Berkeley
170: *> \author Univ. of Colorado Denver
171: *> \author NAG Ltd.
172: *
173: *> \ingroup auxOTHERcomputational
174: *
175: *> \par Contributors:
176: * ==================
177: *>
178: *> Jeff Rutter, Computer Science Division, University of California
179: *> at Berkeley, USA \n
180: *> Modified by Francoise Tisseur, University of Tennessee
181: *>
182: * =====================================================================
183: SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
184: $ CTOT, W, S, INFO )
185: *
186: * -- LAPACK computational routine --
187: * -- LAPACK is a software package provided by Univ. of Tennessee, --
188: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189: *
190: * .. Scalar Arguments ..
191: INTEGER INFO, K, LDQ, N, N1
192: DOUBLE PRECISION RHO
193: * ..
194: * .. Array Arguments ..
195: INTEGER CTOT( * ), INDX( * )
196: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
197: $ S( * ), W( * )
198: * ..
199: *
200: * =====================================================================
201: *
202: * .. Parameters ..
203: DOUBLE PRECISION ONE, ZERO
204: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
205: * ..
206: * .. Local Scalars ..
207: INTEGER I, II, IQ2, J, N12, N2, N23
208: DOUBLE PRECISION TEMP
209: * ..
210: * .. External Functions ..
211: DOUBLE PRECISION DLAMC3, DNRM2
212: EXTERNAL DLAMC3, DNRM2
213: * ..
214: * .. External Subroutines ..
215: EXTERNAL DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA
216: * ..
217: * .. Intrinsic Functions ..
218: INTRINSIC MAX, SIGN, SQRT
219: * ..
220: * .. Executable Statements ..
221: *
222: * Test the input parameters.
223: *
224: INFO = 0
225: *
226: IF( K.LT.0 ) THEN
227: INFO = -1
228: ELSE IF( N.LT.K ) THEN
229: INFO = -2
230: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
231: INFO = -6
232: END IF
233: IF( INFO.NE.0 ) THEN
234: CALL XERBLA( 'DLAED3', -INFO )
235: RETURN
236: END IF
237: *
238: * Quick return if possible
239: *
240: IF( K.EQ.0 )
241: $ RETURN
242: *
243: * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
244: * be computed with high relative accuracy (barring over/underflow).
245: * This is a problem on machines without a guard digit in
246: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
247: * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
248: * which on any of these machines zeros out the bottommost
249: * bit of DLAMDA(I) if it is 1; this makes the subsequent
250: * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
251: * occurs. On binary machines with a guard digit (almost all
252: * machines) it does not change DLAMDA(I) at all. On hexadecimal
253: * and decimal machines with a guard digit, it slightly
254: * changes the bottommost bits of DLAMDA(I). It does not account
255: * for hexadecimal or decimal machines without guard digits
256: * (we know of none). We use a subroutine call to compute
257: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
258: * this code.
259: *
260: DO 10 I = 1, K
261: DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
262: 10 CONTINUE
263: *
264: DO 20 J = 1, K
265: CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
266: *
267: * If the zero finder fails, the computation is terminated.
268: *
269: IF( INFO.NE.0 )
270: $ GO TO 120
271: 20 CONTINUE
272: *
273: IF( K.EQ.1 )
274: $ GO TO 110
275: IF( K.EQ.2 ) THEN
276: DO 30 J = 1, K
277: W( 1 ) = Q( 1, J )
278: W( 2 ) = Q( 2, J )
279: II = INDX( 1 )
280: Q( 1, J ) = W( II )
281: II = INDX( 2 )
282: Q( 2, J ) = W( II )
283: 30 CONTINUE
284: GO TO 110
285: END IF
286: *
287: * Compute updated W.
288: *
289: CALL DCOPY( K, W, 1, S, 1 )
290: *
291: * Initialize W(I) = Q(I,I)
292: *
293: CALL DCOPY( K, Q, LDQ+1, W, 1 )
294: DO 60 J = 1, K
295: DO 40 I = 1, J - 1
296: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
297: 40 CONTINUE
298: DO 50 I = J + 1, K
299: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
300: 50 CONTINUE
301: 60 CONTINUE
302: DO 70 I = 1, K
303: W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
304: 70 CONTINUE
305: *
306: * Compute eigenvectors of the modified rank-1 modification.
307: *
308: DO 100 J = 1, K
309: DO 80 I = 1, K
310: S( I ) = W( I ) / Q( I, J )
311: 80 CONTINUE
312: TEMP = DNRM2( K, S, 1 )
313: DO 90 I = 1, K
314: II = INDX( I )
315: Q( I, J ) = S( II ) / TEMP
316: 90 CONTINUE
317: 100 CONTINUE
318: *
319: * Compute the updated eigenvectors.
320: *
321: 110 CONTINUE
322: *
323: N2 = N - N1
324: N12 = CTOT( 1 ) + CTOT( 2 )
325: N23 = CTOT( 2 ) + CTOT( 3 )
326: *
327: CALL DLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
328: IQ2 = N1*N12 + 1
329: IF( N23.NE.0 ) THEN
330: CALL DGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
331: $ ZERO, Q( N1+1, 1 ), LDQ )
332: ELSE
333: CALL DLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
334: END IF
335: *
336: CALL DLACPY( 'A', N12, K, Q, LDQ, S, N12 )
337: IF( N12.NE.0 ) THEN
338: CALL DGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
339: $ LDQ )
340: ELSE
341: CALL DLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
342: END IF
343: *
344: *
345: 120 CONTINUE
346: RETURN
347: *
348: * End of DLAED3
349: *
350: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack