LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.

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]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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
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.4.2) --
216 * -- LAPACK is a software package provided by Univ. of Tennessee, --
217 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218 * September 2012
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 *
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
idamax
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
dlapy2
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65

Here is the call graph for this function:

Here is the caller graph for this function: