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