X-Git-Url: https://vcs.maemo.org/git/?a=blobdiff_plain;f=3rdparty%2Flapack%2Fdlaed5.c;fp=3rdparty%2Flapack%2Fdlaed5.c;h=a4be3eaa4257aea8864f243db4b60bd48c01c1fe;hb=e4c14cdbdf2fe805e79cd96ded236f57e7b89060;hp=0000000000000000000000000000000000000000;hpb=454138ff8a20f6edb9b65a910101403d8b520643;p=opencv

diff --git a/3rdparty/lapack/dlaed5.c b/3rdparty/lapack/dlaed5.c
new file mode 100644
index 0000000..a4be3ea
--- /dev/null
+++ b/3rdparty/lapack/dlaed5.c
@@ -0,0 +1,135 @@
+#include "clapack.h"
+
+/* Subroutine */ int dlaed5_(integer *i__, doublereal *d__, doublereal *z__, 
+	doublereal *delta, doublereal *rho, doublereal *dlam)
+{
+    /* System generated locals */
+    doublereal d__1;
+
+    /* Builtin functions */
+    double sqrt(doublereal);
+
+    /* Local variables */
+    doublereal b, c__, w, del, tau, temp;
+
+
+/*  -- LAPACK routine (version 3.1) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/*     November 2006 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+/*  Purpose */
+/*  ======= */
+
+/*  This subroutine computes the I-th eigenvalue of a symmetric rank-one */
+/*  modification of a 2-by-2 diagonal matrix */
+
+/*             diag( D )  +  RHO *  Z * transpose(Z) . */
+
+/*  The diagonal elements in the array D are assumed to satisfy */
+
+/*             D(i) < D(j)  for  i < j . */
+
+/*  We also assume RHO > 0 and that the Euclidean norm of the vector */
+/*  Z is one. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  I      (input) INTEGER */
+/*         The index of the eigenvalue to be computed.  I = 1 or I = 2. */
+
+/*  D      (input) DOUBLE PRECISION array, dimension (2) */
+/*         The original eigenvalues.  We assume D(1) < D(2). */
+
+/*  Z      (input) DOUBLE PRECISION array, dimension (2) */
+/*         The components of the updating vector. */
+
+/*  DELTA  (output) DOUBLE PRECISION array, dimension (2) */
+/*         The vector DELTA contains the information necessary */
+/*         to construct the eigenvectors. */
+
+/*  RHO    (input) DOUBLE PRECISION */
+/*         The scalar in the symmetric updating formula. */
+
+/*  DLAM   (output) DOUBLE PRECISION */
+/*         The computed lambda_I, the I-th updated eigenvalue. */
+
+/*  Further Details */
+/*  =============== */
+
+/*  Based on contributions by */
+/*     Ren-Cang Li, Computer Science Division, University of California */
+/*     at Berkeley, USA */
+
+/*  ===================================================================== */
+
+/*     .. Parameters .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+    /* Parameter adjustments */
+    --delta;
+    --z__;
+    --d__;
+
+    /* Function Body */
+    del = d__[2] - d__[1];
+    if (*i__ == 1) {
+	w = *rho * 2. * (z__[2] * z__[2] - z__[1] * z__[1]) / del + 1.;
+	if (w > 0.) {
+	    b = del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
+	    c__ = *rho * z__[1] * z__[1] * del;
+
+/*           B > ZERO, always */
+
+	    tau = c__ * 2. / (b + sqrt((d__1 = b * b - c__ * 4., abs(d__1))));
+	    *dlam = d__[1] + tau;
+	    delta[1] = -z__[1] / tau;
+	    delta[2] = z__[2] / (del - tau);
+	} else {
+	    b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
+	    c__ = *rho * z__[2] * z__[2] * del;
+	    if (b > 0.) {
+		tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
+	    } else {
+		tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
+	    }
+	    *dlam = d__[2] + tau;
+	    delta[1] = -z__[1] / (del + tau);
+	    delta[2] = -z__[2] / tau;
+	}
+	temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
+	delta[1] /= temp;
+	delta[2] /= temp;
+    } else {
+
+/*     Now I=2 */
+
+	b = -del + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
+	c__ = *rho * z__[2] * z__[2] * del;
+	if (b > 0.) {
+	    tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
+	} else {
+	    tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
+	}
+	*dlam = d__[2] + tau;
+	delta[1] = -z__[1] / (del + tau);
+	delta[2] = -z__[2] / tau;
+	temp = sqrt(delta[1] * delta[1] + delta[2] * delta[2]);
+	delta[1] /= temp;
+	delta[2] /= temp;
+    }
+    return 0;
+
+/*     End OF DLAED5 */
+
+} /* dlaed5_ */