Annotation of rpl/lapack/lapack/dlaed9.f, revision 1.6
1.1 bertrand 1: SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
2: $ S, LDS, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
11: DOUBLE PRECISION RHO
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
15: $ W( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLAED9 finds the roots of the secular equation, as defined by the
22: * values in D, Z, and RHO, between KSTART and KSTOP. It makes the
23: * appropriate calls to DLAED4 and then stores the new matrix of
24: * eigenvectors for use in calculating the next level of Z vectors.
25: *
26: * Arguments
27: * =========
28: *
29: * K (input) INTEGER
30: * The number of terms in the rational function to be solved by
31: * DLAED4. K >= 0.
32: *
33: * KSTART (input) INTEGER
34: * KSTOP (input) INTEGER
35: * The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
36: * are to be computed. 1 <= KSTART <= KSTOP <= K.
37: *
38: * N (input) INTEGER
39: * The number of rows and columns in the Q matrix.
40: * N >= K (delation may result in N > K).
41: *
42: * D (output) DOUBLE PRECISION array, dimension (N)
43: * D(I) contains the updated eigenvalues
44: * for KSTART <= I <= KSTOP.
45: *
46: * Q (workspace) DOUBLE PRECISION array, dimension (LDQ,N)
47: *
48: * LDQ (input) INTEGER
49: * The leading dimension of the array Q. LDQ >= max( 1, N ).
50: *
51: * RHO (input) DOUBLE PRECISION
52: * The value of the parameter in the rank one update equation.
53: * RHO >= 0 required.
54: *
55: * DLAMDA (input) DOUBLE PRECISION array, dimension (K)
56: * The first K elements of this array contain the old roots
57: * of the deflated updating problem. These are the poles
58: * of the secular equation.
59: *
60: * W (input) DOUBLE PRECISION array, dimension (K)
61: * The first K elements of this array contain the components
62: * of the deflation-adjusted updating vector.
63: *
64: * S (output) DOUBLE PRECISION array, dimension (LDS, K)
65: * Will contain the eigenvectors of the repaired matrix which
66: * will be stored for subsequent Z vector calculation and
67: * multiplied by the previously accumulated eigenvectors
68: * to update the system.
69: *
70: * LDS (input) INTEGER
71: * The leading dimension of S. LDS >= max( 1, K ).
72: *
73: * INFO (output) INTEGER
74: * = 0: successful exit.
75: * < 0: if INFO = -i, the i-th argument had an illegal value.
76: * > 0: if INFO = 1, an eigenvalue did not converge
77: *
78: * Further Details
79: * ===============
80: *
81: * Based on contributions by
82: * Jeff Rutter, Computer Science Division, University of California
83: * at Berkeley, USA
84: *
85: * =====================================================================
86: *
87: * .. Local Scalars ..
88: INTEGER I, J
89: DOUBLE PRECISION TEMP
90: * ..
91: * .. External Functions ..
92: DOUBLE PRECISION DLAMC3, DNRM2
93: EXTERNAL DLAMC3, DNRM2
94: * ..
95: * .. External Subroutines ..
96: EXTERNAL DCOPY, DLAED4, XERBLA
97: * ..
98: * .. Intrinsic Functions ..
99: INTRINSIC MAX, SIGN, SQRT
100: * ..
101: * .. Executable Statements ..
102: *
103: * Test the input parameters.
104: *
105: INFO = 0
106: *
107: IF( K.LT.0 ) THEN
108: INFO = -1
109: ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
110: INFO = -2
111: ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
112: $ THEN
113: INFO = -3
114: ELSE IF( N.LT.K ) THEN
115: INFO = -4
116: ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
117: INFO = -7
118: ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
119: INFO = -12
120: END IF
121: IF( INFO.NE.0 ) THEN
122: CALL XERBLA( 'DLAED9', -INFO )
123: RETURN
124: END IF
125: *
126: * Quick return if possible
127: *
128: IF( K.EQ.0 )
129: $ RETURN
130: *
131: * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
132: * be computed with high relative accuracy (barring over/underflow).
133: * This is a problem on machines without a guard digit in
134: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
135: * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
136: * which on any of these machines zeros out the bottommost
137: * bit of DLAMDA(I) if it is 1; this makes the subsequent
138: * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
139: * occurs. On binary machines with a guard digit (almost all
140: * machines) it does not change DLAMDA(I) at all. On hexadecimal
141: * and decimal machines with a guard digit, it slightly
142: * changes the bottommost bits of DLAMDA(I). It does not account
143: * for hexadecimal or decimal machines without guard digits
144: * (we know of none). We use a subroutine call to compute
145: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
146: * this code.
147: *
148: DO 10 I = 1, N
149: DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
150: 10 CONTINUE
151: *
152: DO 20 J = KSTART, KSTOP
153: CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
154: *
155: * If the zero finder fails, the computation is terminated.
156: *
157: IF( INFO.NE.0 )
158: $ GO TO 120
159: 20 CONTINUE
160: *
161: IF( K.EQ.1 .OR. K.EQ.2 ) THEN
162: DO 40 I = 1, K
163: DO 30 J = 1, K
164: S( J, I ) = Q( J, I )
165: 30 CONTINUE
166: 40 CONTINUE
167: GO TO 120
168: END IF
169: *
170: * Compute updated W.
171: *
172: CALL DCOPY( K, W, 1, S, 1 )
173: *
174: * Initialize W(I) = Q(I,I)
175: *
176: CALL DCOPY( K, Q, LDQ+1, W, 1 )
177: DO 70 J = 1, K
178: DO 50 I = 1, J - 1
179: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
180: 50 CONTINUE
181: DO 60 I = J + 1, K
182: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
183: 60 CONTINUE
184: 70 CONTINUE
185: DO 80 I = 1, K
186: W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
187: 80 CONTINUE
188: *
189: * Compute eigenvectors of the modified rank-1 modification.
190: *
191: DO 110 J = 1, K
192: DO 90 I = 1, K
193: Q( I, J ) = W( I ) / Q( I, J )
194: 90 CONTINUE
195: TEMP = DNRM2( K, Q( 1, J ), 1 )
196: DO 100 I = 1, K
197: S( I, J ) = Q( I, J ) / TEMP
198: 100 CONTINUE
199: 110 CONTINUE
200: *
201: 120 CONTINUE
202: RETURN
203: *
204: * End of DLAED9
205: *
206: 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