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

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