LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlaed8()

subroutine zlaed8 ( integer k,
integer n,
integer qsiz,
complex*16, dimension( ldq, * ) q,
integer ldq,
double precision, dimension( * ) d,
double precision rho,
integer cutpnt,
double precision, dimension( * ) z,
double precision, dimension( * ) dlambda,
complex*16, dimension( ldq2, * ) q2,
integer ldq2,
double precision, dimension( * ) w,
integer, dimension( * ) indxp,
integer, dimension( * ) indx,
integer, dimension( * ) indxq,
integer, dimension( * ) perm,
integer givptr,
integer, dimension( 2, * ) givcol,
double precision, dimension( 2, * ) givnum,
integer info )

ZLAED8 used by ZSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is dense.

Download ZLAED8 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZLAED8 merges the two sets of eigenvalues together into a single
!> sorted set.  Then it tries to deflate the size of the problem.
!> There are two ways in which deflation can occur:  when two or more
!> eigenvalues are close together or if there is a tiny element in the
!> Z vector.  For each such occurrence the order of the related secular
!> equation problem is reduced by one.
!>
Parameters
[out]K
!>          K is INTEGER
!>         Contains the number of non-deflated eigenvalues.
!>         This is the order of the related secular equation.
!>
[in]N
!>          N is INTEGER
!>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
!>
[in]QSIZ
!>          QSIZ is INTEGER
!>         The dimension of the unitary matrix used to reduce
!>         the dense or band matrix to tridiagonal form.
!>         QSIZ >= N if ICOMPQ = 1.
!>
[in,out]Q
!>          Q is COMPLEX*16 array, dimension (LDQ,N)
!>         On entry, Q contains the eigenvectors of the partially solved
!>         system which has been previously updated in matrix
!>         multiplies with other partially solved eigensystems.
!>         On exit, Q contains the trailing (N-K) updated eigenvectors
!>         (those which were deflated) in its last N-K columns.
!>
[in]LDQ
!>          LDQ is INTEGER
!>         The leading dimension of the array Q.  LDQ >= max( 1, N ).
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>         On entry, D contains the eigenvalues of the two submatrices to
!>         be combined.  On exit, D contains the trailing (N-K) updated
!>         eigenvalues (those which were deflated) sorted into increasing
!>         order.
!>
[in,out]RHO
!>          RHO is DOUBLE PRECISION
!>         Contains the off diagonal element associated with the rank-1
!>         cut which originally split the two submatrices which are now
!>         being recombined. RHO is modified during the computation to
!>         the value required by DLAED3.
!>
[in]CUTPNT
!>          CUTPNT is INTEGER
!>         Contains the location of the last eigenvalue in the leading
!>         sub-matrix.  MIN(1,N) <= CUTPNT <= N.
!>
[in]Z
!>          Z is DOUBLE PRECISION array, dimension (N)
!>         On input this vector contains the updating vector (the last
!>         row of the first sub-eigenvector matrix and the first row of
!>         the second sub-eigenvector matrix).  The contents of Z are
!>         destroyed during the updating process.
!>
[out]DLAMBDA
!>          DLAMBDA is DOUBLE PRECISION array, dimension (N)
!>         Contains a copy of the first K eigenvalues which will be used
!>         by DLAED3 to form the secular equation.
!>
[out]Q2
!>          Q2 is COMPLEX*16 array, dimension (LDQ2,N)
!>         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
!>         Contains a copy of the first K eigenvectors which will be used
!>         by DLAED7 in a matrix multiply (DGEMM) to update the new
!>         eigenvectors.
!>
[in]LDQ2
!>          LDQ2 is INTEGER
!>         The leading dimension of the array Q2.  LDQ2 >= max( 1, N ).
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>         This will hold the first k values of the final
!>         deflation-altered z-vector and will be passed to DLAED3.
!>
[out]INDXP
!>          INDXP is INTEGER array, dimension (N)
!>         This will contain the permutation used to place deflated
!>         values of D at the end of the array. On output INDXP(1:K)
!>         points to the nondeflated D-values and INDXP(K+1:N)
!>         points to the deflated eigenvalues.
!>
[out]INDX
!>          INDX is INTEGER array, dimension (N)
!>         This will contain the permutation used to sort the contents of
!>         D into ascending order.
!>
[in]INDXQ
!>          INDXQ is INTEGER array, dimension (N)
!>         This contains the permutation which separately sorts the two
!>         sub-problems in D into ascending order.  Note that elements in
!>         the second half of this permutation must first have CUTPNT
!>         added to their values in order to be accurate.
!>
[out]PERM
!>          PERM is INTEGER array, dimension (N)
!>         Contains the permutations (from deflation and sorting) to be
!>         applied to each eigenblock.
!>
[out]GIVPTR
!>          GIVPTR is INTEGER
!>         Contains the number of Givens rotations which took place in
!>         this subproblem.
!>
[out]GIVCOL
!>          GIVCOL is INTEGER array, dimension (2, N)
!>         Each pair of numbers indicates a pair of columns to take place
!>         in a Givens rotation.
!>
[out]GIVNUM
!>          GIVNUM is DOUBLE PRECISION array, dimension (2, N)
!>         Each number indicates the S value to be used in the
!>         corresponding Givens rotation.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file zlaed8.f.

227*
228* -- LAPACK computational routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
234 DOUBLE PRECISION RHO
235* ..
236* .. Array Arguments ..
237 INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
238 $ INDXQ( * ), PERM( * )
239 DOUBLE PRECISION D( * ), DLAMBDA( * ), GIVNUM( 2, * ), W( * ),
240 $ Z( * )
241 COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
242* ..
243*
244* =====================================================================
245*
246* .. Parameters ..
247 DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
248 parameter( mone = -1.0d0, zero = 0.0d0, one = 1.0d0,
249 $ two = 2.0d0, eight = 8.0d0 )
250* ..
251* .. Local Scalars ..
252 INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
253 DOUBLE PRECISION C, EPS, S, T, TAU, TOL
254* ..
255* .. External Functions ..
256 INTEGER IDAMAX
257 DOUBLE PRECISION DLAMCH, DLAPY2
258 EXTERNAL idamax, dlamch, dlapy2
259* ..
260* .. External Subroutines ..
261 EXTERNAL dcopy, dlamrg, dscal, xerbla, zcopy,
262 $ zdrot,
263 $ zlacpy
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( 'ZLAED8', -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 dscal( 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 dscal( 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 dlamrg( 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 = idamax( n, z, 1 )
339 jmax = idamax( n, d, 1 )
340 eps = dlamch( '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 zcopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
352 50 CONTINUE
353 CALL zlacpy( 'A', qsiz, n, q2( 1, 1 ), ldq2, q( 1, 1 ),
354 $ ldq )
355 RETURN
356 END IF
357*
358* If there are multiple eigenvalues then the problem deflates. Here
359* the number of equal eigenvalues are found. As each equal
360* eigenvalue is found, an elementary reflector is computed to rotate
361* the corresponding eigensubspace so that the corresponding
362* components of Z are zero in this new basis.
363*
364 k = 0
365 k2 = n + 1
366 DO 60 j = 1, n
367 IF( rho*abs( z( j ) ).LE.tol ) THEN
368*
369* Deflate due to small z component.
370*
371 k2 = k2 - 1
372 indxp( k2 ) = j
373 IF( j.EQ.n )
374 $ GO TO 100
375 ELSE
376 jlam = j
377 GO TO 70
378 END IF
379 60 CONTINUE
380 70 CONTINUE
381 j = j + 1
382 IF( j.GT.n )
383 $ GO TO 90
384 IF( rho*abs( z( j ) ).LE.tol ) THEN
385*
386* Deflate due to small z component.
387*
388 k2 = k2 - 1
389 indxp( k2 ) = j
390 ELSE
391*
392* Check if eigenvalues are close enough to allow deflation.
393*
394 s = z( jlam )
395 c = z( j )
396*
397* Find sqrt(a**2+b**2) without overflow or
398* destructive underflow.
399*
400 tau = dlapy2( c, s )
401 t = d( j ) - d( jlam )
402 c = c / tau
403 s = -s / tau
404 IF( abs( t*c*s ).LE.tol ) THEN
405*
406* Deflation is possible.
407*
408 z( j ) = tau
409 z( jlam ) = zero
410*
411* Record the appropriate Givens rotation
412*
413 givptr = givptr + 1
414 givcol( 1, givptr ) = indxq( indx( jlam ) )
415 givcol( 2, givptr ) = indxq( indx( j ) )
416 givnum( 1, givptr ) = c
417 givnum( 2, givptr ) = s
418 CALL zdrot( qsiz, q( 1, indxq( indx( jlam ) ) ), 1,
419 $ q( 1, indxq( indx( j ) ) ), 1, c, s )
420 t = d( jlam )*c*c + d( j )*s*s
421 d( j ) = d( jlam )*s*s + d( j )*c*c
422 d( jlam ) = t
423 k2 = k2 - 1
424 i = 1
425 80 CONTINUE
426 IF( k2+i.LE.n ) THEN
427 IF( d( jlam ).LT.d( indxp( k2+i ) ) ) THEN
428 indxp( k2+i-1 ) = indxp( k2+i )
429 indxp( k2+i ) = jlam
430 i = i + 1
431 GO TO 80
432 ELSE
433 indxp( k2+i-1 ) = jlam
434 END IF
435 ELSE
436 indxp( k2+i-1 ) = jlam
437 END IF
438 jlam = j
439 ELSE
440 k = k + 1
441 w( k ) = z( jlam )
442 dlambda( k ) = d( jlam )
443 indxp( k ) = jlam
444 jlam = j
445 END IF
446 END IF
447 GO TO 70
448 90 CONTINUE
449*
450* Record the last eigenvalue.
451*
452 k = k + 1
453 w( k ) = z( jlam )
454 dlambda( k ) = d( jlam )
455 indxp( k ) = jlam
456*
457 100 CONTINUE
458*
459* Sort the eigenvalues and corresponding eigenvectors into DLAMBDA
460* and Q2 respectively. The eigenvalues/vectors which were not
461* deflated go into the first K slots of DLAMBDA and Q2 respectively,
462* while those which were deflated go into the last N - K slots.
463*
464 DO 110 j = 1, n
465 jp = indxp( j )
466 dlambda( j ) = d( jp )
467 perm( j ) = indxq( indx( jp ) )
468 CALL zcopy( qsiz, q( 1, perm( j ) ), 1, q2( 1, j ), 1 )
469 110 CONTINUE
470*
471* The deflated eigenvalues and their corresponding vectors go back
472* into the last N - K slots of D and Q respectively.
473*
474 IF( k.LT.n ) THEN
475 CALL dcopy( n-k, dlambda( k+1 ), 1, d( k+1 ), 1 )
476 CALL zlacpy( 'A', qsiz, n-k, q2( 1, k+1 ), ldq2, q( 1,
477 $ k+1 ),
478 $ ldq )
479 END IF
480*
481 RETURN
482*
483* End of ZLAED8
484*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlamrg
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:97
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:61
zdrot
subroutine zdrot(n, zx, incx, zy, incy, c, s)
ZDROT
Definition zdrot.f:98
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: