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

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