--- /dev/null
+#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_ */