--- /dev/null
+#include "clapack.h"
+
+/* Subroutine */ int dlaed5_(integer *i__, doublereal *d__, doublereal *z__,
+ doublereal *delta, doublereal *rho, doublereal *dlam)
+{
+ /* System generated locals */
+ doublereal d__1;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ doublereal b, c__, w, del, tau, temp;
+
+
+/* -- LAPACK routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+/* .. Array Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* This subroutine computes the I-th eigenvalue of a symmetric rank-one */
+/* modification of a 2-by-2 diagonal matrix */
+
+/* diag( D ) + RHO * Z * transpose(Z) . */
+
+/* The diagonal elements in the array D are assumed to satisfy */
+
+/* D(i) < D(j) for i < j . */
+
+/* We also assume RHO > 0 and that the Euclidean norm of the vector */
+/* Z is one. */
+
+/* Arguments */
+/* ========= */
+
+/* I (input) INTEGER */
+/* The index of the eigenvalue to be computed. I = 1 or I = 2. */
+
+/* D (input) DOUBLE PRECISION array, dimension (2) */
+/* The original eigenvalues. We assume D(1) < D(2). */
+
+/* Z (input) DOUBLE PRECISION array, dimension (2) */
+/* The components of the updating vector. */
+
+/* DELTA (output) DOUBLE PRECISION array, dimension (2) */
+/* The vector DELTA contains the information necessary */
+/* to construct the eigenvectors. */
+
+/* RHO (input) DOUBLE PRECISION */
+/* The scalar in the symmetric updating formula. */
+
+/* DLAM (output) DOUBLE PRECISION */
+/* The computed lambda_I, the I-th updated eigenvalue. */
+
+/* Further Details */
+/* =============== */
+
+/* Based on contributions by */
+/* Ren-Cang Li, Computer Science Division, University of California */
+/* at Berkeley, USA */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+ /* Parameter adjustments */
+ --delta;
+ --z__;
+ --d__;
+
+ /* Function Body */
+ del = d__[2] - d__[1];
+ if (*i__ == 1) {
+ w = *rho * 2. * (z__[2] * z__[2] - z__[1] * z__[1]) / del + 1.;
+ if (w > 0.) {
+ b = del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
+ c__ = *rho * z__[1] * z__[1] * del;
+
+/* B > ZERO, always */
+
+ tau = c__ * 2. / (b + sqrt((d__1 = b * b - c__ * 4., abs(d__1))));
+ *dlam = d__[1] + tau;
+ delta[1] = -z__[1] / tau;
+ delta[2] = z__[2] / (del - tau);
+ } else {
+ b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
+ c__ = *rho * z__[2] * z__[2] * del;
+ if (b > 0.) {
+ tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
+ } else {
+ tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
+ }
+ *dlam = d__[2] + tau;
+ delta[1] = -z__[1] / (del + tau);
+ delta[2] = -z__[2] / tau;
+ }
+ temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
+ delta[1] /= temp;
+ delta[2] /= temp;
+ } else {
+
+/* Now I=2 */
+
+ b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
+ c__ = *rho * z__[2] * z__[2] * del;
+ if (b > 0.) {
+ tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
+ } else {
+ tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
+ }
+ *dlam = d__[2] + tau;
+ delta[1] = -z__[1] / (del + tau);
+ delta[2] = -z__[2] / tau;
+ temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
+ delta[1] /= temp;
+ delta[2] /= temp;
+ }
+ return 0;
+
+/* End OF DLAED5 */
+
+} /* dlaed5_ */