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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlarrk_(integer *n, integer *iw, doublereal *gl,
   4         doublereal *gu, doublereal *d__, doublereal *e2, doublereal *pivmin,
   5         doublereal *reltol, doublereal *w, doublereal *werr, integer *info)
   6 {
   7     /* System generated locals */
   8     integer i__1;
   9     doublereal d__1, d__2;
  10
  11     /* Builtin functions */
  12     double log(doublereal);
  13
  14     /* Local variables */
  15     integer i__, it;
  16     doublereal mid, eps, tmp1, tmp2, left, atoli, right;
  17     integer itmax;
  18     doublereal rtoli, tnorm;
  19     extern doublereal dlamch_(char *);
  20     integer negcnt;
  21
  22
  23 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  24 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  25 /*     November 2006 */
  26
  27 /*     .. Scalar Arguments .. */
  28 /*     .. */
  29 /*     .. Array Arguments .. */
  30 /*     .. */
  31
  32 /*  Purpose */
  33 /*  ======= */
  34
  35 /*  DLARRK computes one eigenvalue of a symmetric tridiagonal */
  36 /*  matrix T to suitable accuracy. This is an auxiliary code to be */
  37 /*  called from DSTEMR. */
  38
  39 /*  To avoid overflow, the matrix must be scaled so that its */
  40 /*  largest element is no greater than overflow**(1/2) * */
  41 /*  underflow**(1/4) in absolute value, and for greatest */
  42 /*  accuracy, it should not be much smaller than that. */
  43
  44 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
  45 /*  Matrix", Report CS41, Computer Science Dept., Stanford */
  46 /*  University, July 21, 1966. */
  47
  48 /*  Arguments */
  49 /*  ========= */
  50
  51 /*  N       (input) INTEGER */
  52 /*          The order of the tridiagonal matrix T.  N >= 0. */
  53
  54 /*  IW      (input) INTEGER */
  55 /*          The index of the eigenvalues to be returned. */
  56
  57 /*  GL      (input) DOUBLE PRECISION */
  58 /*  GU      (input) DOUBLE PRECISION */
  59 /*          An upper and a lower bound on the eigenvalue. */
  60
  61 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
  62 /*          The n diagonal elements of the tridiagonal matrix T. */
  63
  64 /*  E2      (input) DOUBLE PRECISION array, dimension (N-1) */
  65 /*          The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
  66
  67 /*  PIVMIN  (input) DOUBLE PRECISION */
  68 /*          The minimum pivot allowed in the Sturm sequence for T. */
  69
  70 /*  RELTOL  (input) DOUBLE PRECISION */
  71 /*          The minimum relative width of an interval.  When an interval */
  72 /*          is narrower than RELTOL times the larger (in */
  73 /*          magnitude) endpoint, then it is considered to be */
  74 /*          sufficiently small, i.e., converged.  Note: this should */
  75 /*          always be at least radix*machine epsilon. */
  76
  77 /*  W       (output) DOUBLE PRECISION */
  78
  79 /*  WERR    (output) DOUBLE PRECISION */
  80 /*          The error bound on the corresponding eigenvalue approximation */
  81 /*          in W. */
  82
  83 /*  INFO    (output) INTEGER */
  84 /*          = 0:       Eigenvalue converged */
  85 /*          = -1:      Eigenvalue did NOT converge */
  86
  87 /*  Internal Parameters */
  88 /*  =================== */
  89
  90 /*  FUDGE   DOUBLE PRECISION, default = 2 */
  91 /*          A "fudge factor" to widen the Gershgorin intervals. */
  92
  93 /*  ===================================================================== */
  94
  95 /*     .. Parameters .. */
  96 /*     .. */
  97 /*     .. Local Scalars .. */
  98 /*     .. */
  99 /*     .. External Functions .. */
 100 /*     .. */
 101 /*     .. Intrinsic Functions .. */
 102 /*     .. */
 103 /*     .. Executable Statements .. */
 104
 105 /*     Get machine constants */
 106     /* Parameter adjustments */
 107     --e2;
 108     --d__;
 109
 110     /* Function Body */
 111     eps = dlamch_("P");
 112 /* Computing MAX */
 113     d__1 = abs(*gl), d__2 = abs(*gu);
 114     tnorm = max(d__1,d__2);
 115     rtoli = *reltol;
 116     atoli = *pivmin * 4.;
 117     itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.)) + 2;
 118     *info = -1;
 119     left = *gl - tnorm * 2. * eps * *n - *pivmin * 4.;
 120     right = *gu + tnorm * 2. * eps * *n + *pivmin * 4.;
 121     it = 0;
 122 L10:
 123
 124 /*     Check if interval converged or maximum number of iterations reached */
 125
 126     tmp1 = (d__1 = right - left, abs(d__1));
 127 /* Computing MAX */
 128     d__1 = abs(right), d__2 = abs(left);
 129     tmp2 = max(d__1,d__2);
 130 /* Computing MAX */
 131     d__1 = max(atoli,*pivmin), d__2 = rtoli * tmp2;
 132     if (tmp1 < max(d__1,d__2)) {
 133         *info = 0;
 134         goto L30;
 135     }
 136     if (it > itmax) {
 137         goto L30;
 138     }
 139
 140 /*     Count number of negative pivots for mid-point */
 141
 142     ++it;
 143     mid = (left + right) * .5;
 144     negcnt = 0;
 145     tmp1 = d__[1] - mid;
 146     if (abs(tmp1) < *pivmin) {
 147         tmp1 = -(*pivmin);
 148     }
 149     if (tmp1 <= 0.) {
 150         ++negcnt;
 151     }
 152
 153     i__1 = *n;
 154     for (i__ = 2; i__ <= i__1; ++i__) {
 155         tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid;
 156         if (abs(tmp1) < *pivmin) {
 157             tmp1 = -(*pivmin);
 158         }
 159         if (tmp1 <= 0.) {
 160             ++negcnt;
 161         }
 162 /* L20: */
 163     }
 164     if (negcnt >= *iw) {
 165         right = mid;
 166     } else {
 167         left = mid;
 168     }
 169     goto L10;
 170 L30:
 171
 172 /*     Converged or maximum number of iterations reached */
 173
 174     *w = (left + right) * .5;
 175     *werr = (d__1 = right - left, abs(d__1)) * .5;
 176     return 0;
 177
 178 /*     End of DLARRK */
 179
 180 } /* dlarrk_ */