001:       SUBROUTINE SLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            I
009:       REAL               DSIGMA, RHO
010: *     ..
011: *     .. Array Arguments ..
012:       REAL               D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  This subroutine computes the square root of the I-th eigenvalue
019: *  of a positive symmetric rank-one modification of a 2-by-2 diagonal
020: *  matrix
021: *
022: *             diag( D ) * diag( D ) +  RHO *  Z * transpose(Z) .
023: *
024: *  The diagonal entries in the array D are assumed to satisfy
025: *
026: *             0 <= D(i) < D(j)  for  i < j .
027: *
028: *  We also assume RHO > 0 and that the Euclidean norm of the vector
029: *  Z is one.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  I      (input) INTEGER
035: *         The index of the eigenvalue to be computed.  I = 1 or I = 2.
036: *
037: *  D      (input) REAL array, dimension (2)
038: *         The original eigenvalues.  We assume 0 <= D(1) < D(2).
039: *
040: *  Z      (input) REAL array, dimension (2)
041: *         The components of the updating vector.
042: *
043: *  DELTA  (output) REAL array, dimension (2)
044: *         Contains (D(j) - sigma_I) in its  j-th component.
045: *         The vector DELTA contains the information necessary
046: *         to construct the eigenvectors.
047: *
048: *  RHO    (input) REAL
049: *         The scalar in the symmetric updating formula.
050: *
051: *  DSIGMA (output) REAL
052: *         The computed sigma_I, the I-th updated eigenvalue.
053: *
054: *  WORK   (workspace) REAL array, dimension (2)
055: *         WORK contains (D(j) + sigma_I) in its  j-th component.
056: *
057: *  Further Details
058: *  ===============
059: *
060: *  Based on contributions by
061: *     Ren-Cang Li, Computer Science Division, University of California
062: *     at Berkeley, USA
063: *
064: *  =====================================================================
065: *
066: *     .. Parameters ..
067:       REAL               ZERO, ONE, TWO, THREE, FOUR
068:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
069:      $                   THREE = 3.0E+0, FOUR = 4.0E+0 )
070: *     ..
071: *     .. Local Scalars ..
072:       REAL               B, C, DEL, DELSQ, TAU, W
073: *     ..
074: *     .. Intrinsic Functions ..
075:       INTRINSIC          ABS, SQRT
076: *     ..
077: *     .. Executable Statements ..
078: *
079:       DEL = D( 2 ) - D( 1 )
080:       DELSQ = DEL*( D( 2 )+D( 1 ) )
081:       IF( I.EQ.1 ) THEN
082:          W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )-
083:      $       Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL
084:          IF( W.GT.ZERO ) THEN
085:             B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
086:             C = RHO*Z( 1 )*Z( 1 )*DELSQ
087: *
088: *           B > ZERO, always
089: *
090: *           The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 )
091: *
092:             TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) )
093: *
094: *           The following TAU is DSIGMA - D( 1 )
095: *
096:             TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) )
097:             DSIGMA = D( 1 ) + TAU
098:             DELTA( 1 ) = -TAU
099:             DELTA( 2 ) = DEL - TAU
100:             WORK( 1 ) = TWO*D( 1 ) + TAU
101:             WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 )
102: *           DELTA( 1 ) = -Z( 1 ) / TAU
103: *           DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )
104:          ELSE
105:             B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
106:             C = RHO*Z( 2 )*Z( 2 )*DELSQ
107: *
108: *           The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
109: *
110:             IF( B.GT.ZERO ) THEN
111:                TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) )
112:             ELSE
113:                TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO
114:             END IF
115: *
116: *           The following TAU is DSIGMA - D( 2 )
117: *
118:             TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) )
119:             DSIGMA = D( 2 ) + TAU
120:             DELTA( 1 ) = -( DEL+TAU )
121:             DELTA( 2 ) = -TAU
122:             WORK( 1 ) = D( 1 ) + TAU + D( 2 )
123:             WORK( 2 ) = TWO*D( 2 ) + TAU
124: *           DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
125: *           DELTA( 2 ) = -Z( 2 ) / TAU
126:          END IF
127: *        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
128: *        DELTA( 1 ) = DELTA( 1 ) / TEMP
129: *        DELTA( 2 ) = DELTA( 2 ) / TEMP
130:       ELSE
131: *
132: *        Now I=2
133: *
134:          B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
135:          C = RHO*Z( 2 )*Z( 2 )*DELSQ
136: *
137: *        The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
138: *
139:          IF( B.GT.ZERO ) THEN
140:             TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO
141:          ELSE
142:             TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) )
143:          END IF
144: *
145: *        The following TAU is DSIGMA - D( 2 )
146: *
147:          TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) )
148:          DSIGMA = D( 2 ) + TAU
149:          DELTA( 1 ) = -( DEL+TAU )
150:          DELTA( 2 ) = -TAU
151:          WORK( 1 ) = D( 1 ) + TAU + D( 2 )
152:          WORK( 2 ) = TWO*D( 2 ) + TAU
153: *        DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
154: *        DELTA( 2 ) = -Z( 2 ) / TAU
155: *        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
156: *        DELTA( 1 ) = DELTA( 1 ) / TEMP
157: *        DELTA( 2 ) = DELTA( 2 ) / TEMP
158:       END IF
159:       RETURN
160: *
161: *     End of SLASD5
162: *
163:       END
164: