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

   1 #include "clapack.h"
   2
   3 integer slaneg_(integer *n, real *d__, real *lld, real *sigma, real *pivmin,
   4         integer *r__)
   5 {
   6     /* System generated locals */
   7     integer ret_val, i__1, i__2, i__3, i__4;
   8
   9     /* Local variables */
  10     integer j;
  11     real p, t;
  12     integer bj;
  13     real tmp;
  14     integer neg1, neg2;
  15     real bsav, gamma, dplus;
  16     integer negcnt;
  17     logical sawnan;
  18     extern logical sisnan_(real *);
  19     real dminus;
  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 /*  Purpose */
  32 /*  ======= */
  33
  34 /*  SLANEG computes the Sturm count, the number of negative pivots */
  35 /*  encountered while factoring tridiagonal T - sigma I = L D L^T. */
  36 /*  This implementation works directly on the factors without forming */
  37 /*  the tridiagonal matrix T.  The Sturm count is also the number of */
  38 /*  eigenvalues of T less than sigma. */
  39
  40 /*  This routine is called from SLARRB. */
  41
  42 /*  The current routine does not use the PIVMIN parameter but rather */
  43 /*  requires IEEE-754 propagation of Infinities and NaNs.  This */
  44 /*  routine also has no input range restrictions but does require */
  45 /*  default exception handling such that x/0 produces Inf when x is */
  46 /*  non-zero, and Inf/Inf produces NaN.  For more information, see: */
  47
  48 /*    Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in */
  49 /*    Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on */
  50 /*    Scientific Computing, v28, n5, 2006.  DOI 10.1137/050641624 */
  51 /*    (Tech report version in LAWN 172 with the same title.) */
  52
  53 /*  Arguments */
  54 /*  ========= */
  55
  56 /*  N       (input) INTEGER */
  57 /*          The order of the matrix. */
  58
  59 /*  D       (input) REAL             array, dimension (N) */
  60 /*          The N diagonal elements of the diagonal matrix D. */
  61
  62 /*  LLD     (input) REAL             array, dimension (N-1) */
  63 /*          The (N-1) elements L(i)*L(i)*D(i). */
  64
  65 /*  SIGMA   (input) REAL */
  66 /*          Shift amount in T - sigma I = L D L^T. */
  67
  68 /*  PIVMIN  (input) REAL */
  69 /*          The minimum pivot in the Sturm sequence.  May be used */
  70 /*          when zero pivots are encountered on non-IEEE-754 */
  71 /*          architectures. */
  72
  73 /*  R       (input) INTEGER */
  74 /*          The twist index for the twisted factorization that is used */
  75 /*          for the negcount. */
  76
  77 /*  Further Details */
  78 /*  =============== */
  79
  80 /*  Based on contributions by */
  81 /*     Osni Marques, LBNL/NERSC, USA */
  82 /*     Christof Voemel, University of California, Berkeley, USA */
  83 /*     Jason Riedy, University of California, Berkeley, USA */
  84
  85 /*  ===================================================================== */
  86
  87 /*     .. Parameters .. */
  88 /*     Some architectures propagate Infinities and NaNs very slowly, so */
  89 /*     the code computes counts in BLKLEN chunks.  Then a NaN can */
  90 /*     propagate at most BLKLEN columns before being detected.  This is */
  91 /*     not a general tuning parameter; it needs only to be just large */
  92 /*     enough that the overhead is tiny in common cases. */
  93 /*     .. */
  94 /*     .. Local Scalars .. */
  95 /*     .. */
  96 /*     .. Intrinsic Functions .. */
  97 /*     .. */
  98 /*     .. External Functions .. */
  99 /*     .. */
 100 /*     .. Executable Statements .. */
 101     /* Parameter adjustments */
 102     --lld;
 103     --d__;
 104
 105     /* Function Body */
 106     negcnt = 0;
 107 /*     I) upper part: L D L^T - SIGMA I = L+ D+ L+^T */
 108     t = -(*sigma);
 109     i__1 = *r__ - 1;
 110     for (bj = 1; bj <= i__1; bj += 128) {
 111         neg1 = 0;
 112         bsav = t;
 113 /* Computing MIN */
 114         i__3 = bj + 127, i__4 = *r__ - 1;
 115         i__2 = min(i__3,i__4);
 116         for (j = bj; j <= i__2; ++j) {
 117             dplus = d__[j] + t;
 118             if (dplus < 0.f) {
 119                 ++neg1;
 120             }
 121             tmp = t / dplus;
 122             t = tmp * lld[j] - *sigma;
 123 /* L21: */
 124         }
 125         sawnan = sisnan_(&t);
 126 /*     Run a slower version of the above loop if a NaN is detected. */
 127 /*     A NaN should occur only with a zero pivot after an infinite */
 128 /*     pivot.  In that case, substituting 1 for T/DPLUS is the */
 129 /*     correct limit. */
 130         if (sawnan) {
 131             neg1 = 0;
 132             t = bsav;
 133 /* Computing MIN */
 134             i__3 = bj + 127, i__4 = *r__ - 1;
 135             i__2 = min(i__3,i__4);
 136             for (j = bj; j <= i__2; ++j) {
 137                 dplus = d__[j] + t;
 138                 if (dplus < 0.f) {
 139                     ++neg1;
 140                 }
 141                 tmp = t / dplus;
 142                 if (sisnan_(&tmp)) {
 143                     tmp = 1.f;
 144                 }
 145                 t = tmp * lld[j] - *sigma;
 146 /* L22: */
 147             }
 148         }
 149         negcnt += neg1;
 150 /* L210: */
 151     }
 152
 153 /*     II) lower part: L D L^T - SIGMA I = U- D- U-^T */
 154     p = d__[*n] - *sigma;
 155     i__1 = *r__;
 156     for (bj = *n - 1; bj >= i__1; bj += -128) {
 157         neg2 = 0;
 158         bsav = p;
 159 /* Computing MAX */
 160         i__3 = bj - 127;
 161         i__2 = max(i__3,*r__);
 162         for (j = bj; j >= i__2; --j) {
 163             dminus = lld[j] + p;
 164             if (dminus < 0.f) {
 165                 ++neg2;
 166             }
 167             tmp = p / dminus;
 168             p = tmp * d__[j] - *sigma;
 169 /* L23: */
 170         }
 171         sawnan = sisnan_(&p);
 172 /*     As above, run a slower version that substitutes 1 for Inf/Inf. */
 173
 174         if (sawnan) {
 175             neg2 = 0;
 176             p = bsav;
 177 /* Computing MAX */
 178             i__3 = bj - 127;
 179             i__2 = max(i__3,*r__);
 180             for (j = bj; j >= i__2; --j) {
 181                 dminus = lld[j] + p;
 182                 if (dminus < 0.f) {
 183                     ++neg2;
 184                 }
 185                 tmp = p / dminus;
 186                 if (sisnan_(&tmp)) {
 187                     tmp = 1.f;
 188                 }
 189                 p = tmp * d__[j] - *sigma;
 190 /* L24: */
 191             }
 192         }
 193         negcnt += neg2;
 194 /* L230: */
 195     }
 196
 197 /*     III) Twist index */
 198 /*       T was shifted by SIGMA initially. */
 199     gamma = t + *sigma + p;
 200     if (gamma < 0.f) {
 201         ++negcnt;
 202     }
 203     ret_val = negcnt;
 204     return ret_val;
 205 } /* slaneg_ */