LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlaed3.f
Go to the documentation of this file.
1 *> \brief \b DLAED3 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is tridiagonal.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLAED3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
22 * CTOT, W, S, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, K, LDQ, N, N1
26 * DOUBLE PRECISION RHO
27 * ..
28 * .. Array Arguments ..
29 * INTEGER CTOT( * ), INDX( * )
30 * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
31 * $ S( * ), W( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> DLAED3 finds the roots of the secular equation, as defined by the
41 *> values in D, W, and RHO, between 1 and K. It makes the
42 *> appropriate calls to DLAED4 and then updates the eigenvectors by
43 *> multiplying the matrix of eigenvectors of the pair of eigensystems
44 *> being combined by the matrix of eigenvectors of the K-by-K system
45 *> which is solved here.
46 *>
47 *> This code makes very mild assumptions about floating point
48 *> arithmetic. It will work on machines with a guard digit in
49 *> add/subtract, or on those binary machines without guard digits
50 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
51 *> It could conceivably fail on hexadecimal or decimal machines
52 *> without guard digits, but we know of none.
53 *> \endverbatim
54 *
55 * Arguments:
56 * ==========
57 *
58 *> \param[in] K
59 *> \verbatim
60 *> K is INTEGER
61 *> The number of terms in the rational function to be solved by
62 *> DLAED4. K >= 0.
63 *> \endverbatim
64 *>
65 *> \param[in] N
66 *> \verbatim
67 *> N is INTEGER
68 *> The number of rows and columns in the Q matrix.
69 *> N >= K (deflation may result in N>K).
70 *> \endverbatim
71 *>
72 *> \param[in] N1
73 *> \verbatim
74 *> N1 is INTEGER
75 *> The location of the last eigenvalue in the leading submatrix.
76 *> min(1,N) <= N1 <= N/2.
77 *> \endverbatim
78 *>
79 *> \param[out] D
80 *> \verbatim
81 *> D is DOUBLE PRECISION array, dimension (N)
82 *> D(I) contains the updated eigenvalues for
83 *> 1 <= I <= K.
84 *> \endverbatim
85 *>
86 *> \param[out] Q
87 *> \verbatim
88 *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
89 *> Initially the first K columns are used as workspace.
90 *> On output the columns 1 to K contain
91 *> the updated eigenvectors.
92 *> \endverbatim
93 *>
94 *> \param[in] LDQ
95 *> \verbatim
96 *> LDQ is INTEGER
97 *> The leading dimension of the array Q. LDQ >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[in] RHO
101 *> \verbatim
102 *> RHO is DOUBLE PRECISION
103 *> The value of the parameter in the rank one update equation.
104 *> RHO >= 0 required.
105 *> \endverbatim
106 *>
107 *> \param[in,out] DLAMDA
108 *> \verbatim
109 *> DLAMDA is DOUBLE PRECISION array, dimension (K)
110 *> The first K elements of this array contain the old roots
111 *> of the deflated updating problem. These are the poles
112 *> of the secular equation. May be changed on output by
113 *> having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
114 *> Cray-2, or Cray C-90, as described above.
115 *> \endverbatim
116 *>
117 *> \param[in] Q2
118 *> \verbatim
119 *> Q2 is DOUBLE PRECISION array, dimension (LDQ2, N)
120 *> The first K columns of this matrix contain the non-deflated
121 *> eigenvectors for the split problem.
122 *> \endverbatim
123 *>
124 *> \param[in] INDX
125 *> \verbatim
126 *> INDX is INTEGER array, dimension (N)
127 *> The permutation used to arrange the columns of the deflated
128 *> Q matrix into three groups (see DLAED2).
129 *> The rows of the eigenvectors found by DLAED4 must be likewise
130 *> permuted before the matrix multiply can take place.
131 *> \endverbatim
132 *>
133 *> \param[in] CTOT
134 *> \verbatim
135 *> CTOT is INTEGER array, dimension (4)
136 *> A count of the total number of the various types of columns
137 *> in Q, as described in INDX. The fourth column type is any
138 *> column which has been deflated.
139 *> \endverbatim
140 *>
141 *> \param[in,out] W
142 *> \verbatim
143 *> W is DOUBLE PRECISION array, dimension (K)
144 *> The first K elements of this array contain the components
145 *> of the deflation-adjusted updating vector. Destroyed on
146 *> output.
147 *> \endverbatim
148 *>
149 *> \param[out] S
150 *> \verbatim
151 *> S is DOUBLE PRECISION array, dimension (N1 + 1)*K
152 *> Will contain the eigenvectors of the repaired matrix which
153 *> will be multiplied by the previously accumulated eigenvectors
154 *> to update the system.
155 *> \endverbatim
156 *>
157 *> \param[out] INFO
158 *> \verbatim
159 *> INFO is INTEGER
160 *> = 0: successful exit.
161 *> < 0: if INFO = -i, the i-th argument had an illegal value.
162 *> > 0: if INFO = 1, an eigenvalue did not converge
163 *> \endverbatim
164 *
165 * Authors:
166 * ========
167 *
168 *> \author Univ. of Tennessee
169 *> \author Univ. of California Berkeley
170 *> \author Univ. of Colorado Denver
171 *> \author NAG Ltd.
172 *
173 *> \date September 2012
174 *
175 *> \ingroup auxOTHERcomputational
176 *
177 *> \par Contributors:
178 * ==================
179 *>
180 *> Jeff Rutter, Computer Science Division, University of California
181 *> at Berkeley, USA \n
182 *> Modified by Francoise Tisseur, University of Tennessee
183 *>
184 * =====================================================================
185  SUBROUTINE dlaed3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
186  $ ctot, w, s, info )
187 *
188 * -- LAPACK computational routine (version 3.4.2) --
189 * -- LAPACK is a software package provided by Univ. of Tennessee, --
190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191 * September 2012
192 *
193 * .. Scalar Arguments ..
194  INTEGER info, k, ldq, n, n1
195  DOUBLE PRECISION rho
196 * ..
197 * .. Array Arguments ..
198  INTEGER ctot( * ), indx( * )
199  DOUBLE PRECISION d( * ), dlamda( * ), q( ldq, * ), q2( * ),
200  $ s( * ), w( * )
201 * ..
202 *
203 * =====================================================================
204 *
205 * .. Parameters ..
206  DOUBLE PRECISION one, zero
207  parameter( one = 1.0d0, zero = 0.0d0 )
208 * ..
209 * .. Local Scalars ..
210  INTEGER i, ii, iq2, j, n12, n2, n23
211  DOUBLE PRECISION temp
212 * ..
213 * .. External Functions ..
214  DOUBLE PRECISION dlamc3, dnrm2
215  EXTERNAL dlamc3, dnrm2
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla
219 * ..
220 * .. Intrinsic Functions ..
221  INTRINSIC max, sign, sqrt
222 * ..
223 * .. Executable Statements ..
224 *
225 * Test the input parameters.
226 *
227  info = 0
228 *
229  IF( k.LT.0 ) THEN
230  info = -1
231  ELSE IF( n.LT.k ) THEN
232  info = -2
233  ELSE IF( ldq.LT.max( 1, n ) ) THEN
234  info = -6
235  END IF
236  IF( info.NE.0 ) THEN
237  CALL xerbla( 'DLAED3', -info )
238  return
239  END IF
240 *
241 * Quick return if possible
242 *
243  IF( k.EQ.0 )
244  $ return
245 *
246 * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
247 * be computed with high relative accuracy (barring over/underflow).
248 * This is a problem on machines without a guard digit in
249 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
250 * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
251 * which on any of these machines zeros out the bottommost
252 * bit of DLAMDA(I) if it is 1; this makes the subsequent
253 * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
254 * occurs. On binary machines with a guard digit (almost all
255 * machines) it does not change DLAMDA(I) at all. On hexadecimal
256 * and decimal machines with a guard digit, it slightly
257 * changes the bottommost bits of DLAMDA(I). It does not account
258 * for hexadecimal or decimal machines without guard digits
259 * (we know of none). We use a subroutine call to compute
260 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
261 * this code.
262 *
263  DO 10 i = 1, k
264  dlamda( i ) = dlamc3( dlamda( i ), dlamda( i ) ) - dlamda( i )
265  10 continue
266 *
267  DO 20 j = 1, k
268  CALL dlaed4( k, j, dlamda, w, q( 1, j ), rho, d( j ), info )
269 *
270 * If the zero finder fails, the computation is terminated.
271 *
272  IF( info.NE.0 )
273  $ go to 120
274  20 continue
275 *
276  IF( k.EQ.1 )
277  $ go to 110
278  IF( k.EQ.2 ) THEN
279  DO 30 j = 1, k
280  w( 1 ) = q( 1, j )
281  w( 2 ) = q( 2, j )
282  ii = indx( 1 )
283  q( 1, j ) = w( ii )
284  ii = indx( 2 )
285  q( 2, j ) = w( ii )
286  30 continue
287  go to 110
288  END IF
289 *
290 * Compute updated W.
291 *
292  CALL dcopy( k, w, 1, s, 1 )
293 *
294 * Initialize W(I) = Q(I,I)
295 *
296  CALL dcopy( k, q, ldq+1, w, 1 )
297  DO 60 j = 1, k
298  DO 40 i = 1, j - 1
299  w( i ) = w( i )*( q( i, j ) / ( dlamda( i )-dlamda( j ) ) )
300  40 continue
301  DO 50 i = j + 1, k
302  w( i ) = w( i )*( q( i, j ) / ( dlamda( i )-dlamda( j ) ) )
303  50 continue
304  60 continue
305  DO 70 i = 1, k
306  w( i ) = sign( sqrt( -w( i ) ), s( i ) )
307  70 continue
308 *
309 * Compute eigenvectors of the modified rank-1 modification.
310 *
311  DO 100 j = 1, k
312  DO 80 i = 1, k
313  s( i ) = w( i ) / q( i, j )
314  80 continue
315  temp = dnrm2( k, s, 1 )
316  DO 90 i = 1, k
317  ii = indx( i )
318  q( i, j ) = s( ii ) / temp
319  90 continue
320  100 continue
321 *
322 * Compute the updated eigenvectors.
323 *
324  110 continue
325 *
326  n2 = n - n1
327  n12 = ctot( 1 ) + ctot( 2 )
328  n23 = ctot( 2 ) + ctot( 3 )
329 *
330  CALL dlacpy( 'A', n23, k, q( ctot( 1 )+1, 1 ), ldq, s, n23 )
331  iq2 = n1*n12 + 1
332  IF( n23.NE.0 ) THEN
333  CALL dgemm( 'N', 'N', n2, k, n23, one, q2( iq2 ), n2, s, n23,
334  $ zero, q( n1+1, 1 ), ldq )
335  ELSE
336  CALL dlaset( 'A', n2, k, zero, zero, q( n1+1, 1 ), ldq )
337  END IF
338 *
339  CALL dlacpy( 'A', n12, k, q, ldq, s, n12 )
340  IF( n12.NE.0 ) THEN
341  CALL dgemm( 'N', 'N', n1, k, n12, one, q2, n1, s, n12, zero, q,
342  $ ldq )
343  ELSE
344  CALL dlaset( 'A', n1, k, zero, zero, q( 1, 1 ), ldq )
345  END IF
346 *
347 *
348  120 continue
349  return
350 *
351 * End of DLAED3
352 *
353  END