LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
claed8.f
Go to the documentation of this file.
1 *> \brief \b CLAED8 used by CSTEDC. 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 CLAED8 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claed8.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claed8.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claed8.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAED8( 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 * REAL RHO
28 * ..
29 * .. Array Arguments ..
30 * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
31 * $ INDXQ( * ), PERM( * )
32 * REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
33 * $ Z( * )
34 * COMPLEX Q( LDQ, * ), Q2( LDQ2, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CLAED8 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 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 REAL 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 REAL
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 SLAED3.
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 REAL 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 REAL array, dimension (N)
128 *> Contains a copy of the first K eigenvalues which will be used
129 *> by SLAED3 to form the secular equation.
130 *> \endverbatim
131 *>
132 *> \param[out] Q2
133 *> \verbatim
134 *> Q2 is COMPLEX 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 SLAED7 in a matrix multiply (SGEMM) 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 REAL array, dimension (N)
150 *> This will hold the first k values of the final
151 *> deflation-altered z-vector and will be passed to SLAED3.
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 REAL 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 *> \ingroup complexOTHERcomputational
223 *
224 * =====================================================================
225  SUBROUTINE claed8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
226  $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
227  $ GIVCOL, GIVNUM, INFO )
228 *
229 * -- LAPACK computational routine --
230 * -- LAPACK is a software package provided by Univ. of Tennessee, --
231 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232 *
233 * .. Scalar Arguments ..
234  INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
235  REAL RHO
236 * ..
237 * .. Array Arguments ..
238  INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
239  $ INDXQ( * ), PERM( * )
240  REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
241  $ Z( * )
242  COMPLEX Q( LDQ, * ), Q2( LDQ2, * )
243 * ..
244 *
245 * =====================================================================
246 *
247 * .. Parameters ..
248  REAL MONE, ZERO, ONE, TWO, EIGHT
249  PARAMETER ( MONE = -1.0e0, zero = 0.0e0, one = 1.0e0,
250  $ two = 2.0e0, eight = 8.0e0 )
251 * ..
252 * .. Local Scalars ..
253  INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
254  REAL C, EPS, S, T, TAU, TOL
255 * ..
256 * .. External Functions ..
257  INTEGER ISAMAX
258  REAL SLAMCH, SLAPY2
259  EXTERNAL isamax, slamch, slapy2
260 * ..
261 * .. External Subroutines ..
262  EXTERNAL ccopy, clacpy, csrot, scopy, slamrg, sscal,
263  $ xerbla
264 * ..
265 * .. Intrinsic Functions ..
266  INTRINSIC abs, max, min, sqrt
267 * ..
268 * .. Executable Statements ..
269 *
270 * Test the input parameters.
271 *
272  info = 0
273 *
274  IF( n.LT.0 ) THEN
275  info = -2
276  ELSE IF( qsiz.LT.n ) THEN
277  info = -3
278  ELSE IF( ldq.LT.max( 1, n ) ) THEN
279  info = -5
280  ELSE IF( cutpnt.LT.min( 1, n ) .OR. cutpnt.GT.n ) THEN
281  info = -8
282  ELSE IF( ldq2.LT.max( 1, n ) ) THEN
283  info = -12
284  END IF
285  IF( info.NE.0 ) THEN
286  CALL xerbla( 'CLAED8', -info )
287  RETURN
288  END IF
289 *
290 * Need to initialize GIVPTR to O here in case of quick exit
291 * to prevent an unspecified code behavior (usually sigfault)
292 * when IWORK array on entry to *stedc is not zeroed
293 * (or at least some IWORK entries which used in *laed7 for GIVPTR).
294 *
295  givptr = 0
296 *
297 * Quick return if possible
298 *
299  IF( n.EQ.0 )
300  $ RETURN
301 *
302  n1 = cutpnt
303  n2 = n - n1
304  n1p1 = n1 + 1
305 *
306  IF( rho.LT.zero ) THEN
307  CALL sscal( n2, mone, z( n1p1 ), 1 )
308  END IF
309 *
310 * Normalize z so that norm(z) = 1
311 *
312  t = one / sqrt( two )
313  DO 10 j = 1, n
314  indx( j ) = j
315  10 CONTINUE
316  CALL sscal( n, t, z, 1 )
317  rho = abs( two*rho )
318 *
319 * Sort the eigenvalues into increasing order
320 *
321  DO 20 i = cutpnt + 1, n
322  indxq( i ) = indxq( i ) + cutpnt
323  20 CONTINUE
324  DO 30 i = 1, n
325  dlamda( i ) = d( indxq( i ) )
326  w( i ) = z( indxq( i ) )
327  30 CONTINUE
328  i = 1
329  j = cutpnt + 1
330  CALL slamrg( n1, n2, dlamda, 1, 1, indx )
331  DO 40 i = 1, n
332  d( i ) = dlamda( indx( i ) )
333  z( i ) = w( indx( i ) )
334  40 CONTINUE
335 *
336 * Calculate the allowable deflation tolerance
337 *
338  imax = isamax( n, z, 1 )
339  jmax = isamax( n, d, 1 )
340  eps = slamch( 'Epsilon' )
341  tol = eight*eps*abs( d( jmax ) )
342 *
343 * If the rank-1 modifier is small enough, no more needs to be done
344 * -- except to reorganize Q so that its columns correspond with the
345 * elements in D.
346 *
347  IF( rho*abs( z( imax ) ).LE.tol ) THEN
348  k = 0
349  DO 50 j = 1, n
350  perm( j ) = indxq( indx( j ) )
351  CALL ccopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
352  50 CONTINUE
353  CALL clacpy( 'A', qsiz, n, q2( 1, 1 ), ldq2, q( 1, 1 ), ldq )
354  RETURN
355  END IF
356 *
357 * If there are multiple eigenvalues then the problem deflates. Here
358 * the number of equal eigenvalues are found. As each equal
359 * eigenvalue is found, an elementary reflector is computed to rotate
360 * the corresponding eigensubspace so that the corresponding
361 * components of Z are zero in this new basis.
362 *
363  k = 0
364  k2 = n + 1
365  DO 60 j = 1, n
366  IF( rho*abs( z( j ) ).LE.tol ) THEN
367 *
368 * Deflate due to small z component.
369 *
370  k2 = k2 - 1
371  indxp( k2 ) = j
372  IF( j.EQ.n )
373  $ GO TO 100
374  ELSE
375  jlam = j
376  GO TO 70
377  END IF
378  60 CONTINUE
379  70 CONTINUE
380  j = j + 1
381  IF( j.GT.n )
382  $ GO TO 90
383  IF( rho*abs( z( j ) ).LE.tol ) THEN
384 *
385 * Deflate due to small z component.
386 *
387  k2 = k2 - 1
388  indxp( k2 ) = j
389  ELSE
390 *
391 * Check if eigenvalues are close enough to allow deflation.
392 *
393  s = z( jlam )
394  c = z( j )
395 *
396 * Find sqrt(a**2+b**2) without overflow or
397 * destructive underflow.
398 *
399  tau = slapy2( c, s )
400  t = d( j ) - d( jlam )
401  c = c / tau
402  s = -s / tau
403  IF( abs( t*c*s ).LE.tol ) THEN
404 *
405 * Deflation is possible.
406 *
407  z( j ) = tau
408  z( jlam ) = zero
409 *
410 * Record the appropriate Givens rotation
411 *
412  givptr = givptr + 1
413  givcol( 1, givptr ) = indxq( indx( jlam ) )
414  givcol( 2, givptr ) = indxq( indx( j ) )
415  givnum( 1, givptr ) = c
416  givnum( 2, givptr ) = s
417  CALL csrot( qsiz, q( 1, indxq( indx( jlam ) ) ), 1,
418  $ q( 1, indxq( indx( j ) ) ), 1, c, s )
419  t = d( jlam )*c*c + d( j )*s*s
420  d( j ) = d( jlam )*s*s + d( j )*c*c
421  d( jlam ) = t
422  k2 = k2 - 1
423  i = 1
424  80 CONTINUE
425  IF( k2+i.LE.n ) THEN
426  IF( d( jlam ).LT.d( indxp( k2+i ) ) ) THEN
427  indxp( k2+i-1 ) = indxp( k2+i )
428  indxp( k2+i ) = jlam
429  i = i + 1
430  GO TO 80
431  ELSE
432  indxp( k2+i-1 ) = jlam
433  END IF
434  ELSE
435  indxp( k2+i-1 ) = jlam
436  END IF
437  jlam = j
438  ELSE
439  k = k + 1
440  w( k ) = z( jlam )
441  dlamda( k ) = d( jlam )
442  indxp( k ) = jlam
443  jlam = j
444  END IF
445  END IF
446  GO TO 70
447  90 CONTINUE
448 *
449 * Record the last eigenvalue.
450 *
451  k = k + 1
452  w( k ) = z( jlam )
453  dlamda( k ) = d( jlam )
454  indxp( k ) = jlam
455 *
456  100 CONTINUE
457 *
458 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
459 * and Q2 respectively. The eigenvalues/vectors which were not
460 * deflated go into the first K slots of DLAMDA and Q2 respectively,
461 * while those which were deflated go into the last N - K slots.
462 *
463  DO 110 j = 1, n
464  jp = indxp( j )
465  dlamda( j ) = d( jp )
466  perm( j ) = indxq( indx( jp ) )
467  CALL ccopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
468  110 CONTINUE
469 *
470 * The deflated eigenvalues and their corresponding vectors go back
471 * into the last N - K slots of D and Q respectively.
472 *
473  IF( k.LT.n ) THEN
474  CALL scopy( n-k, dlamda( k+1 ), 1, d( k+1 ), 1 )
475  CALL clacpy( 'A', qsiz, n-k, q2( 1, k+1 ), ldq2, q( 1, k+1 ),
476  $ ldq )
477  END IF
478 *
479  RETURN
480 *
481 * End of CLAED8
482 *
483  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slamrg(N1, N2, A, STRD1, STRD2, INDEX)
SLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: slamrg.f:99
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
Definition: csrot.f:98
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine claed8(K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, GIVCOL, GIVNUM, INFO)
CLAED8 used by CSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matri...
Definition: claed8.f:228
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79