3 /* Subroutine */ int dlaed5_(integer *i__, doublereal *d__, doublereal *z__,
4 doublereal *delta, doublereal *rho, doublereal *dlam)
6 /* System generated locals */
9 /* Builtin functions */
10 double sqrt(doublereal);
13 doublereal b, c__, w, del, tau, temp;
16 /* -- LAPACK routine (version 3.1) -- */
17 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
20 /* .. Scalar Arguments .. */
22 /* .. Array Arguments .. */
28 /* This subroutine computes the I-th eigenvalue of a symmetric rank-one */
29 /* modification of a 2-by-2 diagonal matrix */
31 /* diag( D ) + RHO * Z * transpose(Z) . */
33 /* The diagonal elements in the array D are assumed to satisfy */
35 /* D(i) < D(j) for i < j . */
37 /* We also assume RHO > 0 and that the Euclidean norm of the vector */
43 /* I (input) INTEGER */
44 /* The index of the eigenvalue to be computed. I = 1 or I = 2. */
46 /* D (input) DOUBLE PRECISION array, dimension (2) */
47 /* The original eigenvalues. We assume D(1) < D(2). */
49 /* Z (input) DOUBLE PRECISION array, dimension (2) */
50 /* The components of the updating vector. */
52 /* DELTA (output) DOUBLE PRECISION array, dimension (2) */
53 /* The vector DELTA contains the information necessary */
54 /* to construct the eigenvectors. */
56 /* RHO (input) DOUBLE PRECISION */
57 /* The scalar in the symmetric updating formula. */
59 /* DLAM (output) DOUBLE PRECISION */
60 /* The computed lambda_I, the I-th updated eigenvalue. */
65 /* Based on contributions by */
66 /* Ren-Cang Li, Computer Science Division, University of California */
67 /* at Berkeley, USA */
69 /* ===================================================================== */
71 /* .. Parameters .. */
73 /* .. Local Scalars .. */
75 /* .. Intrinsic Functions .. */
77 /* .. Executable Statements .. */
79 /* Parameter adjustments */
85 del = d__[2] - d__[1];
87 w = *rho * 2. * (z__[2] * z__[2] - z__[1] * z__[1]) / del + 1.;
89 b = del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
90 c__ = *rho * z__[1] * z__[1] * del;
92 /* B > ZERO, always */
94 tau = c__ * 2. / (b + sqrt((d__1 = b * b - c__ * 4., abs(d__1))));
96 delta[1] = -z__[1] / tau;
97 delta[2] = z__[2] / (del - tau);
99 b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
100 c__ = *rho * z__[2] * z__[2] * del;
102 tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
104 tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
106 *dlam = d__[2] + tau;
107 delta[1] = -z__[1] / (del + tau);
108 delta[2] = -z__[2] / tau;
110 temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
117 b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
118 c__ = *rho * z__[2] * z__[2] * del;
120 tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
122 tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
124 *dlam = d__[2] + tau;
125 delta[1] = -z__[1] / (del + tau);
126 delta[2] = -z__[2] / tau;
127 temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);