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
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