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

◆ dlaed2()

subroutine dlaed2 ( integer k,
integer n,
integer n1,
double precision, dimension( * ) d,
double precision, dimension( ldq, * ) q,
integer ldq,
integer, dimension( * ) indxq,
double precision rho,
double precision, dimension( * ) z,
double precision, dimension( * ) dlambda,
double precision, dimension( * ) w,
double precision, dimension( * ) q2,
integer, dimension( * ) indx,
integer, dimension( * ) indxc,
integer, dimension( * ) indxp,
integer, dimension( * ) coltyp,
integer info )

DLAED2 used by DSTEDC. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal.

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

Purpose:
!>
!> DLAED2 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 entry 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
!>         The number of non-deflated eigenvalues, and the order of the
!>         related secular equation. 0 <= K <=N.
!>
[in]N
!>          N is INTEGER
!>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
!>
[in]N1
!>          N1 is INTEGER
!>         The location of the last eigenvalue in the leading sub-matrix.
!>         min(1,N) <= N1 <= N/2.
!>
[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]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
!>         On entry, Q contains the eigenvectors of two submatrices in
!>         the two square blocks with corners at (1,1), (N1,N1)
!>         and (N1+1, N1+1), (N,N).
!>         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]INDXQ
!>          INDXQ is INTEGER array, dimension (N)
!>         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 N1 added to their
!>         values. Destroyed on exit.
!>
[in,out]RHO
!>          RHO is DOUBLE PRECISION
!>         On entry, the off-diagonal element associated with the rank-1
!>         cut which originally split the two submatrices which are now
!>         being recombined.
!>         On exit, RHO has been modified to the value required by
!>         DLAED3.
!>
[in]Z
!>          Z is DOUBLE PRECISION array, dimension (N)
!>         On entry, Z contains the updating vector (the last
!>         row of the first sub-eigenvector matrix and the first row of
!>         the second sub-eigenvector matrix).
!>         On exit, the contents of Z have been destroyed by the updating
!>         process.
!>
[out]DLAMBDA
!>          DLAMBDA is DOUBLE PRECISION array, dimension (N)
!>         A copy of the first K eigenvalues which will be used by
!>         DLAED3 to form the secular equation.
!>
[out]W
!>          W is DOUBLE PRECISION array, dimension (N)
!>         The first k values of the final deflation-altered z-vector
!>         which will be passed to DLAED3.
!>
[out]Q2
!>          Q2 is DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2)
!>         A copy of the first K eigenvectors which will be used by
!>         DLAED3 in a matrix multiply (DGEMM) to solve for the new
!>         eigenvectors.
!>
[out]INDX
!>          INDX is INTEGER array, dimension (N)
!>         The permutation used to sort the contents of DLAMBDA into
!>         ascending order.
!>
[out]INDXC
!>          INDXC is INTEGER array, dimension (N)
!>         The permutation used to arrange the columns of the deflated
!>         Q matrix into three groups:  the first group contains non-zero
!>         elements only at and above N1, the second contains
!>         non-zero elements only below N1, and the third is dense.
!>
[out]INDXP
!>          INDXP is INTEGER array, dimension (N)
!>         The permutation used to place deflated values of D at the end
!>         of the array.  INDXP(1:K) points to the nondeflated D-values
!>         and INDXP(K+1:N) points to the deflated eigenvalues.
!>
[out]COLTYP
!>          COLTYP is INTEGER array, dimension (N)
!>         During execution, a label which will indicate which of the
!>         following types a column in the Q2 matrix is:
!>         1 : non-zero in the upper half only;
!>         2 : dense;
!>         3 : non-zero in the lower half only;
!>         4 : deflated.
!>         On exit, COLTYP(i) is the number of columns of type i,
!>         for i=1 to 4 only.
!>
[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.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 208 of file dlaed2.f.

211*
212* -- LAPACK computational routine --
213* -- LAPACK is a software package provided by Univ. of Tennessee, --
214* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215*
216* .. Scalar Arguments ..
217 INTEGER INFO, K, LDQ, N, N1
218 DOUBLE PRECISION RHO
219* ..
220* .. Array Arguments ..
221 INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ),
222 $ INDXQ( * )
223 DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ),
224 $ W( * ), Z( * )
225* ..
226*
227* =====================================================================
228*
229* .. Parameters ..
230 DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
231 parameter( mone = -1.0d0, zero = 0.0d0, one = 1.0d0,
232 $ two = 2.0d0, eight = 8.0d0 )
233* ..
234* .. Local Arrays ..
235 INTEGER CTOT( 4 ), PSM( 4 )
236* ..
237* .. Local Scalars ..
238 INTEGER CT, I, IMAX, IQ1, IQ2, J, JMAX, JS, K2, N1P1,
239 $ N2, NJ, PJ
240 DOUBLE PRECISION C, EPS, S, T, TAU, TOL
241* ..
242* .. External Functions ..
243 INTEGER IDAMAX
244 DOUBLE PRECISION DLAMCH, DLAPY2
245 EXTERNAL idamax, dlamch, dlapy2
246* ..
247* .. External Subroutines ..
248 EXTERNAL dcopy, dlacpy, dlamrg, drot, dscal,
249 $ xerbla
250* ..
251* .. Intrinsic Functions ..
252 INTRINSIC abs, max, min, sqrt
253* ..
254* .. Executable Statements ..
255*
256* Test the input parameters.
257*
258 info = 0
259*
260 IF( n.LT.0 ) THEN
261 info = -2
262 ELSE IF( ldq.LT.max( 1, n ) ) THEN
263 info = -6
264 ELSE IF( min( 1, ( n / 2 ) ).GT.n1 .OR. ( n / 2 ).LT.n1 ) THEN
265 info = -3
266 END IF
267 IF( info.NE.0 ) THEN
268 CALL xerbla( 'DLAED2', -info )
269 RETURN
270 END IF
271*
272* Quick return if possible
273*
274 IF( n.EQ.0 )
275 $ RETURN
276*
277 n2 = n - n1
278 n1p1 = n1 + 1
279*
280 IF( rho.LT.zero ) THEN
281 CALL dscal( n2, mone, z( n1p1 ), 1 )
282 END IF
283*
284* Normalize z so that norm(z) = 1. Since z is the concatenation of
285* two normalized vectors, norm2(z) = sqrt(2).
286*
287 t = one / sqrt( two )
288 CALL dscal( n, t, z, 1 )
289*
290* RHO = ABS( norm(z)**2 * RHO )
291*
292 rho = abs( two*rho )
293*
294* Sort the eigenvalues into increasing order
295*
296 DO 10 i = n1p1, n
297 indxq( i ) = indxq( i ) + n1
298 10 CONTINUE
299*
300* re-integrate the deflated parts from the last pass
301*
302 DO 20 i = 1, n
303 dlambda( i ) = d( indxq( i ) )
304 20 CONTINUE
305 CALL dlamrg( n1, n2, dlambda, 1, 1, indxc )
306 DO 30 i = 1, n
307 indx( i ) = indxq( indxc( i ) )
308 30 CONTINUE
309*
310* Calculate the allowable deflation tolerance
311*
312 imax = idamax( n, z, 1 )
313 jmax = idamax( n, d, 1 )
314 eps = dlamch( 'Epsilon' )
315 tol = eight*eps*max( abs( d( jmax ) ), abs( z( imax ) ) )
316*
317* If the rank-1 modifier is small enough, no more needs to be done
318* except to reorganize Q so that its columns correspond with the
319* elements in D.
320*
321 IF( rho*abs( z( imax ) ).LE.tol ) THEN
322 k = 0
323 iq2 = 1
324 DO 40 j = 1, n
325 i = indx( j )
326 CALL dcopy( n, q( 1, i ), 1, q2( iq2 ), 1 )
327 dlambda( j ) = d( i )
328 iq2 = iq2 + n
329 40 CONTINUE
330 CALL dlacpy( 'A', n, n, q2, n, q, ldq )
331 CALL dcopy( n, dlambda, 1, d, 1 )
332 GO TO 190
333 END IF
334*
335* If there are multiple eigenvalues then the problem deflates. Here
336* the number of equal eigenvalues are found. As each equal
337* eigenvalue is found, an elementary reflector is computed to rotate
338* the corresponding eigensubspace so that the corresponding
339* components of Z are zero in this new basis.
340*
341 DO 50 i = 1, n1
342 coltyp( i ) = 1
343 50 CONTINUE
344 DO 60 i = n1p1, n
345 coltyp( i ) = 3
346 60 CONTINUE
347*
348*
349 k = 0
350 k2 = n + 1
351 DO 70 j = 1, n
352 nj = indx( j )
353 IF( rho*abs( z( nj ) ).LE.tol ) THEN
354*
355* Deflate due to small z component.
356*
357 k2 = k2 - 1
358 coltyp( nj ) = 4
359 indxp( k2 ) = nj
360 IF( j.EQ.n )
361 $ GO TO 100
362 ELSE
363 pj = nj
364 GO TO 80
365 END IF
366 70 CONTINUE
367 80 CONTINUE
368 j = j + 1
369 nj = indx( j )
370 IF( j.GT.n )
371 $ GO TO 100
372 IF( rho*abs( z( nj ) ).LE.tol ) THEN
373*
374* Deflate due to small z component.
375*
376 k2 = k2 - 1
377 coltyp( nj ) = 4
378 indxp( k2 ) = nj
379 ELSE
380*
381* Check if eigenvalues are close enough to allow deflation.
382*
383 s = z( pj )
384 c = z( nj )
385*
386* Find sqrt(a**2+b**2) without overflow or
387* destructive underflow.
388*
389 tau = dlapy2( c, s )
390 t = d( nj ) - d( pj )
391 c = c / tau
392 s = -s / tau
393 IF( abs( t*c*s ).LE.tol ) THEN
394*
395* Deflation is possible.
396*
397 z( nj ) = tau
398 z( pj ) = zero
399 IF( coltyp( nj ).NE.coltyp( pj ) )
400 $ coltyp( nj ) = 2
401 coltyp( pj ) = 4
402 CALL drot( n, q( 1, pj ), 1, q( 1, nj ), 1, c, s )
403 t = d( pj )*c**2 + d( nj )*s**2
404 d( nj ) = d( pj )*s**2 + d( nj )*c**2
405 d( pj ) = t
406 k2 = k2 - 1
407 i = 1
408 90 CONTINUE
409 IF( k2+i.LE.n ) THEN
410 IF( d( pj ).LT.d( indxp( k2+i ) ) ) THEN
411 indxp( k2+i-1 ) = indxp( k2+i )
412 indxp( k2+i ) = pj
413 i = i + 1
414 GO TO 90
415 ELSE
416 indxp( k2+i-1 ) = pj
417 END IF
418 ELSE
419 indxp( k2+i-1 ) = pj
420 END IF
421 pj = nj
422 ELSE
423 k = k + 1
424 dlambda( k ) = d( pj )
425 w( k ) = z( pj )
426 indxp( k ) = pj
427 pj = nj
428 END IF
429 END IF
430 GO TO 80
431 100 CONTINUE
432*
433* Record the last eigenvalue.
434*
435 k = k + 1
436 dlambda( k ) = d( pj )
437 w( k ) = z( pj )
438 indxp( k ) = pj
439*
440* Count up the total number of the various types of columns, then
441* form a permutation which positions the four column types into
442* four uniform groups (although one or more of these groups may be
443* empty).
444*
445 DO 110 j = 1, 4
446 ctot( j ) = 0
447 110 CONTINUE
448 DO 120 j = 1, n
449 ct = coltyp( j )
450 ctot( ct ) = ctot( ct ) + 1
451 120 CONTINUE
452*
453* PSM(*) = Position in SubMatrix (of types 1 through 4)
454*
455 psm( 1 ) = 1
456 psm( 2 ) = 1 + ctot( 1 )
457 psm( 3 ) = psm( 2 ) + ctot( 2 )
458 psm( 4 ) = psm( 3 ) + ctot( 3 )
459 k = n - ctot( 4 )
460*
461* Fill out the INDXC array so that the permutation which it induces
462* will place all type-1 columns first, all type-2 columns next,
463* then all type-3's, and finally all type-4's.
464*
465 DO 130 j = 1, n
466 js = indxp( j )
467 ct = coltyp( js )
468 indx( psm( ct ) ) = js
469 indxc( psm( ct ) ) = j
470 psm( ct ) = psm( ct ) + 1
471 130 CONTINUE
472*
473* Sort the eigenvalues and corresponding eigenvectors into DLAMBDA
474* and Q2 respectively. The eigenvalues/vectors which were not
475* deflated go into the first K slots of DLAMBDA and Q2 respectively,
476* while those which were deflated go into the last N - K slots.
477*
478 i = 1
479 iq1 = 1
480 iq2 = 1 + ( ctot( 1 )+ctot( 2 ) )*n1
481 DO 140 j = 1, ctot( 1 )
482 js = indx( i )
483 CALL dcopy( n1, q( 1, js ), 1, q2( iq1 ), 1 )
484 z( i ) = d( js )
485 i = i + 1
486 iq1 = iq1 + n1
487 140 CONTINUE
488*
489 DO 150 j = 1, ctot( 2 )
490 js = indx( i )
491 CALL dcopy( n1, q( 1, js ), 1, q2( iq1 ), 1 )
492 CALL dcopy( n2, q( n1+1, js ), 1, q2( iq2 ), 1 )
493 z( i ) = d( js )
494 i = i + 1
495 iq1 = iq1 + n1
496 iq2 = iq2 + n2
497 150 CONTINUE
498*
499 DO 160 j = 1, ctot( 3 )
500 js = indx( i )
501 CALL dcopy( n2, q( n1+1, js ), 1, q2( iq2 ), 1 )
502 z( i ) = d( js )
503 i = i + 1
504 iq2 = iq2 + n2
505 160 CONTINUE
506*
507 iq1 = iq2
508 DO 170 j = 1, ctot( 4 )
509 js = indx( i )
510 CALL dcopy( n, q( 1, js ), 1, q2( iq2 ), 1 )
511 iq2 = iq2 + n
512 z( i ) = d( js )
513 i = i + 1
514 170 CONTINUE
515*
516* The deflated eigenvalues and their corresponding vectors go back
517* into the last N - K slots of D and Q respectively.
518*
519 IF( k.LT.n ) THEN
520 CALL dlacpy( 'A', n, ctot( 4 ), q2( iq1 ), n,
521 $ q( 1, k+1 ), ldq )
522 CALL dcopy( n-k, z( k+1 ), 1, d( k+1 ), 1 )
523 END IF
524*
525* Copy CTOT into COLTYP for referencing in DLAED3.
526*
527 DO 180 j = 1, 4
528 coltyp( j ) = ctot( j )
529 180 CONTINUE
530*
531 190 CONTINUE
532 RETURN
533*
534* End of DLAED2
535*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.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
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
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: