LAPACK  3.10.0 LAPACK: Linear Algebra PACKage
dlaed2.f
Go to the documentation of this file.
1 *> \brief \b DLAED2 used by DSTEDC. 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 *> \ingroup auxOTHERcomputational
201 *
202 *> \par Contributors:
203 * ==================
204 *>
205 *> Jeff Rutter, Computer Science Division, University of California
206 *> at Berkeley, USA \n
207 *> Modified by Francoise Tisseur, University of Tennessee
208 *>
209 * =====================================================================
210  SUBROUTINE dlaed2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W,
211  \$ Q2, INDX, INDXC, INDXP, COLTYP, INFO )
212 *
213 * -- LAPACK computational routine --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 *
217 * .. Scalar Arguments ..
218  INTEGER INFO, K, LDQ, N, N1
219  DOUBLE PRECISION RHO
220 * ..
221 * .. Array Arguments ..
222  INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ),
223  \$ indxq( * )
224  DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
225  \$ w( * ), z( * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
232  parameter( mone = -1.0d0, zero = 0.0d0, one = 1.0d0,
233  \$ two = 2.0d0, eight = 8.0d0 )
234 * ..
235 * .. Local Arrays ..
236  INTEGER CTOT( 4 ), PSM( 4 )
237 * ..
238 * .. Local Scalars ..
239  INTEGER CT, I, IMAX, IQ1, IQ2, J, JMAX, JS, K2, N1P1,
240  \$ n2, nj, pj
241  DOUBLE PRECISION C, EPS, S, T, TAU, TOL
242 * ..
243 * .. External Functions ..
244  INTEGER IDAMAX
245  DOUBLE PRECISION DLAMCH, DLAPY2
246  EXTERNAL idamax, dlamch, dlapy2
247 * ..
248 * .. External Subroutines ..
249  EXTERNAL dcopy, dlacpy, dlamrg, drot, dscal, xerbla
250 * ..
251 * .. Intrinsic Functions ..
252  INTRINSIC abs, max, min, sqrt
253 * ..
254 * .. Executable Statements ..
255 *
256 * Test the input parameters.
257 *
258  info = 0
259 *
260  IF( n.LT.0 ) THEN
261  info = -2
262  ELSE IF( ldq.LT.max( 1, n ) ) THEN
263  info = -6
264  ELSE IF( min( 1, ( n / 2 ) ).GT.n1 .OR. ( n / 2 ).LT.n1 ) THEN
265  info = -3
266  END IF
267  IF( info.NE.0 ) THEN
268  CALL xerbla( 'DLAED2', -info )
269  RETURN
270  END IF
271 *
272 * Quick return if possible
273 *
274  IF( n.EQ.0 )
275  \$ RETURN
276 *
277  n2 = n - n1
278  n1p1 = n1 + 1
279 *
280  IF( rho.LT.zero ) THEN
281  CALL dscal( n2, mone, z( n1p1 ), 1 )
282  END IF
283 *
284 * Normalize z so that norm(z) = 1. Since z is the concatenation of
285 * two normalized vectors, norm2(z) = sqrt(2).
286 *
287  t = one / sqrt( two )
288  CALL dscal( n, t, z, 1 )
289 *
290 * RHO = ABS( norm(z)**2 * RHO )
291 *
292  rho = abs( two*rho )
293 *
294 * Sort the eigenvalues into increasing order
295 *
296  DO 10 i = n1p1, n
297  indxq( i ) = indxq( i ) + n1
298  10 CONTINUE
299 *
300 * re-integrate the deflated parts from the last pass
301 *
302  DO 20 i = 1, n
303  dlamda( i ) = d( indxq( i ) )
304  20 CONTINUE
305  CALL dlamrg( n1, n2, dlamda, 1, 1, indxc )
306  DO 30 i = 1, n
307  indx( i ) = indxq( indxc( i ) )
308  30 CONTINUE
309 *
310 * Calculate the allowable deflation tolerance
311 *
312  imax = idamax( n, z, 1 )
313  jmax = idamax( n, d, 1 )
314  eps = dlamch( 'Epsilon' )
315  tol = eight*eps*max( abs( d( jmax ) ), abs( z( imax ) ) )
316 *
317 * If the rank-1 modifier is small enough, no more needs to be done
318 * except to reorganize Q so that its columns correspond with the
319 * elements in D.
320 *
321  IF( rho*abs( z( imax ) ).LE.tol ) THEN
322  k = 0
323  iq2 = 1
324  DO 40 j = 1, n
325  i = indx( j )
326  CALL dcopy( n, q( 1, i ), 1, q2( iq2 ), 1 )
327  dlamda( j ) = d( i )
328  iq2 = iq2 + n
329  40 CONTINUE
330  CALL dlacpy( 'A', n, n, q2, n, q, ldq )
331  CALL dcopy( n, dlamda, 1, d, 1 )
332  GO TO 190
333  END IF
334 *
335 * If there are multiple eigenvalues then the problem deflates. Here
336 * the number of equal eigenvalues are found. As each equal
337 * eigenvalue is found, an elementary reflector is computed to rotate
338 * the corresponding eigensubspace so that the corresponding
339 * components of Z are zero in this new basis.
340 *
341  DO 50 i = 1, n1
342  coltyp( i ) = 1
343  50 CONTINUE
344  DO 60 i = n1p1, n
345  coltyp( i ) = 3
346  60 CONTINUE
347 *
348 *
349  k = 0
350  k2 = n + 1
351  DO 70 j = 1, n
352  nj = indx( j )
353  IF( rho*abs( z( nj ) ).LE.tol ) THEN
354 *
355 * Deflate due to small z component.
356 *
357  k2 = k2 - 1
358  coltyp( nj ) = 4
359  indxp( k2 ) = nj
360  IF( j.EQ.n )
361  \$ GO TO 100
362  ELSE
363  pj = nj
364  GO TO 80
365  END IF
366  70 CONTINUE
367  80 CONTINUE
368  j = j + 1
369  nj = indx( j )
370  IF( j.GT.n )
371  \$ GO TO 100
372  IF( rho*abs( z( nj ) ).LE.tol ) THEN
373 *
374 * Deflate due to small z component.
375 *
376  k2 = k2 - 1
377  coltyp( nj ) = 4
378  indxp( k2 ) = nj
379  ELSE
380 *
381 * Check if eigenvalues are close enough to allow deflation.
382 *
383  s = z( pj )
384  c = z( nj )
385 *
386 * Find sqrt(a**2+b**2) without overflow or
387 * destructive underflow.
388 *
389  tau = dlapy2( c, s )
390  t = d( nj ) - d( pj )
391  c = c / tau
392  s = -s / tau
393  IF( abs( t*c*s ).LE.tol ) THEN
394 *
395 * Deflation is possible.
396 *
397  z( nj ) = tau
398  z( pj ) = zero
399  IF( coltyp( nj ).NE.coltyp( pj ) )
400  \$ coltyp( nj ) = 2
401  coltyp( pj ) = 4
402  CALL drot( n, q( 1, pj ), 1, q( 1, nj ), 1, c, s )
403  t = d( pj )*c**2 + d( nj )*s**2
404  d( nj ) = d( pj )*s**2 + d( nj )*c**2
405  d( pj ) = t
406  k2 = k2 - 1
407  i = 1
408  90 CONTINUE
409  IF( k2+i.LE.n ) THEN
410  IF( d( pj ).LT.d( indxp( k2+i ) ) ) THEN
411  indxp( k2+i-1 ) = indxp( k2+i )
412  indxp( k2+i ) = pj
413  i = i + 1
414  GO TO 90
415  ELSE
416  indxp( k2+i-1 ) = pj
417  END IF
418  ELSE
419  indxp( k2+i-1 ) = pj
420  END IF
421  pj = nj
422  ELSE
423  k = k + 1
424  dlamda( k ) = d( pj )
425  w( k ) = z( pj )
426  indxp( k ) = pj
427  pj = nj
428  END IF
429  END IF
430  GO TO 80
431  100 CONTINUE
432 *
433 * Record the last eigenvalue.
434 *
435  k = k + 1
436  dlamda( k ) = d( pj )
437  w( k ) = z( pj )
438  indxp( k ) = pj
439 *
440 * Count up the total number of the various types of columns, then
441 * form a permutation which positions the four column types into
442 * four uniform groups (although one or more of these groups may be
443 * empty).
444 *
445  DO 110 j = 1, 4
446  ctot( j ) = 0
447  110 CONTINUE
448  DO 120 j = 1, n
449  ct = coltyp( j )
450  ctot( ct ) = ctot( ct ) + 1
451  120 CONTINUE
452 *
453 * PSM(*) = Position in SubMatrix (of types 1 through 4)
454 *
455  psm( 1 ) = 1
456  psm( 2 ) = 1 + ctot( 1 )
457  psm( 3 ) = psm( 2 ) + ctot( 2 )
458  psm( 4 ) = psm( 3 ) + ctot( 3 )
459  k = n - ctot( 4 )
460 *
461 * Fill out the INDXC array so that the permutation which it induces
462 * will place all type-1 columns first, all type-2 columns next,
463 * then all type-3's, and finally all type-4's.
464 *
465  DO 130 j = 1, n
466  js = indxp( j )
467  ct = coltyp( js )
468  indx( psm( ct ) ) = js
469  indxc( psm( ct ) ) = j
470  psm( ct ) = psm( ct ) + 1
471  130 CONTINUE
472 *
473 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
474 * and Q2 respectively. The eigenvalues/vectors which were not
475 * deflated go into the first K slots of DLAMDA and Q2 respectively,
476 * while those which were deflated go into the last N - K slots.
477 *
478  i = 1
479  iq1 = 1
480  iq2 = 1 + ( ctot( 1 )+ctot( 2 ) )*n1
481  DO 140 j = 1, ctot( 1 )
482  js = indx( i )
483  CALL dcopy( n1, q( 1, js ), 1, q2( iq1 ), 1 )
484  z( i ) = d( js )
485  i = i + 1
486  iq1 = iq1 + n1
487  140 CONTINUE
488 *
489  DO 150 j = 1, ctot( 2 )
490  js = indx( i )
491  CALL dcopy( n1, q( 1, js ), 1, q2( iq1 ), 1 )
492  CALL dcopy( n2, q( n1+1, js ), 1, q2( iq2 ), 1 )
493  z( i ) = d( js )
494  i = i + 1
495  iq1 = iq1 + n1
496  iq2 = iq2 + n2
497  150 CONTINUE
498 *
499  DO 160 j = 1, ctot( 3 )
500  js = indx( i )
501  CALL dcopy( n2, q( n1+1, js ), 1, q2( iq2 ), 1 )
502  z( i ) = d( js )
503  i = i + 1
504  iq2 = iq2 + n2
505  160 CONTINUE
506 *
507  iq1 = iq2
508  DO 170 j = 1, ctot( 4 )
509  js = indx( i )
510  CALL dcopy( n, q( 1, js ), 1, q2( iq2 ), 1 )
511  iq2 = iq2 + n
512  z( i ) = d( js )
513  i = i + 1
514  170 CONTINUE
515 *
516 * The deflated eigenvalues and their corresponding vectors go back
517 * into the last N - K slots of D and Q respectively.
518 *
519  IF( k.LT.n ) THEN
520  CALL dlacpy( 'A', n, ctot( 4 ), q2( iq1 ), n,
521  \$ q( 1, k+1 ), ldq )
522  CALL dcopy( n-k, z( k+1 ), 1, d( k+1 ), 1 )
523  END IF
524 *
525 * Copy CTOT into COLTYP for referencing in DLAED3.
526 *
527  DO 180 j = 1, 4
528  coltyp( j ) = ctot( j )
529  180 CONTINUE
530 *
531  190 CONTINUE
532  RETURN
533 *
534 * End of DLAED2
535 *
536  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlamrg(N1, N2, A, DTRD1, DTRD2, INDEX)
DLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: dlamrg.f:99
subroutine dlaed2(K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, Q2, INDX, INDXC, INDXP, COLTYP, INFO)
DLAED2 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matri...
Definition: dlaed2.f:212
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79