LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlaed2.f
Go to the documentation of this file.
1 *> \brief \b DLAED2 used by sstedc. Merges eigenvalues and deflates secular equation. 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 DLAED2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W,
22 * Q2, INDX, INDXC, INDXP, COLTYP, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER INFO, K, LDQ, N, N1
26 * DOUBLE PRECISION RHO
27 * ..
28 * .. Array Arguments ..
29 * INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ),
30 * $ INDXQ( * )
31 * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
32 * $ W( * ), Z( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DLAED2 merges the two sets of eigenvalues together into a single
42 *> sorted set. Then it tries to deflate the size of the problem.
43 *> There are two ways in which deflation can occur: when two or more
44 *> eigenvalues are close together or if there is a tiny entry in the
45 *> Z vector. For each such occurrence the order of the related secular
46 *> equation problem is reduced by one.
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \param[out] K
53 *> \verbatim
54 *> K is INTEGER
55 *> The number of non-deflated eigenvalues, and the order of the
56 *> related secular equation. 0 <= K <=N.
57 *> \endverbatim
58 *>
59 *> \param[in] N
60 *> \verbatim
61 *> N is INTEGER
62 *> The dimension of the symmetric tridiagonal matrix. N >= 0.
63 *> \endverbatim
64 *>
65 *> \param[in] N1
66 *> \verbatim
67 *> N1 is INTEGER
68 *> The location of the last eigenvalue in the leading sub-matrix.
69 *> min(1,N) <= N1 <= N/2.
70 *> \endverbatim
71 *>
72 *> \param[in,out] D
73 *> \verbatim
74 *> D is DOUBLE PRECISION array, dimension (N)
75 *> On entry, D contains the eigenvalues of the two submatrices to
76 *> be combined.
77 *> On exit, D contains the trailing (N-K) updated eigenvalues
78 *> (those which were deflated) sorted into increasing order.
79 *> \endverbatim
80 *>
81 *> \param[in,out] Q
82 *> \verbatim
83 *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
84 *> On entry, Q contains the eigenvectors of two submatrices in
85 *> the two square blocks with corners at (1,1), (N1,N1)
86 *> and (N1+1, N1+1), (N,N).
87 *> On exit, Q contains the trailing (N-K) updated eigenvectors
88 *> (those which were deflated) in its last N-K columns.
89 *> \endverbatim
90 *>
91 *> \param[in] LDQ
92 *> \verbatim
93 *> LDQ is INTEGER
94 *> The leading dimension of the array Q. LDQ >= max(1,N).
95 *> \endverbatim
96 *>
97 *> \param[in,out] INDXQ
98 *> \verbatim
99 *> INDXQ is INTEGER array, dimension (N)
100 *> The permutation which separately sorts the two sub-problems
101 *> in D into ascending order. Note that elements in the second
102 *> half of this permutation must first have N1 added to their
103 *> values. Destroyed on exit.
104 *> \endverbatim
105 *>
106 *> \param[in,out] RHO
107 *> \verbatim
108 *> RHO is DOUBLE PRECISION
109 *> On entry, the off-diagonal element associated with the rank-1
110 *> cut which originally split the two submatrices which are now
111 *> being recombined.
112 *> On exit, RHO has been modified to the value required by
113 *> DLAED3.
114 *> \endverbatim
115 *>
116 *> \param[in] Z
117 *> \verbatim
118 *> Z is DOUBLE PRECISION array, dimension (N)
119 *> On entry, Z contains the updating vector (the last
120 *> row of the first sub-eigenvector matrix and the first row of
121 *> the second sub-eigenvector matrix).
122 *> On exit, the contents of Z have been destroyed by the updating
123 *> process.
124 *> \endverbatim
125 *>
126 *> \param[out] DLAMDA
127 *> \verbatim
128 *> DLAMDA is DOUBLE PRECISION array, dimension (N)
129 *> A copy of the first K eigenvalues which will be used by
130 *> DLAED3 to form the secular equation.
131 *> \endverbatim
132 *>
133 *> \param[out] W
134 *> \verbatim
135 *> W is DOUBLE PRECISION array, dimension (N)
136 *> The first k values of the final deflation-altered z-vector
137 *> which will be passed to DLAED3.
138 *> \endverbatim
139 *>
140 *> \param[out] Q2
141 *> \verbatim
142 *> Q2 is DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2)
143 *> A copy of the first K eigenvectors which will be used by
144 *> DLAED3 in a matrix multiply (DGEMM) to solve for the new
145 *> eigenvectors.
146 *> \endverbatim
147 *>
148 *> \param[out] INDX
149 *> \verbatim
150 *> INDX is INTEGER array, dimension (N)
151 *> The permutation used to sort the contents of DLAMDA into
152 *> ascending order.
153 *> \endverbatim
154 *>
155 *> \param[out] INDXC
156 *> \verbatim
157 *> INDXC is INTEGER array, dimension (N)
158 *> The permutation used to arrange the columns of the deflated
159 *> Q matrix into three groups: the first group contains non-zero
160 *> elements only at and above N1, the second contains
161 *> non-zero elements only below N1, and the third is dense.
162 *> \endverbatim
163 *>
164 *> \param[out] INDXP
165 *> \verbatim
166 *> INDXP is INTEGER array, dimension (N)
167 *> The permutation used to place deflated values of D at the end
168 *> of the array. INDXP(1:K) points to the nondeflated D-values
169 *> and INDXP(K+1:N) points to the deflated eigenvalues.
170 *> \endverbatim
171 *>
172 *> \param[out] COLTYP
173 *> \verbatim
174 *> COLTYP is INTEGER array, dimension (N)
175 *> During execution, a label which will indicate which of the
176 *> following types a column in the Q2 matrix is:
177 *> 1 : non-zero in the upper half only;
178 *> 2 : dense;
179 *> 3 : non-zero in the lower half only;
180 *> 4 : deflated.
181 *> On exit, COLTYP(i) is the number of columns of type i,
182 *> for i=1 to 4 only.
183 *> \endverbatim
184 *>
185 *> \param[out] INFO
186 *> \verbatim
187 *> INFO is INTEGER
188 *> = 0: successful exit.
189 *> < 0: if INFO = -i, the i-th argument had an illegal value.
190 *> \endverbatim
191 *
192 * Authors:
193 * ========
194 *
195 *> \author Univ. of Tennessee
196 *> \author Univ. of California Berkeley
197 *> \author Univ. of Colorado Denver
198 *> \author NAG Ltd.
199 *
200 *> \date September 2012
201 *
202 *> \ingroup auxOTHERcomputational
203 *
204 *> \par Contributors:
205 * ==================
206 *>
207 *> Jeff Rutter, Computer Science Division, University of California
208 *> at Berkeley, USA \n
209 *> Modified by Francoise Tisseur, University of Tennessee
210 *>
211 * =====================================================================
212  SUBROUTINE dlaed2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W,
213  $ q2, indx, indxc, indxp, coltyp, info )
214 *
215 * -- LAPACK computational routine (version 3.4.2) --
216 * -- LAPACK is a software package provided by Univ. of Tennessee, --
217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218 * September 2012
219 *
220 * .. Scalar Arguments ..
221  INTEGER info, k, ldq, n, n1
222  DOUBLE PRECISION rho
223 * ..
224 * .. Array Arguments ..
225  INTEGER coltyp( * ), indx( * ), indxc( * ), indxp( * ),
226  $ indxq( * )
227  DOUBLE PRECISION d( * ), dlamda( * ), q( ldq, * ), q2( * ),
228  $ w( * ), z( * )
229 * ..
230 *
231 * =====================================================================
232 *
233 * .. Parameters ..
234  DOUBLE PRECISION mone, zero, one, two, eight
235  parameter( mone = -1.0d0, zero = 0.0d0, one = 1.0d0,
236  $ two = 2.0d0, eight = 8.0d0 )
237 * ..
238 * .. Local Arrays ..
239  INTEGER ctot( 4 ), psm( 4 )
240 * ..
241 * .. Local Scalars ..
242  INTEGER ct, i, imax, iq1, iq2, j, jmax, js, k2, n1p1,
243  $ n2, nj, pj
244  DOUBLE PRECISION c, eps, s, t, tau, tol
245 * ..
246 * .. External Functions ..
247  INTEGER idamax
248  DOUBLE PRECISION dlamch, dlapy2
249  EXTERNAL idamax, dlamch, dlapy2
250 * ..
251 * .. External Subroutines ..
252  EXTERNAL dcopy, dlacpy, dlamrg, drot, dscal, xerbla
253 * ..
254 * .. Intrinsic Functions ..
255  INTRINSIC abs, max, min, sqrt
256 * ..
257 * .. Executable Statements ..
258 *
259 * Test the input parameters.
260 *
261  info = 0
262 *
263  IF( n.LT.0 ) THEN
264  info = -2
265  ELSE IF( ldq.LT.max( 1, n ) ) THEN
266  info = -6
267  ELSE IF( min( 1, ( n / 2 ) ).GT.n1 .OR. ( n / 2 ).LT.n1 ) THEN
268  info = -3
269  END IF
270  IF( info.NE.0 ) THEN
271  CALL xerbla( 'DLAED2', -info )
272  return
273  END IF
274 *
275 * Quick return if possible
276 *
277  IF( n.EQ.0 )
278  $ return
279 *
280  n2 = n - n1
281  n1p1 = n1 + 1
282 *
283  IF( rho.LT.zero ) THEN
284  CALL dscal( n2, mone, z( n1p1 ), 1 )
285  END IF
286 *
287 * Normalize z so that norm(z) = 1. Since z is the concatenation of
288 * two normalized vectors, norm2(z) = sqrt(2).
289 *
290  t = one / sqrt( two )
291  CALL dscal( n, t, z, 1 )
292 *
293 * RHO = ABS( norm(z)**2 * RHO )
294 *
295  rho = abs( two*rho )
296 *
297 * Sort the eigenvalues into increasing order
298 *
299  DO 10 i = n1p1, n
300  indxq( i ) = indxq( i ) + n1
301  10 continue
302 *
303 * re-integrate the deflated parts from the last pass
304 *
305  DO 20 i = 1, n
306  dlamda( i ) = d( indxq( i ) )
307  20 continue
308  CALL dlamrg( n1, n2, dlamda, 1, 1, indxc )
309  DO 30 i = 1, n
310  indx( i ) = indxq( indxc( i ) )
311  30 continue
312 *
313 * Calculate the allowable deflation tolerance
314 *
315  imax = idamax( n, z, 1 )
316  jmax = idamax( n, d, 1 )
317  eps = dlamch( 'Epsilon' )
318  tol = eight*eps*max( abs( d( jmax ) ), abs( z( imax ) ) )
319 *
320 * If the rank-1 modifier is small enough, no more needs to be done
321 * except to reorganize Q so that its columns correspond with the
322 * elements in D.
323 *
324  IF( rho*abs( z( imax ) ).LE.tol ) THEN
325  k = 0
326  iq2 = 1
327  DO 40 j = 1, n
328  i = indx( j )
329  CALL dcopy( n, q( 1, i ), 1, q2( iq2 ), 1 )
330  dlamda( j ) = d( i )
331  iq2 = iq2 + n
332  40 continue
333  CALL dlacpy( 'A', n, n, q2, n, q, ldq )
334  CALL dcopy( n, dlamda, 1, d, 1 )
335  go to 190
336  END IF
337 *
338 * If there are multiple eigenvalues then the problem deflates. Here
339 * the number of equal eigenvalues are found. As each equal
340 * eigenvalue is found, an elementary reflector is computed to rotate
341 * the corresponding eigensubspace so that the corresponding
342 * components of Z are zero in this new basis.
343 *
344  DO 50 i = 1, n1
345  coltyp( i ) = 1
346  50 continue
347  DO 60 i = n1p1, n
348  coltyp( i ) = 3
349  60 continue
350 *
351 *
352  k = 0
353  k2 = n + 1
354  DO 70 j = 1, n
355  nj = indx( j )
356  IF( rho*abs( z( nj ) ).LE.tol ) THEN
357 *
358 * Deflate due to small z component.
359 *
360  k2 = k2 - 1
361  coltyp( nj ) = 4
362  indxp( k2 ) = nj
363  IF( j.EQ.n )
364  $ go to 100
365  ELSE
366  pj = nj
367  go to 80
368  END IF
369  70 continue
370  80 continue
371  j = j + 1
372  nj = indx( j )
373  IF( j.GT.n )
374  $ go to 100
375  IF( rho*abs( z( nj ) ).LE.tol ) THEN
376 *
377 * Deflate due to small z component.
378 *
379  k2 = k2 - 1
380  coltyp( nj ) = 4
381  indxp( k2 ) = nj
382  ELSE
383 *
384 * Check if eigenvalues are close enough to allow deflation.
385 *
386  s = z( pj )
387  c = z( nj )
388 *
389 * Find sqrt(a**2+b**2) without overflow or
390 * destructive underflow.
391 *
392  tau = dlapy2( c, s )
393  t = d( nj ) - d( pj )
394  c = c / tau
395  s = -s / tau
396  IF( abs( t*c*s ).LE.tol ) THEN
397 *
398 * Deflation is possible.
399 *
400  z( nj ) = tau
401  z( pj ) = zero
402  IF( coltyp( nj ).NE.coltyp( pj ) )
403  $ coltyp( nj ) = 2
404  coltyp( pj ) = 4
405  CALL drot( n, q( 1, pj ), 1, q( 1, nj ), 1, c, s )
406  t = d( pj )*c**2 + d( nj )*s**2
407  d( nj ) = d( pj )*s**2 + d( nj )*c**2
408  d( pj ) = t
409  k2 = k2 - 1
410  i = 1
411  90 continue
412  IF( k2+i.LE.n ) THEN
413  IF( d( pj ).LT.d( indxp( k2+i ) ) ) THEN
414  indxp( k2+i-1 ) = indxp( k2+i )
415  indxp( k2+i ) = pj
416  i = i + 1
417  go to 90
418  ELSE
419  indxp( k2+i-1 ) = pj
420  END IF
421  ELSE
422  indxp( k2+i-1 ) = pj
423  END IF
424  pj = nj
425  ELSE
426  k = k + 1
427  dlamda( k ) = d( pj )
428  w( k ) = z( pj )
429  indxp( k ) = pj
430  pj = nj
431  END IF
432  END IF
433  go to 80
434  100 continue
435 *
436 * Record the last eigenvalue.
437 *
438  k = k + 1
439  dlamda( k ) = d( pj )
440  w( k ) = z( pj )
441  indxp( k ) = pj
442 *
443 * Count up the total number of the various types of columns, then
444 * form a permutation which positions the four column types into
445 * four uniform groups (although one or more of these groups may be
446 * empty).
447 *
448  DO 110 j = 1, 4
449  ctot( j ) = 0
450  110 continue
451  DO 120 j = 1, n
452  ct = coltyp( j )
453  ctot( ct ) = ctot( ct ) + 1
454  120 continue
455 *
456 * PSM(*) = Position in SubMatrix (of types 1 through 4)
457 *
458  psm( 1 ) = 1
459  psm( 2 ) = 1 + ctot( 1 )
460  psm( 3 ) = psm( 2 ) + ctot( 2 )
461  psm( 4 ) = psm( 3 ) + ctot( 3 )
462  k = n - ctot( 4 )
463 *
464 * Fill out the INDXC array so that the permutation which it induces
465 * will place all type-1 columns first, all type-2 columns next,
466 * then all type-3's, and finally all type-4's.
467 *
468  DO 130 j = 1, n
469  js = indxp( j )
470  ct = coltyp( js )
471  indx( psm( ct ) ) = js
472  indxc( psm( ct ) ) = j
473  psm( ct ) = psm( ct ) + 1
474  130 continue
475 *
476 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
477 * and Q2 respectively. The eigenvalues/vectors which were not
478 * deflated go into the first K slots of DLAMDA and Q2 respectively,
479 * while those which were deflated go into the last N - K slots.
480 *
481  i = 1
482  iq1 = 1
483  iq2 = 1 + ( ctot( 1 )+ctot( 2 ) )*n1
484  DO 140 j = 1, ctot( 1 )
485  js = indx( i )
486  CALL dcopy( n1, q( 1, js ), 1, q2( iq1 ), 1 )
487  z( i ) = d( js )
488  i = i + 1
489  iq1 = iq1 + n1
490  140 continue
491 *
492  DO 150 j = 1, ctot( 2 )
493  js = indx( i )
494  CALL dcopy( n1, q( 1, js ), 1, q2( iq1 ), 1 )
495  CALL dcopy( n2, q( n1+1, js ), 1, q2( iq2 ), 1 )
496  z( i ) = d( js )
497  i = i + 1
498  iq1 = iq1 + n1
499  iq2 = iq2 + n2
500  150 continue
501 *
502  DO 160 j = 1, ctot( 3 )
503  js = indx( i )
504  CALL dcopy( n2, q( n1+1, js ), 1, q2( iq2 ), 1 )
505  z( i ) = d( js )
506  i = i + 1
507  iq2 = iq2 + n2
508  160 continue
509 *
510  iq1 = iq2
511  DO 170 j = 1, ctot( 4 )
512  js = indx( i )
513  CALL dcopy( n, q( 1, js ), 1, q2( iq2 ), 1 )
514  iq2 = iq2 + n
515  z( i ) = d( js )
516  i = i + 1
517  170 continue
518 *
519 * The deflated eigenvalues and their corresponding vectors go back
520 * into the last N - K slots of D and Q respectively.
521 *
522  IF( k.LT.n ) THEN
523  CALL dlacpy( 'A', n, ctot( 4 ), q2( iq1 ), n,
524  $ q( 1, k+1 ), ldq )
525  CALL dcopy( n-k, z( k+1 ), 1, d( k+1 ), 1 )
526  END IF
527 *
528 * Copy CTOT into COLTYP for referencing in DLAED3.
529 *
530  DO 180 j = 1, 4
531  coltyp( j ) = ctot( j )
532  180 continue
533 *
534  190 continue
535  return
536 *
537 * End of DLAED2
538 *
539  END