LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
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
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, DLAMBDA,
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( * ), DLAMBDA( * ), 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] DLAMBDA
126*> \verbatim
127*> DLAMBDA 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 laed8
223*
224* =====================================================================
225 SUBROUTINE claed8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMBDA,
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( * ), DLAMBDA( * ), 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 dlambda( 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, dlambda, 1, 1, indx )
331 DO 40 i = 1, n
332 d( i ) = dlambda( 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 dlambda( 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 dlambda( k ) = d( jlam )
454 indxp( k ) = jlam
455*
456 100 CONTINUE
457*
458* Sort the eigenvalues and corresponding eigenvectors into DLAMBDA
459* and Q2 respectively. The eigenvalues/vectors which were not
460* deflated go into the first K slots of DLAMBDA 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 dlambda( 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, dlambda( 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)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
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, dlambda, 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 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 csrot(n, cx, incx, cy, incy, c, s)
CSROT
Definition csrot.f:98
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79