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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int slasd5_(integer *i__, real *d__, real *z__, real *delta,
   4         real *rho, real *dsigma, real *work)
   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, delsq;
  14
  15
  16 /*  -- LAPACK auxiliary 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 square root of the I-th eigenvalue */
  29 /*  of a positive symmetric rank-one modification of a 2-by-2 diagonal */
  30 /*  matrix */
  31
  32 /*             diag( D ) * diag( D ) +  RHO *  Z * transpose(Z) . */
  33
  34 /*  The diagonal entries in the array D are assumed to satisfy */
  35
  36 /*             0 <= D(i) < D(j)  for  i < j . */
  37
  38 /*  We also assume RHO > 0 and that the Euclidean norm of the vector */
  39 /*  Z is one. */
  40
  41 /*  Arguments */
  42 /*  ========= */
  43
  44 /*  I      (input) INTEGER */
  45 /*         The index of the eigenvalue to be computed.  I = 1 or I = 2. */
  46
  47 /*  D      (input) REAL array, dimension (2) */
  48 /*         The original eigenvalues.  We assume 0 <= D(1) < D(2). */
  49
  50 /*  Z      (input) REAL array, dimension (2) */
  51 /*         The components of the updating vector. */
  52
  53 /*  DELTA  (output) REAL array, dimension (2) */
  54 /*         Contains (D(j) - sigma_I) in its  j-th component. */
  55 /*         The vector DELTA contains the information necessary */
  56 /*         to construct the eigenvectors. */
  57
  58 /*  RHO    (input) REAL */
  59 /*         The scalar in the symmetric updating formula. */
  60
  61 /*  DSIGMA (output) REAL */
  62 /*         The computed sigma_I, the I-th updated eigenvalue. */
  63
  64 /*  WORK   (workspace) REAL array, dimension (2) */
  65 /*         WORK contains (D(j) + sigma_I) in its  j-th component. */
  66
  67 /*  Further Details */
  68 /*  =============== */
  69
  70 /*  Based on contributions by */
  71 /*     Ren-Cang Li, Computer Science Division, University of California */
  72 /*     at Berkeley, USA */
  73
  74 /*  ===================================================================== */
  75
  76 /*     .. Parameters .. */
  77 /*     .. */
  78 /*     .. Local Scalars .. */
  79 /*     .. */
  80 /*     .. Intrinsic Functions .. */
  81 /*     .. */
  82 /*     .. Executable Statements .. */
  83
  84     /* Parameter adjustments */
  85     --work;
  86     --delta;
  87     --z__;
  88     --d__;
  89
  90     /* Function Body */
  91     del = d__[2] - d__[1];
  92     delsq = del * (d__[2] + d__[1]);
  93     if (*i__ == 1) {
  94         w = *rho * 4.f * (z__[2] * z__[2] / (d__[1] + d__[2] * 3.f) - z__[1] *
  95                  z__[1] / (d__[1] * 3.f + d__[2])) / del + 1.f;
  96         if (w > 0.f) {
  97             b = delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  98             c__ = *rho * z__[1] * z__[1] * delsq;
  99
 100 /*           B > ZERO, always */
 101
 102 /*           The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 ) */
 103
 104             tau = c__ * 2.f / (b + sqrt((r__1 = b * b - c__ * 4.f, dabs(r__1))
 105                     ));
 106
 107 /*           The following TAU is DSIGMA - D( 1 ) */
 108
 109             tau /= d__[1] + sqrt(d__[1] * d__[1] + tau);
 110             *dsigma = d__[1] + tau;
 111             delta[1] = -tau;
 112             delta[2] = del - tau;
 113             work[1] = d__[1] * 2.f + tau;
 114             work[2] = d__[1] + tau + d__[2];
 115 /*           DELTA( 1 ) = -Z( 1 ) / TAU */
 116 /*           DELTA( 2 ) = Z( 2 ) / ( DEL-TAU ) */
 117         } else {
 118             b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
 119             c__ = *rho * z__[2] * z__[2] * delsq;
 120
 121 /*           The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) */
 122
 123             if (b > 0.f) {
 124                 tau = c__ * -2.f / (b + sqrt(b * b + c__ * 4.f));
 125             } else {
 126                 tau = (b - sqrt(b * b + c__ * 4.f)) / 2.f;
 127             }
 128
 129 /*           The following TAU is DSIGMA - D( 2 ) */
 130
 131             tau /= d__[2] + sqrt((r__1 = d__[2] * d__[2] + tau, dabs(r__1)));
 132             *dsigma = d__[2] + tau;
 133             delta[1] = -(del + tau);
 134             delta[2] = -tau;
 135             work[1] = d__[1] + tau + d__[2];
 136             work[2] = d__[2] * 2.f + tau;
 137 /*           DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) */
 138 /*           DELTA( 2 ) = -Z( 2 ) / TAU */
 139         }
 140 /*        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) */
 141 /*        DELTA( 1 ) = DELTA( 1 ) / TEMP */
 142 /*        DELTA( 2 ) = DELTA( 2 ) / TEMP */
 143     } else {
 144
 145 /*        Now I=2 */
 146
 147         b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
 148         c__ = *rho * z__[2] * z__[2] * delsq;
 149
 150 /*        The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 ) */
 151
 152         if (b > 0.f) {
 153             tau = (b + sqrt(b * b + c__ * 4.f)) / 2.f;
 154         } else {
 155             tau = c__ * 2.f / (-b + sqrt(b * b + c__ * 4.f));
 156         }
 157
 158 /*        The following TAU is DSIGMA - D( 2 ) */
 159
 160         tau /= d__[2] + sqrt(d__[2] * d__[2] + tau);
 161         *dsigma = d__[2] + tau;
 162         delta[1] = -(del + tau);
 163         delta[2] = -tau;
 164         work[1] = d__[1] + tau + d__[2];
 165         work[2] = d__[2] * 2.f + tau;
 166 /*        DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU ) */
 167 /*        DELTA( 2 ) = -Z( 2 ) / TAU */
 168 /*        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) ) */
 169 /*        DELTA( 1 ) = DELTA( 1 ) / TEMP */
 170 /*        DELTA( 2 ) = DELTA( 2 ) / TEMP */
 171     }
 172     return 0;
 173
 174 /*     End of SLASD5 */
 175
 176 } /* slasd5_ */