1: *> \brief \b DLAED9 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAED9 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
22: * S, LDS, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
26: * DOUBLE PRECISION RHO
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
30: * $ W( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLAED9 finds the roots of the secular equation, as defined by the
40: *> values in D, Z, and RHO, between KSTART and KSTOP. It makes the
41: *> appropriate calls to DLAED4 and then stores the new matrix of
42: *> eigenvectors for use in calculating the next level of Z vectors.
43: *> \endverbatim
44: *
45: * Arguments:
46: * ==========
47: *
48: *> \param[in] K
49: *> \verbatim
50: *> K is INTEGER
51: *> The number of terms in the rational function to be solved by
52: *> DLAED4. K >= 0.
53: *> \endverbatim
54: *>
55: *> \param[in] KSTART
56: *> \verbatim
57: *> KSTART is INTEGER
58: *> \endverbatim
59: *>
60: *> \param[in] KSTOP
61: *> \verbatim
62: *> KSTOP is INTEGER
63: *> The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
64: *> are to be computed. 1 <= KSTART <= KSTOP <= K.
65: *> \endverbatim
66: *>
67: *> \param[in] N
68: *> \verbatim
69: *> N is INTEGER
70: *> The number of rows and columns in the Q matrix.
71: *> N >= K (delation may result in N > K).
72: *> \endverbatim
73: *>
74: *> \param[out] D
75: *> \verbatim
76: *> D is DOUBLE PRECISION array, dimension (N)
77: *> D(I) contains the updated eigenvalues
78: *> for KSTART <= I <= KSTOP.
79: *> \endverbatim
80: *>
81: *> \param[out] Q
82: *> \verbatim
83: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
84: *> \endverbatim
85: *>
86: *> \param[in] LDQ
87: *> \verbatim
88: *> LDQ is INTEGER
89: *> The leading dimension of the array Q. LDQ >= max( 1, N ).
90: *> \endverbatim
91: *>
92: *> \param[in] RHO
93: *> \verbatim
94: *> RHO is DOUBLE PRECISION
95: *> The value of the parameter in the rank one update equation.
96: *> RHO >= 0 required.
97: *> \endverbatim
98: *>
99: *> \param[in] DLAMDA
100: *> \verbatim
101: *> DLAMDA is DOUBLE PRECISION array, dimension (K)
102: *> The first K elements of this array contain the old roots
103: *> of the deflated updating problem. These are the poles
104: *> of the secular equation.
105: *> \endverbatim
106: *>
107: *> \param[in] W
108: *> \verbatim
109: *> W is DOUBLE PRECISION array, dimension (K)
110: *> The first K elements of this array contain the components
111: *> of the deflation-adjusted updating vector.
112: *> \endverbatim
113: *>
114: *> \param[out] S
115: *> \verbatim
116: *> S is DOUBLE PRECISION array, dimension (LDS, K)
117: *> Will contain the eigenvectors of the repaired matrix which
118: *> will be stored for subsequent Z vector calculation and
119: *> multiplied by the previously accumulated eigenvectors
120: *> to update the system.
121: *> \endverbatim
122: *>
123: *> \param[in] LDS
124: *> \verbatim
125: *> LDS is INTEGER
126: *> The leading dimension of S. LDS >= max( 1, K ).
127: *> \endverbatim
128: *>
129: *> \param[out] INFO
130: *> \verbatim
131: *> INFO is INTEGER
132: *> = 0: successful exit.
133: *> < 0: if INFO = -i, the i-th argument had an illegal value.
134: *> > 0: if INFO = 1, an eigenvalue did not converge
135: *> \endverbatim
136: *
137: * Authors:
138: * ========
139: *
140: *> \author Univ. of Tennessee
141: *> \author Univ. of California Berkeley
142: *> \author Univ. of Colorado Denver
143: *> \author NAG Ltd.
144: *
145: *> \ingroup auxOTHERcomputational
146: *
147: *> \par Contributors:
148: * ==================
149: *>
150: *> Jeff Rutter, Computer Science Division, University of California
151: *> at Berkeley, USA
152: *
153: * =====================================================================
154: SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
155: $ S, LDS, INFO )
156: *
157: * -- LAPACK computational routine --
158: * -- LAPACK is a software package provided by Univ. of Tennessee, --
159: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160: *
161: * .. Scalar Arguments ..
162: INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
163: DOUBLE PRECISION RHO
164: * ..
165: * .. Array Arguments ..
166: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
167: $ W( * )
168: * ..
169: *
170: * =====================================================================
171: *
172: * .. Local Scalars ..
173: INTEGER I, J
174: DOUBLE PRECISION TEMP
175: * ..
176: * .. External Functions ..
177: DOUBLE PRECISION DLAMC3, DNRM2
178: EXTERNAL DLAMC3, DNRM2
179: * ..
180: * .. External Subroutines ..
181: EXTERNAL DCOPY, DLAED4, XERBLA
182: * ..
183: * .. Intrinsic Functions ..
184: INTRINSIC MAX, SIGN, SQRT
185: * ..
186: * .. Executable Statements ..
187: *
188: * Test the input parameters.
189: *
190: INFO = 0
191: *
192: IF( K.LT.0 ) THEN
193: INFO = -1
194: ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
195: INFO = -2
196: ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
197: $ THEN
198: INFO = -3
199: ELSE IF( N.LT.K ) THEN
200: INFO = -4
201: ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
202: INFO = -7
203: ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
204: INFO = -12
205: END IF
206: IF( INFO.NE.0 ) THEN
207: CALL XERBLA( 'DLAED9', -INFO )
208: RETURN
209: END IF
210: *
211: * Quick return if possible
212: *
213: IF( K.EQ.0 )
214: $ RETURN
215: *
216: * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
217: * be computed with high relative accuracy (barring over/underflow).
218: * This is a problem on machines without a guard digit in
219: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
220: * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
221: * which on any of these machines zeros out the bottommost
222: * bit of DLAMDA(I) if it is 1; this makes the subsequent
223: * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
224: * occurs. On binary machines with a guard digit (almost all
225: * machines) it does not change DLAMDA(I) at all. On hexadecimal
226: * and decimal machines with a guard digit, it slightly
227: * changes the bottommost bits of DLAMDA(I). It does not account
228: * for hexadecimal or decimal machines without guard digits
229: * (we know of none). We use a subroutine call to compute
230: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
231: * this code.
232: *
233: DO 10 I = 1, N
234: DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
235: 10 CONTINUE
236: *
237: DO 20 J = KSTART, KSTOP
238: CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
239: *
240: * If the zero finder fails, the computation is terminated.
241: *
242: IF( INFO.NE.0 )
243: $ GO TO 120
244: 20 CONTINUE
245: *
246: IF( K.EQ.1 .OR. K.EQ.2 ) THEN
247: DO 40 I = 1, K
248: DO 30 J = 1, K
249: S( J, I ) = Q( J, I )
250: 30 CONTINUE
251: 40 CONTINUE
252: GO TO 120
253: END IF
254: *
255: * Compute updated W.
256: *
257: CALL DCOPY( K, W, 1, S, 1 )
258: *
259: * Initialize W(I) = Q(I,I)
260: *
261: CALL DCOPY( K, Q, LDQ+1, W, 1 )
262: DO 70 J = 1, K
263: DO 50 I = 1, J - 1
264: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
265: 50 CONTINUE
266: DO 60 I = J + 1, K
267: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
268: 60 CONTINUE
269: 70 CONTINUE
270: DO 80 I = 1, K
271: W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
272: 80 CONTINUE
273: *
274: * Compute eigenvectors of the modified rank-1 modification.
275: *
276: DO 110 J = 1, K
277: DO 90 I = 1, K
278: Q( I, J ) = W( I ) / Q( I, J )
279: 90 CONTINUE
280: TEMP = DNRM2( K, Q( 1, J ), 1 )
281: DO 100 I = 1, K
282: S( I, J ) = Q( I, J ) / TEMP
283: 100 CONTINUE
284: 110 CONTINUE
285: *
286: 120 CONTINUE
287: RETURN
288: *
289: * End of DLAED9
290: *
291: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed9.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack