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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlarrc_(char *jobt, integer *n, doublereal *vl,
   4         doublereal *vu, doublereal *d__, doublereal *e, doublereal *pivmin,
   5         integer *eigcnt, integer *lcnt, integer *rcnt, integer *info)
   6 {
   7     /* System generated locals */
   8     integer i__1;
   9     doublereal d__1;
  10
  11     /* Local variables */
  12     integer i__;
  13     doublereal sl, su, tmp, tmp2;
  14     logical matt;
  15     extern logical lsame_(char *, char *);
  16     doublereal lpivot, rpivot;
  17
  18
  19 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  20 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  21 /*     November 2006 */
  22
  23 /*     .. Scalar Arguments .. */
  24 /*     .. */
  25 /*     .. Array Arguments .. */
  26 /*     .. */
  27
  28 /*  Purpose */
  29 /*  ======= */
  30
  31 /*  Find the number of eigenvalues of the symmetric tridiagonal matrix T */
  32 /*  that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T */
  33 /*  if JOBT = 'L'. */
  34
  35 /*  Arguments */
  36 /*  ========= */
  37
  38 /*  JOBT    (input) CHARACTER*1 */
  39 /*          = 'T':  Compute Sturm count for matrix T. */
  40 /*          = 'L':  Compute Sturm count for matrix L D L^T. */
  41
  42 /*  N       (input) INTEGER */
  43 /*          The order of the matrix. N > 0. */
  44
  45 /*  VL      (input) DOUBLE PRECISION */
  46 /*  VU      (input) DOUBLE PRECISION */
  47 /*          The lower and upper bounds for the eigenvalues. */
  48
  49 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
  50 /*          JOBT = 'T': The N diagonal elements of the tridiagonal matrix T. */
  51 /*          JOBT = 'L': The N diagonal elements of the diagonal matrix D. */
  52
  53 /*  E       (input) DOUBLE PRECISION array, dimension (N) */
  54 /*          JOBT = 'T': The N-1 offdiagonal elements of the matrix T. */
  55 /*          JOBT = 'L': The N-1 offdiagonal elements of the matrix L. */
  56
  57 /*  PIVMIN  (input) DOUBLE PRECISION */
  58 /*          The minimum pivot in the Sturm sequence for T. */
  59
  60 /*  EIGCNT  (output) INTEGER */
  61 /*          The number of eigenvalues of the symmetric tridiagonal matrix T */
  62 /*          that are in the interval (VL,VU] */
  63
  64 /*  LCNT    (output) INTEGER */
  65 /*  RCNT    (output) INTEGER */
  66 /*          The left and right negcounts of the interval. */
  67
  68 /*  INFO    (output) INTEGER */
  69
  70 /*  Further Details */
  71 /*  =============== */
  72
  73 /*  Based on contributions by */
  74 /*     Beresford Parlett, University of California, Berkeley, USA */
  75 /*     Jim Demmel, University of California, Berkeley, USA */
  76 /*     Inderjit Dhillon, University of Texas, Austin, USA */
  77 /*     Osni Marques, LBNL/NERSC, USA */
  78 /*     Christof Voemel, University of California, Berkeley, USA */
  79
  80 /*  ===================================================================== */
  81
  82 /*     .. Parameters .. */
  83 /*     .. */
  84 /*     .. Local Scalars .. */
  85 /*     .. */
  86 /*     .. External Functions .. */
  87 /*     .. */
  88 /*     .. Executable Statements .. */
  89
  90     /* Parameter adjustments */
  91     --e;
  92     --d__;
  93
  94     /* Function Body */
  95     *info = 0;
  96     *lcnt = 0;
  97     *rcnt = 0;
  98     *eigcnt = 0;
  99     matt = lsame_(jobt, "T");
 100     if (matt) {
 101 /*        Sturm sequence count on T */
 102         lpivot = d__[1] - *vl;
 103         rpivot = d__[1] - *vu;
 104         if (lpivot <= 0.) {
 105             ++(*lcnt);
 106         }
 107         if (rpivot <= 0.) {
 108             ++(*rcnt);
 109         }
 110         i__1 = *n - 1;
 111         for (i__ = 1; i__ <= i__1; ++i__) {
 112 /* Computing 2nd power */
 113             d__1 = e[i__];
 114             tmp = d__1 * d__1;
 115             lpivot = d__[i__ + 1] - *vl - tmp / lpivot;
 116             rpivot = d__[i__ + 1] - *vu - tmp / rpivot;
 117             if (lpivot <= 0.) {
 118                 ++(*lcnt);
 119             }
 120             if (rpivot <= 0.) {
 121                 ++(*rcnt);
 122             }
 123 /* L10: */
 124         }
 125     } else {
 126 /*        Sturm sequence count on L D L^T */
 127         sl = -(*vl);
 128         su = -(*vu);
 129         i__1 = *n - 1;
 130         for (i__ = 1; i__ <= i__1; ++i__) {
 131             lpivot = d__[i__] + sl;
 132             rpivot = d__[i__] + su;
 133             if (lpivot <= 0.) {
 134                 ++(*lcnt);
 135             }
 136             if (rpivot <= 0.) {
 137                 ++(*rcnt);
 138             }
 139             tmp = e[i__] * d__[i__] * e[i__];
 140
 141             tmp2 = tmp / lpivot;
 142             if (tmp2 == 0.) {
 143                 sl = tmp - *vl;
 144             } else {
 145                 sl = sl * tmp2 - *vl;
 146             }
 147
 148             tmp2 = tmp / rpivot;
 149             if (tmp2 == 0.) {
 150                 su = tmp - *vu;
 151             } else {
 152                 su = su * tmp2 - *vu;
 153             }
 154 /* L20: */
 155         }
 156         lpivot = d__[*n] + sl;
 157         rpivot = d__[*n] + su;
 158         if (lpivot <= 0.) {
 159             ++(*lcnt);
 160         }
 161         if (rpivot <= 0.) {
 162             ++(*rcnt);
 163         }
 164     }
 165     *eigcnt = *rcnt - *lcnt;
 166     return 0;
 167
 168 /*     end of DLARRC */
 169
 170 } /* dlarrc_ */