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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int slarrr_(integer *n, real *d__, real *e, integer *info)
   4 {
   5     /* System generated locals */
   6     integer i__1;
   7     real r__1;
   8
   9     /* Builtin functions */
  10     double sqrt(doublereal);
  11
  12     /* Local variables */
  13     integer i__;
  14     real eps, tmp, tmp2, rmin, offdig;
  15     extern doublereal slamch_(char *);
  16     real safmin;
  17     logical yesrel;
  18     real smlnum, offdig2;
  19
  20
  21 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  22 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  23 /*     November 2006 */
  24
  25 /*     .. Scalar Arguments .. */
  26 /*     .. */
  27 /*     .. Array Arguments .. */
  28 /*     .. */
  29
  30
  31 /*  Purpose */
  32 /*  ======= */
  33
  34 /*  Perform tests to decide whether the symmetric tridiagonal matrix T */
  35 /*  warrants expensive computations which guarantee high relative accuracy */
  36 /*  in the eigenvalues. */
  37
  38 /*  Arguments */
  39 /*  ========= */
  40
  41 /*  N       (input) INTEGER */
  42 /*          The order of the matrix. N > 0. */
  43
  44 /*  D       (input) REAL             array, dimension (N) */
  45 /*          The N diagonal elements of the tridiagonal matrix T. */
  46
  47 /*  E       (input/output) REAL             array, dimension (N) */
  48 /*          On entry, the first (N-1) entries contain the subdiagonal */
  49 /*          elements of the tridiagonal matrix T; E(N) is set to ZERO. */
  50
  51 /*  INFO    (output) INTEGER */
  52 /*          INFO = 0(default) : the matrix warrants computations preserving */
  53 /*                              relative accuracy. */
  54 /*          INFO = 1          : the matrix warrants computations guaranteeing */
  55 /*                              only absolute accuracy. */
  56
  57 /*  Further Details */
  58 /*  =============== */
  59
  60 /*  Based on contributions by */
  61 /*     Beresford Parlett, University of California, Berkeley, USA */
  62 /*     Jim Demmel, University of California, Berkeley, USA */
  63 /*     Inderjit Dhillon, University of Texas, Austin, USA */
  64 /*     Osni Marques, LBNL/NERSC, USA */
  65 /*     Christof Voemel, University of California, Berkeley, USA */
  66
  67 /*  ===================================================================== */
  68
  69 /*     .. Parameters .. */
  70 /*     .. */
  71 /*     .. Local Scalars .. */
  72 /*     .. */
  73 /*     .. External Functions .. */
  74 /*     .. */
  75 /*     .. Intrinsic Functions .. */
  76 /*     .. */
  77 /*     .. Executable Statements .. */
  78
  79 /*     As a default, do NOT go for relative-accuracy preserving computations. */
  80     /* Parameter adjustments */
  81     --e;
  82     --d__;
  83
  84     /* Function Body */
  85     *info = 1;
  86     safmin = slamch_("Safe minimum");
  87     eps = slamch_("Precision");
  88     smlnum = safmin / eps;
  89     rmin = sqrt(smlnum);
  90 /*     Tests for relative accuracy */
  91
  92 /*     Test for scaled diagonal dominance */
  93 /*     Scale the diagonal entries to one and check whether the sum of the */
  94 /*     off-diagonals is less than one */
  95
  96 /*     The sdd relative error bounds have a 1/(1- 2*x) factor in them, */
  97 /*     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative */
  98 /*     accuracy is promised.  In the notation of the code fragment below, */
  99 /*     1/(1 - (OFFDIG + OFFDIG2)) is the condition number. */
 100 /*     We don't think it is worth going into "sdd mode" unless the relative */
 101 /*     condition number is reasonable, not 1/macheps. */
 102 /*     The threshold should be compatible with other thresholds used in the */
 103 /*     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds */
 104 /*     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000 */
 105 /*     instead of the current OFFDIG + OFFDIG2 < 1 */
 106
 107     yesrel = TRUE_;
 108     offdig = 0.f;
 109     tmp = sqrt((dabs(d__[1])));
 110     if (tmp < rmin) {
 111         yesrel = FALSE_;
 112     }
 113     if (! yesrel) {
 114         goto L11;
 115     }
 116     i__1 = *n;
 117     for (i__ = 2; i__ <= i__1; ++i__) {
 118         tmp2 = sqrt((r__1 = d__[i__], dabs(r__1)));
 119         if (tmp2 < rmin) {
 120             yesrel = FALSE_;
 121         }
 122         if (! yesrel) {
 123             goto L11;
 124         }
 125         offdig2 = (r__1 = e[i__ - 1], dabs(r__1)) / (tmp * tmp2);
 126         if (offdig + offdig2 >= .999f) {
 127             yesrel = FALSE_;
 128         }
 129         if (! yesrel) {
 130             goto L11;
 131         }
 132         tmp = tmp2;
 133         offdig = offdig2;
 134 /* L10: */
 135     }
 136 L11:
 137     if (yesrel) {
 138         *info = 0;
 139         return 0;
 140     } else {
 141     }
 142
 143
 144 /*     *** MORE TO BE IMPLEMENTED *** */
 145
 146
 147 /*     Test if the lower bidiagonal matrix L from T = L D L^T */
 148 /*     (zero shift facto) is well conditioned */
 149
 150
 151 /*     Test if the upper bidiagonal matrix U from T = U D U^T */
 152 /*     (zero shift facto) is well conditioned. */
 153 /*     In this case, the matrix needs to be flipped and, at the end */
 154 /*     of the eigenvector computation, the flip needs to be applied */
 155 /*     to the computed eigenvectors (and the support) */
 156
 157
 158     return 0;
 159
 160 /*     END OF SLARRR */
 161
 162 } /* slarrr_ */