File:  [local] / rpl / lapack / lapack / dlaed5.f
Revision 1.5: download - view: text, annotated - select for diffs - revision graph
Sat Aug 7 13:22:16 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Mise à jour globale de Lapack 3.2.2.

    1:       SUBROUTINE DLAED5( I, D, Z, DELTA, RHO, DLAM )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            I
   10:       DOUBLE PRECISION   DLAM, RHO
   11: *     ..
   12: *     .. Array Arguments ..
   13:       DOUBLE PRECISION   D( 2 ), DELTA( 2 ), Z( 2 )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  This subroutine computes the I-th eigenvalue of a symmetric rank-one
   20: *  modification of a 2-by-2 diagonal matrix
   21: *
   22: *             diag( D )  +  RHO *  Z * transpose(Z) .
   23: *
   24: *  The diagonal elements in the array D are assumed to satisfy
   25: *
   26: *             D(i) < D(j)  for  i < j .
   27: *
   28: *  We also assume RHO > 0 and that the Euclidean norm of the vector
   29: *  Z is one.
   30: *
   31: *  Arguments
   32: *  =========
   33: *
   34: *  I      (input) INTEGER
   35: *         The index of the eigenvalue to be computed.  I = 1 or I = 2.
   36: *
   37: *  D      (input) DOUBLE PRECISION array, dimension (2)
   38: *         The original eigenvalues.  We assume D(1) < D(2).
   39: *
   40: *  Z      (input) DOUBLE PRECISION array, dimension (2)
   41: *         The components of the updating vector.
   42: *
   43: *  DELTA  (output) DOUBLE PRECISION array, dimension (2)
   44: *         The vector DELTA contains the information necessary
   45: *         to construct the eigenvectors.
   46: *
   47: *  RHO    (input) DOUBLE PRECISION
   48: *         The scalar in the symmetric updating formula.
   49: *
   50: *  DLAM   (output) DOUBLE PRECISION
   51: *         The computed lambda_I, the I-th updated eigenvalue.
   52: *
   53: *  Further Details
   54: *  ===============
   55: *
   56: *  Based on contributions by
   57: *     Ren-Cang Li, Computer Science Division, University of California
   58: *     at Berkeley, USA
   59: *
   60: *  =====================================================================
   61: *
   62: *     .. Parameters ..
   63:       DOUBLE PRECISION   ZERO, ONE, TWO, FOUR
   64:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
   65:      $                   FOUR = 4.0D0 )
   66: *     ..
   67: *     .. Local Scalars ..
   68:       DOUBLE PRECISION   B, C, DEL, TAU, TEMP, W
   69: *     ..
   70: *     .. Intrinsic Functions ..
   71:       INTRINSIC          ABS, SQRT
   72: *     ..
   73: *     .. Executable Statements ..
   74: *
   75:       DEL = D( 2 ) - D( 1 )
   76:       IF( I.EQ.1 ) THEN
   77:          W = ONE + TWO*RHO*( Z( 2 )*Z( 2 )-Z( 1 )*Z( 1 ) ) / DEL
   78:          IF( W.GT.ZERO ) THEN
   79:             B = DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
   80:             C = RHO*Z( 1 )*Z( 1 )*DEL
   81: *
   82: *           B > ZERO, always
   83: *
   84:             TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) )
   85:             DLAM = D( 1 ) + TAU
   86:             DELTA( 1 ) = -Z( 1 ) / TAU
   87:             DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )
   88:          ELSE
   89:             B = -DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
   90:             C = RHO*Z( 2 )*Z( 2 )*DEL
   91:             IF( B.GT.ZERO ) THEN
   92:                TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) )
   93:             ELSE
   94:                TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO
   95:             END IF
   96:             DLAM = D( 2 ) + TAU
   97:             DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
   98:             DELTA( 2 ) = -Z( 2 ) / TAU
   99:          END IF
  100:          TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
  101:          DELTA( 1 ) = DELTA( 1 ) / TEMP
  102:          DELTA( 2 ) = DELTA( 2 ) / TEMP
  103:       ELSE
  104: *
  105: *     Now I=2
  106: *
  107:          B = -DEL + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
  108:          C = RHO*Z( 2 )*Z( 2 )*DEL
  109:          IF( B.GT.ZERO ) THEN
  110:             TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO
  111:          ELSE
  112:             TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) )
  113:          END IF
  114:          DLAM = D( 2 ) + TAU
  115:          DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
  116:          DELTA( 2 ) = -Z( 2 ) / TAU
  117:          TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
  118:          DELTA( 1 ) = DELTA( 1 ) / TEMP
  119:          DELTA( 2 ) = DELTA( 2 ) / TEMP
  120:       END IF
  121:       RETURN
  122: *
  123: *     End OF DLAED5
  124: *
  125:       END

CVSweb interface <joel.bertrand@systella.fr>