 LAPACK  3.9.0 LAPACK: Linear Algebra PACKage

## ◆ 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( * ) DLAMDA, 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 sstedc. Merges eigenvalues and deflates secular equation. Used when the original matrix is tridiagonal.

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] DLAMDA ``` DLAMDA 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 DLAMDA 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.```
Date
December 2016
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee

Definition at line 214 of file dlaed2.f.

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