Update to 2.0.0 tree from current Fremantle build

[opencv] / 3rdparty / lapack / dlarrr.c

diff --git a/3rdparty/lapack/dlarrr.c b/3rdparty/lapack/dlarrr.c
new file mode 100644 (file)
index 0000000..5aee2a1
--- /dev/null
+++ b/3rdparty/lapack/dlarrr.c
@@ -0,0 +1,163 @@
+#include "clapack.h"
+
+/* Subroutine */ int dlarrr_(integer *n, doublereal *d__, doublereal *e, 
+       integer *info)
+{
+    /* System generated locals */
+    integer i__1;
+    doublereal d__1;
+
+    /* Builtin functions */
+    double sqrt(doublereal);
+
+    /* Local variables */
+    integer i__;
+    doublereal eps, tmp, tmp2, rmin;
+    extern doublereal dlamch_(char *);
+    doublereal offdig, safmin;
+    logical yesrel;
+    doublereal smlnum, offdig2;
+
+
+/*  -- LAPACK auxiliary routine (version 3.1) -- */
+/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/*     November 2006 */
+
+/*     .. Scalar Arguments .. */
+/*     .. */
+/*     .. Array Arguments .. */
+/*     .. */
+
+
+/*  Purpose */
+/*  ======= */
+
+/*  Perform tests to decide whether the symmetric tridiagonal matrix T */
+/*  warrants expensive computations which guarantee high relative accuracy */
+/*  in the eigenvalues. */
+
+/*  Arguments */
+/*  ========= */
+
+/*  N       (input) INTEGER */
+/*          The order of the matrix. N > 0. */
+
+/*  D       (input) DOUBLE PRECISION array, dimension (N) */
+/*          The N diagonal elements of the tridiagonal matrix T. */
+
+/*  E       (input/output) DOUBLE PRECISION array, dimension (N) */
+/*          On entry, the first (N-1) entries contain the subdiagonal */
+/*          elements of the tridiagonal matrix T; E(N) is set to ZERO. */
+
+/*  INFO    (output) INTEGER */
+/*          INFO = 0(default) : the matrix warrants computations preserving */
+/*                              relative accuracy. */
+/*          INFO = 1          : the matrix warrants computations guaranteeing */
+/*                              only absolute accuracy. */
+
+/*  Further Details */
+/*  =============== */
+
+/*  Based on contributions by */
+/*     Beresford Parlett, University of California, Berkeley, USA */
+/*     Jim Demmel, University of California, Berkeley, USA */
+/*     Inderjit Dhillon, University of Texas, Austin, USA */
+/*     Osni Marques, LBNL/NERSC, USA */
+/*     Christof Voemel, University of California, Berkeley, USA */
+
+/*  ===================================================================== */
+
+/*     .. Parameters .. */
+/*     .. */
+/*     .. Local Scalars .. */
+/*     .. */
+/*     .. External Functions .. */
+/*     .. */
+/*     .. Intrinsic Functions .. */
+/*     .. */
+/*     .. Executable Statements .. */
+
+/*     As a default, do NOT go for relative-accuracy preserving computations. */
+    /* Parameter adjustments */
+    --e;
+    --d__;
+
+    /* Function Body */
+    *info = 1;
+    safmin = dlamch_("Safe minimum");
+    eps = dlamch_("Precision");
+    smlnum = safmin / eps;
+    rmin = sqrt(smlnum);
+/*     Tests for relative accuracy */
+
+/*     Test for scaled diagonal dominance */
+/*     Scale the diagonal entries to one and check whether the sum of the */
+/*     off-diagonals is less than one */
+
+/*     The sdd relative error bounds have a 1/(1- 2*x) factor in them, */
+/*     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative */
+/*     accuracy is promised.  In the notation of the code fragment below, */
+/*     1/(1 - (OFFDIG + OFFDIG2)) is the condition number. */
+/*     We don't think it is worth going into "sdd mode" unless the relative */
+/*     condition number is reasonable, not 1/macheps. */
+/*     The threshold should be compatible with other thresholds used in the */
+/*     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds */
+/*     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000 */
+/*     instead of the current OFFDIG + OFFDIG2 < 1 */
+
+    yesrel = TRUE_;
+    offdig = 0.;
+    tmp = sqrt((abs(d__[1])));
+    if (tmp < rmin) {
+       yesrel = FALSE_;
+    }
+    if (! yesrel) {
+       goto L11;
+    }
+    i__1 = *n;
+    for (i__ = 2; i__ <= i__1; ++i__) {
+       tmp2 = sqrt((d__1 = d__[i__], abs(d__1)));
+       if (tmp2 < rmin) {
+           yesrel = FALSE_;
+       }
+       if (! yesrel) {
+           goto L11;
+       }
+       offdig2 = (d__1 = e[i__ - 1], abs(d__1)) / (tmp * tmp2);
+       if (offdig + offdig2 >= .999) {
+           yesrel = FALSE_;
+       }
+       if (! yesrel) {
+           goto L11;
+       }
+       tmp = tmp2;
+       offdig = offdig2;
+/* L10: */
+    }
+L11:
+    if (yesrel) {
+       *info = 0;
+       return 0;
+    } else {
+    }
+
+
+/*     *** MORE TO BE IMPLEMENTED *** */
+
+
+/*     Test if the lower bidiagonal matrix L from T = L D L^T */
+/*     (zero shift facto) is well conditioned */
+
+
+/*     Test if the upper bidiagonal matrix U from T = U D U^T */
+/*     (zero shift facto) is well conditioned. */
+/*     In this case, the matrix needs to be flipped and, at the end */
+/*     of the eigenvector computation, the flip needs to be applied */
+/*     to the computed eigenvectors (and the support) */
+
+
+    return 0;
+
+/*     END OF DLARRR */
+
+} /* dlarrr_ */