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

Definition at line 210 of file dlaed2.f.

212 *
213 * -- LAPACK computational routine --
214 * -- LAPACK is a software package provided by Univ. of Tennessee, --
215 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216 *
217 * .. Scalar Arguments ..
218  INTEGER INFO, K, LDQ, N, N1
219  DOUBLE PRECISION RHO
220 * ..
221 * .. Array Arguments ..
222  INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ),
223  $ INDXQ( * )
224  DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
225  $ W( * ), Z( * )
226 * ..
227 *
228 * =====================================================================
229 *
230 * .. Parameters ..
231  DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
232  parameter( mone = -1.0d0, zero = 0.0d0, one = 1.0d0,
233  $ two = 2.0d0, eight = 8.0d0 )
234 * ..
235 * .. Local Arrays ..
236  INTEGER CTOT( 4 ), PSM( 4 )
237 * ..
238 * .. Local Scalars ..
239  INTEGER CT, I, IMAX, IQ1, IQ2, J, JMAX, JS, K2, N1P1,
240  $ N2, NJ, PJ
241  DOUBLE PRECISION C, EPS, S, T, TAU, TOL
242 * ..
243 * .. External Functions ..
244  INTEGER IDAMAX
245  DOUBLE PRECISION DLAMCH, DLAPY2
246  EXTERNAL idamax, dlamch, dlapy2
247 * ..
248 * .. External Subroutines ..
249  EXTERNAL dcopy, dlacpy, dlamrg, drot, dscal, 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  dlamda( i ) = d( indxq( i ) )
304  20 CONTINUE
305  CALL dlamrg( n1, n2, dlamda, 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  dlamda( j ) = d( i )
328  iq2 = iq2 + n
329  40 CONTINUE
330  CALL dlacpy( 'A', n, n, q2, n, q, ldq )
331  CALL dcopy( n, dlamda, 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  dlamda( 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  dlamda( 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 DLAMDA
474 * and Q2 respectively. The eigenvalues/vectors which were not
475 * deflated go into the first K slots of DLAMDA 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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 dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
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
Here is the call graph for this function:
Here is the caller graph for this function: