slaed5.c
Go to the documentation of this file.
00001 /* slaed5.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int slaed5_(integer *i__, real *d__, real *z__, real *delta, 
00017         real *rho, real *dlam)
00018 {
00019     /* System generated locals */
00020     real r__1;
00021 
00022     /* Builtin functions */
00023     double sqrt(doublereal);
00024 
00025     /* Local variables */
00026     real b, c__, w, del, tau, temp;
00027 
00028 
00029 /*  -- LAPACK routine (version 3.2) -- */
00030 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00031 /*     November 2006 */
00032 
00033 /*     .. Scalar Arguments .. */
00034 /*     .. */
00035 /*     .. Array Arguments .. */
00036 /*     .. */
00037 
00038 /*  Purpose */
00039 /*  ======= */
00040 
00041 /*  This subroutine computes the I-th eigenvalue of a symmetric rank-one */
00042 /*  modification of a 2-by-2 diagonal matrix */
00043 
00044 /*             diag( D )  +  RHO *  Z * transpose(Z) . */
00045 
00046 /*  The diagonal elements in the array D are assumed to satisfy */
00047 
00048 /*             D(i) < D(j)  for  i < j . */
00049 
00050 /*  We also assume RHO > 0 and that the Euclidean norm of the vector */
00051 /*  Z is one. */
00052 
00053 /*  Arguments */
00054 /*  ========= */
00055 
00056 /*  I      (input) INTEGER */
00057 /*         The index of the eigenvalue to be computed.  I = 1 or I = 2. */
00058 
00059 /*  D      (input) REAL array, dimension (2) */
00060 /*         The original eigenvalues.  We assume D(1) < D(2). */
00061 
00062 /*  Z      (input) REAL array, dimension (2) */
00063 /*         The components of the updating vector. */
00064 
00065 /*  DELTA  (output) REAL array, dimension (2) */
00066 /*         The vector DELTA contains the information necessary */
00067 /*         to construct the eigenvectors. */
00068 
00069 /*  RHO    (input) REAL */
00070 /*         The scalar in the symmetric updating formula. */
00071 
00072 /*  DLAM   (output) REAL */
00073 /*         The computed lambda_I, the I-th updated eigenvalue. */
00074 
00075 /*  Further Details */
00076 /*  =============== */
00077 
00078 /*  Based on contributions by */
00079 /*     Ren-Cang Li, Computer Science Division, University of California */
00080 /*     at Berkeley, USA */
00081 
00082 /*  ===================================================================== */
00083 
00084 /*     .. Parameters .. */
00085 /*     .. */
00086 /*     .. Local Scalars .. */
00087 /*     .. */
00088 /*     .. Intrinsic Functions .. */
00089 /*     .. */
00090 /*     .. Executable Statements .. */
00091 
00092     /* Parameter adjustments */
00093     --delta;
00094     --z__;
00095     --d__;
00096 
00097     /* Function Body */
00098     del = d__[2] - d__[1];
00099     if (*i__ == 1) {
00100         w = *rho * 2.f * (z__[2] * z__[2] - z__[1] * z__[1]) / del + 1.f;
00101         if (w > 0.f) {
00102             b = del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
00103             c__ = *rho * z__[1] * z__[1] * del;
00104 
00105 /*           B > ZERO, always */
00106 
00107             tau = c__ * 2.f / (b + sqrt((r__1 = b * b - c__ * 4.f, dabs(r__1))
00108                     ));
00109             *dlam = d__[1] + tau;
00110             delta[1] = -z__[1] / tau;
00111             delta[2] = z__[2] / (del - tau);
00112         } else {
00113             b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
00114             c__ = *rho * z__[2] * z__[2] * del;
00115             if (b > 0.f) {
00116                 tau = c__ * -2.f / (b + sqrt(b * b + c__ * 4.f));
00117             } else {
00118                 tau = (b - sqrt(b * b + c__ * 4.f)) / 2.f;
00119             }
00120             *dlam = d__[2] + tau;
00121             delta[1] = -z__[1] / (del + tau);
00122             delta[2] = -z__[2] / tau;
00123         }
00124         temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
00125         delta[1] /= temp;
00126         delta[2] /= temp;
00127     } else {
00128 
00129 /*     Now I=2 */
00130 
00131         b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
00132         c__ = *rho * z__[2] * z__[2] * del;
00133         if (b > 0.f) {
00134             tau = (b + sqrt(b * b + c__ * 4.f)) / 2.f;
00135         } else {
00136             tau = c__ * 2.f / (-b + sqrt(b * b + c__ * 4.f));
00137         }
00138         *dlam = d__[2] + tau;
00139         delta[1] = -z__[1] / (del + tau);
00140         delta[2] = -z__[2] / tau;
00141         temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
00142         delta[1] /= temp;
00143         delta[2] /= temp;
00144     }
00145     return 0;
00146 
00147 /*     End OF SLAED5 */
00148 
00149 } /* slaed5_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:09