LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zlaed8.f
Go to the documentation of this file.
1 *> \brief \b ZLAED8 used by sstedc. Merges eigenvalues and deflates secular equation. 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 ZLAED8 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed8.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed8.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed8.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
22 * Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
23 * GIVCOL, GIVNUM, INFO )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
27 * DOUBLE PRECISION RHO
28 * ..
29 * .. Array Arguments ..
30 * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
31 * $ INDXQ( * ), PERM( * )
32 * DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
33 * $ Z( * )
34 * COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZLAED8 merges the two sets of eigenvalues together into a single
44 *> sorted set. Then it tries to deflate the size of the problem.
45 *> There are two ways in which deflation can occur: when two or more
46 *> eigenvalues are close together or if there is a tiny element in the
47 *> Z vector. For each such occurrence the order of the related secular
48 *> equation problem is reduced by one.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[out] K
55 *> \verbatim
56 *> K is INTEGER
57 *> Contains the number of non-deflated eigenvalues.
58 *> This is the order of the related secular equation.
59 *> \endverbatim
60 *>
61 *> \param[in] N
62 *> \verbatim
63 *> N is INTEGER
64 *> The dimension of the symmetric tridiagonal matrix. N >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] QSIZ
68 *> \verbatim
69 *> QSIZ is INTEGER
70 *> The dimension of the unitary matrix used to reduce
71 *> the dense or band matrix to tridiagonal form.
72 *> QSIZ >= N if ICOMPQ = 1.
73 *> \endverbatim
74 *>
75 *> \param[in,out] Q
76 *> \verbatim
77 *> Q is COMPLEX*16 array, dimension (LDQ,N)
78 *> On entry, Q contains the eigenvectors of the partially solved
79 *> system which has been previously updated in matrix
80 *> multiplies with other partially solved eigensystems.
81 *> On exit, Q contains the trailing (N-K) updated eigenvectors
82 *> (those which were deflated) in its last N-K columns.
83 *> \endverbatim
84 *>
85 *> \param[in] LDQ
86 *> \verbatim
87 *> LDQ is INTEGER
88 *> The leading dimension of the array Q. LDQ >= max( 1, N ).
89 *> \endverbatim
90 *>
91 *> \param[in,out] D
92 *> \verbatim
93 *> D is DOUBLE PRECISION array, dimension (N)
94 *> On entry, D contains the eigenvalues of the two submatrices to
95 *> be combined. On exit, D contains the trailing (N-K) updated
96 *> eigenvalues (those which were deflated) sorted into increasing
97 *> order.
98 *> \endverbatim
99 *>
100 *> \param[in,out] RHO
101 *> \verbatim
102 *> RHO is DOUBLE PRECISION
103 *> Contains the off diagonal element associated with the rank-1
104 *> cut which originally split the two submatrices which are now
105 *> being recombined. RHO is modified during the computation to
106 *> the value required by DLAED3.
107 *> \endverbatim
108 *>
109 *> \param[in] CUTPNT
110 *> \verbatim
111 *> CUTPNT is INTEGER
112 *> Contains the location of the last eigenvalue in the leading
113 *> sub-matrix. MIN(1,N) <= CUTPNT <= N.
114 *> \endverbatim
115 *>
116 *> \param[in] Z
117 *> \verbatim
118 *> Z is DOUBLE PRECISION array, dimension (N)
119 *> On input this vector 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). The contents of Z are
122 *> destroyed during the updating process.
123 *> \endverbatim
124 *>
125 *> \param[out] DLAMDA
126 *> \verbatim
127 *> DLAMDA is DOUBLE PRECISION array, dimension (N)
128 *> Contains a copy of the first K eigenvalues which will be used
129 *> by DLAED3 to form the secular equation.
130 *> \endverbatim
131 *>
132 *> \param[out] Q2
133 *> \verbatim
134 *> Q2 is COMPLEX*16 array, dimension (LDQ2,N)
135 *> If ICOMPQ = 0, Q2 is not referenced. Otherwise,
136 *> Contains a copy of the first K eigenvectors which will be used
137 *> by DLAED7 in a matrix multiply (DGEMM) to update the new
138 *> eigenvectors.
139 *> \endverbatim
140 *>
141 *> \param[in] LDQ2
142 *> \verbatim
143 *> LDQ2 is INTEGER
144 *> The leading dimension of the array Q2. LDQ2 >= max( 1, N ).
145 *> \endverbatim
146 *>
147 *> \param[out] W
148 *> \verbatim
149 *> W is DOUBLE PRECISION array, dimension (N)
150 *> This will hold the first k values of the final
151 *> deflation-altered z-vector and will be passed to DLAED3.
152 *> \endverbatim
153 *>
154 *> \param[out] INDXP
155 *> \verbatim
156 *> INDXP is INTEGER array, dimension (N)
157 *> This will contain the permutation used to place deflated
158 *> values of D at the end of the array. On output INDXP(1:K)
159 *> points to the nondeflated D-values and INDXP(K+1:N)
160 *> points to the deflated eigenvalues.
161 *> \endverbatim
162 *>
163 *> \param[out] INDX
164 *> \verbatim
165 *> INDX is INTEGER array, dimension (N)
166 *> This will contain the permutation used to sort the contents of
167 *> D into ascending order.
168 *> \endverbatim
169 *>
170 *> \param[in] INDXQ
171 *> \verbatim
172 *> INDXQ is INTEGER array, dimension (N)
173 *> This contains the permutation which separately sorts the two
174 *> sub-problems in D into ascending order. Note that elements in
175 *> the second half of this permutation must first have CUTPNT
176 *> added to their values in order to be accurate.
177 *> \endverbatim
178 *>
179 *> \param[out] PERM
180 *> \verbatim
181 *> PERM is INTEGER array, dimension (N)
182 *> Contains the permutations (from deflation and sorting) to be
183 *> applied to each eigenblock.
184 *> \endverbatim
185 *>
186 *> \param[out] GIVPTR
187 *> \verbatim
188 *> GIVPTR is INTEGER
189 *> Contains the number of Givens rotations which took place in
190 *> this subproblem.
191 *> \endverbatim
192 *>
193 *> \param[out] GIVCOL
194 *> \verbatim
195 *> GIVCOL is INTEGER array, dimension (2, N)
196 *> Each pair of numbers indicates a pair of columns to take place
197 *> in a Givens rotation.
198 *> \endverbatim
199 *>
200 *> \param[out] GIVNUM
201 *> \verbatim
202 *> GIVNUM is DOUBLE PRECISION array, dimension (2, N)
203 *> Each number indicates the S value to be used in the
204 *> corresponding Givens rotation.
205 *> \endverbatim
206 *>
207 *> \param[out] INFO
208 *> \verbatim
209 *> INFO is INTEGER
210 *> = 0: successful exit.
211 *> < 0: if INFO = -i, the i-th argument had an illegal value.
212 *> \endverbatim
213 *
214 * Authors:
215 * ========
216 *
217 *> \author Univ. of Tennessee
218 *> \author Univ. of California Berkeley
219 *> \author Univ. of Colorado Denver
220 *> \author NAG Ltd.
221 *
222 *> \date September 2012
223 *
224 *> \ingroup complex16OTHERcomputational
225 *
226 * =====================================================================
227  SUBROUTINE zlaed8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
228  $ q2, ldq2, w, indxp, indx, indxq, perm, givptr,
229  $ givcol, givnum, info )
230 *
231 * -- LAPACK computational routine (version 3.4.2) --
232 * -- LAPACK is a software package provided by Univ. of Tennessee, --
233 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234 * September 2012
235 *
236 * .. Scalar Arguments ..
237  INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
238  DOUBLE PRECISION RHO
239 * ..
240 * .. Array Arguments ..
241  INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
242  $ indxq( * ), perm( * )
243  DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
244  $ z( * )
245  COMPLEX*16 Q( ldq, * ), Q2( ldq2, * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
252  parameter ( mone = -1.0d0, zero = 0.0d0, one = 1.0d0,
253  $ two = 2.0d0, eight = 8.0d0 )
254 * ..
255 * .. Local Scalars ..
256  INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
257  DOUBLE PRECISION C, EPS, S, T, TAU, TOL
258 * ..
259 * .. External Functions ..
260  INTEGER IDAMAX
261  DOUBLE PRECISION DLAMCH, DLAPY2
262  EXTERNAL idamax, dlamch, dlapy2
263 * ..
264 * .. External Subroutines ..
265  EXTERNAL dcopy, dlamrg, dscal, xerbla, zcopy, zdrot,
266  $ zlacpy
267 * ..
268 * .. Intrinsic Functions ..
269  INTRINSIC abs, max, min, sqrt
270 * ..
271 * .. Executable Statements ..
272 *
273 * Test the input parameters.
274 *
275  info = 0
276 *
277  IF( n.LT.0 ) THEN
278  info = -2
279  ELSE IF( qsiz.LT.n ) THEN
280  info = -3
281  ELSE IF( ldq.LT.max( 1, n ) ) THEN
282  info = -5
283  ELSE IF( cutpnt.LT.min( 1, n ) .OR. cutpnt.GT.n ) THEN
284  info = -8
285  ELSE IF( ldq2.LT.max( 1, n ) ) THEN
286  info = -12
287  END IF
288  IF( info.NE.0 ) THEN
289  CALL xerbla( 'ZLAED8', -info )
290  RETURN
291  END IF
292 *
293 * Need to initialize GIVPTR to O here in case of quick exit
294 * to prevent an unspecified code behavior (usually sigfault)
295 * when IWORK array on entry to *stedc is not zeroed
296 * (or at least some IWORK entries which used in *laed7 for GIVPTR).
297 *
298  givptr = 0
299 *
300 * Quick return if possible
301 *
302  IF( n.EQ.0 )
303  $ RETURN
304 *
305  n1 = cutpnt
306  n2 = n - n1
307  n1p1 = n1 + 1
308 *
309  IF( rho.LT.zero ) THEN
310  CALL dscal( n2, mone, z( n1p1 ), 1 )
311  END IF
312 *
313 * Normalize z so that norm(z) = 1
314 *
315  t = one / sqrt( two )
316  DO 10 j = 1, n
317  indx( j ) = j
318  10 CONTINUE
319  CALL dscal( n, t, z, 1 )
320  rho = abs( two*rho )
321 *
322 * Sort the eigenvalues into increasing order
323 *
324  DO 20 i = cutpnt + 1, n
325  indxq( i ) = indxq( i ) + cutpnt
326  20 CONTINUE
327  DO 30 i = 1, n
328  dlamda( i ) = d( indxq( i ) )
329  w( i ) = z( indxq( i ) )
330  30 CONTINUE
331  i = 1
332  j = cutpnt + 1
333  CALL dlamrg( n1, n2, dlamda, 1, 1, indx )
334  DO 40 i = 1, n
335  d( i ) = dlamda( indx( i ) )
336  z( i ) = w( indx( i ) )
337  40 CONTINUE
338 *
339 * Calculate the allowable deflation tolerance
340 *
341  imax = idamax( n, z, 1 )
342  jmax = idamax( n, d, 1 )
343  eps = dlamch( 'Epsilon' )
344  tol = eight*eps*abs( d( jmax ) )
345 *
346 * If the rank-1 modifier is small enough, no more needs to be done
347 * -- except to reorganize Q so that its columns correspond with the
348 * elements in D.
349 *
350  IF( rho*abs( z( imax ) ).LE.tol ) THEN
351  k = 0
352  DO 50 j = 1, n
353  perm( j ) = indxq( indx( j ) )
354  CALL zcopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
355  50 CONTINUE
356  CALL zlacpy( 'A', qsiz, n, q2( 1, 1 ), ldq2, q( 1, 1 ), ldq )
357  RETURN
358  END IF
359 *
360 * If there are multiple eigenvalues then the problem deflates. Here
361 * the number of equal eigenvalues are found. As each equal
362 * eigenvalue is found, an elementary reflector is computed to rotate
363 * the corresponding eigensubspace so that the corresponding
364 * components of Z are zero in this new basis.
365 *
366  k = 0
367  k2 = n + 1
368  DO 60 j = 1, n
369  IF( rho*abs( z( j ) ).LE.tol ) THEN
370 *
371 * Deflate due to small z component.
372 *
373  k2 = k2 - 1
374  indxp( k2 ) = j
375  IF( j.EQ.n )
376  $ GO TO 100
377  ELSE
378  jlam = j
379  GO TO 70
380  END IF
381  60 CONTINUE
382  70 CONTINUE
383  j = j + 1
384  IF( j.GT.n )
385  $ GO TO 90
386  IF( rho*abs( z( j ) ).LE.tol ) THEN
387 *
388 * Deflate due to small z component.
389 *
390  k2 = k2 - 1
391  indxp( k2 ) = j
392  ELSE
393 *
394 * Check if eigenvalues are close enough to allow deflation.
395 *
396  s = z( jlam )
397  c = z( j )
398 *
399 * Find sqrt(a**2+b**2) without overflow or
400 * destructive underflow.
401 *
402  tau = dlapy2( c, s )
403  t = d( j ) - d( jlam )
404  c = c / tau
405  s = -s / tau
406  IF( abs( t*c*s ).LE.tol ) THEN
407 *
408 * Deflation is possible.
409 *
410  z( j ) = tau
411  z( jlam ) = zero
412 *
413 * Record the appropriate Givens rotation
414 *
415  givptr = givptr + 1
416  givcol( 1, givptr ) = indxq( indx( jlam ) )
417  givcol( 2, givptr ) = indxq( indx( j ) )
418  givnum( 1, givptr ) = c
419  givnum( 2, givptr ) = s
420  CALL zdrot( qsiz, q( 1, indxq( indx( jlam ) ) ), 1,
421  $ q( 1, indxq( indx( j ) ) ), 1, c, s )
422  t = d( jlam )*c*c + d( j )*s*s
423  d( j ) = d( jlam )*s*s + d( j )*c*c
424  d( jlam ) = t
425  k2 = k2 - 1
426  i = 1
427  80 CONTINUE
428  IF( k2+i.LE.n ) THEN
429  IF( d( jlam ).LT.d( indxp( k2+i ) ) ) THEN
430  indxp( k2+i-1 ) = indxp( k2+i )
431  indxp( k2+i ) = jlam
432  i = i + 1
433  GO TO 80
434  ELSE
435  indxp( k2+i-1 ) = jlam
436  END IF
437  ELSE
438  indxp( k2+i-1 ) = jlam
439  END IF
440  jlam = j
441  ELSE
442  k = k + 1
443  w( k ) = z( jlam )
444  dlamda( k ) = d( jlam )
445  indxp( k ) = jlam
446  jlam = j
447  END IF
448  END IF
449  GO TO 70
450  90 CONTINUE
451 *
452 * Record the last eigenvalue.
453 *
454  k = k + 1
455  w( k ) = z( jlam )
456  dlamda( k ) = d( jlam )
457  indxp( k ) = jlam
458 *
459  100 CONTINUE
460 *
461 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
462 * and Q2 respectively. The eigenvalues/vectors which were not
463 * deflated go into the first K slots of DLAMDA and Q2 respectively,
464 * while those which were deflated go into the last N - K slots.
465 *
466  DO 110 j = 1, n
467  jp = indxp( j )
468  dlamda( j ) = d( jp )
469  perm( j ) = indxq( indx( jp ) )
470  CALL zcopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
471  110 CONTINUE
472 *
473 * The deflated eigenvalues and their corresponding vectors go back
474 * into the last N - K slots of D and Q respectively.
475 *
476  IF( k.LT.n ) THEN
477  CALL dcopy( n-k, dlamda( k+1 ), 1, d( k+1 ), 1 )
478  CALL zlacpy( 'A', qsiz, n-k, q2( 1, k+1 ), ldq2, q( 1, k+1 ),
479  $ ldq )
480  END IF
481 *
482  RETURN
483 *
484 * End of ZLAED8
485 *
486  END
subroutine zlaed8(K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, GIVCOL, GIVNUM, INFO)
ZLAED8 used by sstedc. Merges eigenvalues and deflates secular equation. Used when the original matri...
Definition: zlaed8.f:230
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
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:101
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:100
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55