LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlaed9.f
Go to the documentation of this file.
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlaed4(N, I, D, Z, DELTA, RHO, DLAM, INFO)
DLAED4 used by DSTEDC. Finds a single root of the secular equation.
Definition: dlaed4.f:145
subroutine dlaed9(K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, S, LDS, INFO)
DLAED9 used by DSTEDC. Finds the roots of the secular equation and updates the eigenvectors....
Definition: dlaed9.f:156
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82