Annotation of rpl/lapack/lapack/dlaed9.f, revision 1.13
1.11 bertrand 1: *> \brief \b DLAED9 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.
1.8 bertrand 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: *
1.11 bertrand 145: *> \date September 2012
1.8 bertrand 146: *
147: *> \ingroup auxOTHERcomputational
148: *
149: *> \par Contributors:
150: * ==================
151: *>
152: *> Jeff Rutter, Computer Science Division, University of California
153: *> at Berkeley, USA
154: *
155: * =====================================================================
1.1 bertrand 156: SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
157: $ S, LDS, INFO )
158: *
1.11 bertrand 159: * -- LAPACK computational routine (version 3.4.2) --
1.1 bertrand 160: * -- LAPACK is a software package provided by Univ. of Tennessee, --
161: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 162: * September 2012
1.1 bertrand 163: *
164: * .. Scalar Arguments ..
165: INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
166: DOUBLE PRECISION RHO
167: * ..
168: * .. Array Arguments ..
169: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
170: $ W( * )
171: * ..
172: *
173: * =====================================================================
174: *
175: * .. Local Scalars ..
176: INTEGER I, J
177: DOUBLE PRECISION TEMP
178: * ..
179: * .. External Functions ..
180: DOUBLE PRECISION DLAMC3, DNRM2
181: EXTERNAL DLAMC3, DNRM2
182: * ..
183: * .. External Subroutines ..
184: EXTERNAL DCOPY, DLAED4, XERBLA
185: * ..
186: * .. Intrinsic Functions ..
187: INTRINSIC MAX, SIGN, SQRT
188: * ..
189: * .. Executable Statements ..
190: *
191: * Test the input parameters.
192: *
193: INFO = 0
194: *
195: IF( K.LT.0 ) THEN
196: INFO = -1
197: ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
198: INFO = -2
199: ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
200: $ THEN
201: INFO = -3
202: ELSE IF( N.LT.K ) THEN
203: INFO = -4
204: ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
205: INFO = -7
206: ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
207: INFO = -12
208: END IF
209: IF( INFO.NE.0 ) THEN
210: CALL XERBLA( 'DLAED9', -INFO )
211: RETURN
212: END IF
213: *
214: * Quick return if possible
215: *
216: IF( K.EQ.0 )
217: $ RETURN
218: *
219: * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
220: * be computed with high relative accuracy (barring over/underflow).
221: * This is a problem on machines without a guard digit in
222: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
223: * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
224: * which on any of these machines zeros out the bottommost
225: * bit of DLAMDA(I) if it is 1; this makes the subsequent
226: * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
227: * occurs. On binary machines with a guard digit (almost all
228: * machines) it does not change DLAMDA(I) at all. On hexadecimal
229: * and decimal machines with a guard digit, it slightly
230: * changes the bottommost bits of DLAMDA(I). It does not account
231: * for hexadecimal or decimal machines without guard digits
232: * (we know of none). We use a subroutine call to compute
233: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
234: * this code.
235: *
236: DO 10 I = 1, N
237: DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
238: 10 CONTINUE
239: *
240: DO 20 J = KSTART, KSTOP
241: CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
242: *
243: * If the zero finder fails, the computation is terminated.
244: *
245: IF( INFO.NE.0 )
246: $ GO TO 120
247: 20 CONTINUE
248: *
249: IF( K.EQ.1 .OR. K.EQ.2 ) THEN
250: DO 40 I = 1, K
251: DO 30 J = 1, K
252: S( J, I ) = Q( J, I )
253: 30 CONTINUE
254: 40 CONTINUE
255: GO TO 120
256: END IF
257: *
258: * Compute updated W.
259: *
260: CALL DCOPY( K, W, 1, S, 1 )
261: *
262: * Initialize W(I) = Q(I,I)
263: *
264: CALL DCOPY( K, Q, LDQ+1, W, 1 )
265: DO 70 J = 1, K
266: DO 50 I = 1, J - 1
267: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
268: 50 CONTINUE
269: DO 60 I = J + 1, K
270: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
271: 60 CONTINUE
272: 70 CONTINUE
273: DO 80 I = 1, K
274: W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
275: 80 CONTINUE
276: *
277: * Compute eigenvectors of the modified rank-1 modification.
278: *
279: DO 110 J = 1, K
280: DO 90 I = 1, K
281: Q( I, J ) = W( I ) / Q( I, J )
282: 90 CONTINUE
283: TEMP = DNRM2( K, Q( 1, J ), 1 )
284: DO 100 I = 1, K
285: S( I, J ) = Q( I, J ) / TEMP
286: 100 CONTINUE
287: 110 CONTINUE
288: *
289: 120 CONTINUE
290: RETURN
291: *
292: * End of DLAED9
293: *
294: 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