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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int slaed5_(integer *i__, real *d__, real *z__, real *delta,
   4         real *rho, real *dlam)
   5 {
   6     /* System generated locals */
   7     real r__1;
   8
   9     /* Builtin functions */
  10     double sqrt(doublereal);
  11
  12     /* Local variables */
  13     real 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) REAL array, dimension (2) */
  47 /*         The original eigenvalues.  We assume D(1) < D(2). */
  48
  49 /*  Z      (input) REAL array, dimension (2) */
  50 /*         The components of the updating vector. */
  51
  52 /*  DELTA  (output) REAL array, dimension (2) */
  53 /*         The vector DELTA contains the information necessary */
  54 /*         to construct the eigenvectors. */
  55
  56 /*  RHO    (input) REAL */
  57 /*         The scalar in the symmetric updating formula. */
  58
  59 /*  DLAM   (output) REAL */
  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.f * (z__[2] * z__[2] - z__[1] * z__[1]) / del + 1.f;
  88         if (w > 0.f) {
  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.f / (b + sqrt((r__1 = b * b - c__ * 4.f, dabs(r__1))
  95                     ));
  96             *dlam = d__[1] + tau;
  97             delta[1] = -z__[1] / tau;
  98             delta[2] = z__[2] / (del - tau);
  99         } else {
 100             b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
 101             c__ = *rho * z__[2] * z__[2] * del;
 102             if (b > 0.f) {
 103                 tau = c__ * -2.f / (b + sqrt(b * b + c__ * 4.f));
 104             } else {
 105                 tau = (b - sqrt(b * b + c__ * 4.f)) / 2.f;
 106             }
 107             *dlam = d__[2] + tau;
 108             delta[1] = -z__[1] / (del + tau);
 109             delta[2] = -z__[2] / tau;
 110         }
 111         temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
 112         delta[1] /= temp;
 113         delta[2] /= temp;
 114     } else {
 115
 116 /*     Now I=2 */
 117
 118         b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
 119         c__ = *rho * z__[2] * z__[2] * del;
 120         if (b > 0.f) {
 121             tau = (b + sqrt(b * b + c__ * 4.f)) / 2.f;
 122         } else {
 123             tau = c__ * 2.f / (-b + sqrt(b * b + c__ * 4.f));
 124         }
 125         *dlam = d__[2] + tau;
 126         delta[1] = -z__[1] / (del + tau);
 127         delta[2] = -z__[2] / tau;
 128         temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
 129         delta[1] /= temp;
 130         delta[2] /= temp;
 131     }
 132     return 0;
 133
 134 /*     End OF SLAED5 */
 135
 136 } /* slaed5_ */