/usr/include/rheolef/tqli.h is in librheolef-dev 6.7-6.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 | ///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
///
/// =========================================================================
#include "rheolef/compiler.h"
#include <iterator>
#include <iostream>
#include "rheolef/pythag.h"
template <class T>
inline T SIGN(T a, T b) { return b >= T(0) ? fabs(a) : -fabs(a); }
template <class Iterator1, class Iterator2, class Size>
void tqli (Iterator1 d, Iterator2 e, Size n)
{
typedef typename std::iterator_traits<Iterator1>::value_type T;
static const Size max_iter = 30;
static const T eps = std::numeric_limits<T>::epsilon();
Size m,l,iter,i;
T s,r,p,g,f,dd,c,b;
for (i=2;i<=n;i++) e[i-1]=e[i];
e[n]=0.0;
#ifdef TO_CLEAN
// dump:
std::cout << std::setprecision(31);
for (i=1; i <= n; i++) {
std::cout << "d("<<i<<") = " << d[i] << std::endl;
}
for (i=1; i <= n; i++) {
std::cout << "e("<<i<<") = " << e[i] << std::endl;
}
#endif // TO_CLEAN
for (l=1; l<=n; l++) {
iter=0;
do {
for (m = l; m <= n-1; m++) {
dd = fabs(d[m]) + fabs(d[m+1]);
if (fabs(e[m]) < dd && fabs(e[m])*dd < eps*dd) break; // i.e. e[m] << dd
}
if (m != l) {
iter++;
check_macro (iter < max_iter, "too many ("<<max_iter<<") iterations in tqli");
g = (d[l+1]-d[l])/(2.0*e[l]);
r = pythag (g, T(1));
g = d[m] - d[l] + e[l]/(g+SIGN(r,g));
s = c = 1.0;
p = 0.0;
for (i = m-1; i>=l; i--) {
f = s*e[i];
b = c*e[i];
r=pythag(f,g);
e[i+1] = r;
if (r == T(0)) {
d[i+1] -= p;
e[m] = 0.0;
break;
}
s = f/r;
c = g/r;
g = d[i+1] - p;
r = (d[i] - g)*s + 2.0*c*b;
p=s*r;
d[i+1] = g + p;
g = c*r - b;
}
if (r == T(0) && i >= l) continue;
d[l] -= p;
e[l] = g;
e[m] = 0.0;
}
} while (m != l);
}
}
|