3rdparty/lapack/slaneg.c

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