dlasd5.c

Go to the documentation of this file.

/* dlasd5.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"
#include "blaswrap.h"

/* Subroutine */ int dlasd5_(integer *i__, doublereal *d__, doublereal *z__, 
        doublereal *delta, doublereal *rho, doublereal *dsigma, doublereal *
        work)
{
    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    doublereal b, c__, w, del, tau, delsq;


/*  -- LAPACK auxiliary routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  This subroutine computes the square root of the I-th eigenvalue */
/*  of a positive symmetric rank-one modification of a 2-by-2 diagonal */
/*  matrix */

/*             diag( D ) * diag( D ) +  RHO *  Z * transpose(Z) . */

/*  The diagonal entries in the array D are assumed to satisfy */

/*             0 <= 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 0 <= D(1) < D(2). */

/*  Z      (input) DOUBLE PRECISION array, dimension ( 2 ) */
/*         The components of the updating vector. */

/*  DELTA  (output) DOUBLE PRECISION array, dimension ( 2 ) */
/*         Contains (D(j) - sigma_I) in its  j-th component. */
/*         The vector DELTA contains the information necessary */
/*         to construct the eigenvectors. */

/*  RHO    (input) DOUBLE PRECISION */
/*         The scalar in the symmetric updating formula. */

/*  DSIGMA (output) DOUBLE PRECISION */
/*         The computed sigma_I, the I-th updated eigenvalue. */

/*  WORK   (workspace) DOUBLE PRECISION array, dimension ( 2 ) */
/*         WORK contains (D(j) + sigma_I) in its  j-th component. */

/*  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 */
    --work;
    --delta;
    --z__;
    --d__;

    /* Function Body */
    del = d__[2] - d__[1];
    delsq = del * (d__[2] + d__[1]);
    if (*i__ == 1) {
        w = *rho * 4. * (z__[2] * z__[2] / (d__[1] + d__[2] * 3.) - z__[1] * 
                z__[1] / (d__[1] * 3. + d__[2])) / del + 1.;
        if (w > 0.) {
            b = delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
            c__ = *rho * z__[1] * z__[1] * delsq;

/*           B > ZERO, always */

/*           The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 ) */

            tau = c__ * 2. / (b + sqrt((d__1 = b * b - c__ * 4., abs(d__1))));

/*           The following TAU is DSIGMA - D( 1 ) */

            tau /= d__[1] + sqrt(d__[1] * d__[1] + tau);
            *dsigma = d__[1] + tau;
            delta[1] = -tau;
            delta[2] = del - tau;
            work[1] = d__[1] * 2. + tau;
            work[2] = d__[1] + tau + d__[2];
/*           DELTA( 1 ) = -Z( 1 ) / TAU */
/*           DELTA( 2 ) = Z( 2 ) / ( DEL-TAU ) */
        } else {
            b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
            c__ = *rho * z__[2] * z__[2] * delsq;

/*           The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) */

            if (b > 0.) {
                tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
            } else {
                tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
            }

/*           The following TAU is DSIGMA - D( 2 ) */

            tau /= d__[2] + sqrt((d__1 = d__[2] * d__[2] + tau, abs(d__1)));
            *dsigma = d__[2] + tau;
            delta[1] = -(del + tau);
            delta[2] = -tau;
            work[1] = d__[1] + tau + d__[2];
            work[2] = d__[2] * 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 ) = DELTA( 1 ) / TEMP */
/*        DELTA( 2 ) = DELTA( 2 ) / TEMP */
    } else {

/*        Now I=2 */

        b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
        c__ = *rho * z__[2] * z__[2] * delsq;

/*        The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) */

        if (b > 0.) {
            tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
        } else {
            tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
        }

/*        The following TAU is DSIGMA - D( 2 ) */

        tau /= d__[2] + sqrt(d__[2] * d__[2] + tau);
        *dsigma = d__[2] + tau;
        delta[1] = -(del + tau);
        delta[2] = -tau;
        work[1] = d__[1] + tau + d__[2];
        work[2] = d__[2] * 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 ) = DELTA( 1 ) / TEMP */
/*        DELTA( 2 ) = DELTA( 2 ) / TEMP */
    }
    return 0;

/*     End of DLASD5 */

} /* dlasd5_ */