git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlaed5.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlaed5_(integer *i__, doublereal *d__, doublereal *z__,
   4         doublereal *delta, doublereal *rho, doublereal *dlam)
   5 {
   6     /* System generated locals */
   7     doublereal d__1;
   8
   9     /* Builtin functions */
  10     double sqrt(doublereal);
  11
  12     /* Local variables */
  13     doublereal b, c__, w, del, tau, temp;
  14
  15
  16 /*  -- LAPACK routine (version 3.1) -- */
  17 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  18 /*     November 2006 */
  19
  20 /*     .. Scalar Arguments .. */
  21 /*     .. */
  22 /*     .. Array Arguments .. */
  23 /*     .. */
  24
  25 /*  Purpose */
  26 /*  ======= */
  27
  28 /*  This subroutine computes the I-th eigenvalue of a symmetric rank-one */
  29 /*  modification of a 2-by-2 diagonal matrix */
  30
  31 /*             diag( D )  +  RHO *  Z * transpose(Z) . */
  32
  33 /*  The diagonal elements in the array D are assumed to satisfy */
  34
  35 /*             D(i) < D(j)  for  i < j . */
  36
  37 /*  We also assume RHO > 0 and that the Euclidean norm of the vector */
  38 /*  Z is one. */
  39
  40 /*  Arguments */
  41 /*  ========= */
  42
  43 /*  I      (input) INTEGER */
  44 /*         The index of the eigenvalue to be computed.  I = 1 or I = 2. */
  45
  46 /*  D      (input) DOUBLE PRECISION array, dimension (2) */
  47 /*         The original eigenvalues.  We assume D(1) < D(2). */
  48
  49 /*  Z      (input) DOUBLE PRECISION array, dimension (2) */
  50 /*         The components of the updating vector. */
  51
  52 /*  DELTA  (output) DOUBLE PRECISION array, dimension (2) */
  53 /*         The vector DELTA contains the information necessary */
  54 /*         to construct the eigenvectors. */
  55
  56 /*  RHO    (input) DOUBLE PRECISION */
  57 /*         The scalar in the symmetric updating formula. */
  58
  59 /*  DLAM   (output) DOUBLE PRECISION */
  60 /*         The computed lambda_I, the I-th updated eigenvalue. */
  61
  62 /*  Further Details */
  63 /*  =============== */
  64
  65 /*  Based on contributions by */
  66 /*     Ren-Cang Li, Computer Science Division, University of California */
  67 /*     at Berkeley, USA */
  68
  69 /*  ===================================================================== */
  70
  71 /*     .. Parameters .. */
  72 /*     .. */
  73 /*     .. Local Scalars .. */
  74 /*     .. */
  75 /*     .. Intrinsic Functions .. */
  76 /*     .. */
  77 /*     .. Executable Statements .. */
  78
  79     /* Parameter adjustments */
  80     --delta;
  81     --z__;
  82     --d__;
  83
  84     /* Function Body */
  85     del = d__[2] - d__[1];
  86     if (*i__ == 1) {
  87         w = *rho * 2. * (z__[2] * z__[2] - z__[1] * z__[1]) / del + 1.;
  88         if (w > 0.) {
  89             b = del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  90             c__ = *rho * z__[1] * z__[1] * del;
  91
  92 /*           B > ZERO, always */
  93
  94             tau = c__ * 2. / (b + sqrt((d__1 = b * b - c__ * 4., abs(d__1))));
  95             *dlam = d__[1] + tau;
  96             delta[1] = -z__[1] / tau;
  97             delta[2] = z__[2] / (del - tau);
  98         } else {
  99             b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
 100             c__ = *rho * z__[2] * z__[2] * del;
 101             if (b > 0.) {
 102                 tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
 103             } else {
 104                 tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
 105             }
 106             *dlam = d__[2] + tau;
 107             delta[1] = -z__[1] / (del + tau);
 108             delta[2] = -z__[2] / tau;
 109         }
 110         temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
 111         delta[1] /= temp;
 112         delta[2] /= temp;
 113     } else {
 114
 115 /*     Now I=2 */
 116
 117         b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
 118         c__ = *rho * z__[2] * z__[2] * del;
 119         if (b > 0.) {
 120             tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
 121         } else {
 122             tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
 123         }
 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]);
 128         delta[1] /= temp;
 129         delta[2] /= temp;
 130     }
 131     return 0;
 132
 133 /*     End OF DLAED5 */
 134
 135 } /* dlaed5_ */