kennedy@ARISIA.XEROX.COM (Robert Kennedy) (08/18/89)
g++ version 1.36.0-, sun 4/110, SunOS 4.0
The following file causes the abort() at line 1201 of cplus-decl.c to
be called from function finish_file() in gcc-cc1plus, resulting in
compiler death from signal 6.
--------------- begin lin_solve.c ----------------
/*
routine to solve the d-dimensional linear system cx = b using Gaussian
elimination with scaled column pivoting. The algorithm used here
follows that given in Burden and Faires, pp. 326-328.
*/
/*
Should only include qflt.h and remes.h when variable-sized arrays work
*/
#include "poly.h"
#include "remes.h"
void lin_solve(qflt c[][max_deg + 2], qflt x[], qflt b[], int d)
{
/* Static to avoid construction every time */
static qflt a[max_deg + 2][max_deg + 3], /* augmented matrix */
s[max_deg + 2], /* scale factors */
absa, mx;
int i, j, k, p,
rp[d]; /* row pointers */
/*
Copy the entire mess into an augmented matrix for two reasons:
1) We don't want to sully the originals.
2) We want nice uniformity in the algorithm.
At the same time, we compute the maximal row entries.
*/
for (i = 0; i < d; i++)
{
s[i] = 0.0;
for (j = 0; j < d; j++)
{
a[i][j] = c[i][j];
if ((absa = fabs(a[i][j])) > s[i]) s[i] = absa;
}
if (s[i] == qflt_zero) printf("lin_solve: singular!\n");
a[i][d] = b[i];
rp[i] = i;
}
for (i = 0; i < d - 1; i++)
{
/*
Find pivot row.
*/
mx = qflt_zero;
for (j = i; j < d; j++)
if (mx < (absa = fabs(a[rp[j]][i] / s[rp[i]])))
{
p = j;
mx = absa;
}
if ((mx == qflt_zero) || (a[rp[p]][i] == qflt_zero))
printf("lin_solve: singular!\n");
if (rp[i] != rp[p])
{
k = rp[i];
rp[i] = rp[p];
rp[p] = k;
}
for (j = i + 1; j < d; j++)
{
absa = a[rp[j]][i] / a[rp[i]][i];
for (k = 0; k <= d; k++)
a[rp[j]][k] -= absa * a[rp[i]][k];
}
}
if (a[rp[d-1]][d-1] == qflt_zero) printf("lin_solve: singular!\n");
x[d-1] = a[rp[d-1]][d] / a[rp[d-1]][d-1];
for (i = d - 2; i >= 0; i--)
{
absa = qflt_zero;
for (j = i + 1; j < d; j++)
absa += a[rp[i]][j] * x[j];
x[i] = (a[rp[i]][d] - absa) / a[rp[i]][i];
}
}
-------------- end lin_solve.c -------------
-------------- begin poly.h -------------
/*
poly.h
class declaration for polynomial functions
*/
#include "../qflt.h"
const int max_deg = 3;
class p_func;
class poly
{
int d;
qflt coeff[max_deg + 1];
public:
inline poly(int);
qflt& operator [] (int);
qflt operator () (qflt&);
poly operator * (poly&);
poly operator + (poly&);
poly operator - ();
poly operator - (poly&);
friend ostream& operator << (ostream&, poly&);
friend istream& operator >> (istream&, poly&);
friend poly poly_deriv(poly&);
friend qflt& p_func::operator [] (int);
friend p_func::p_func(int, int);
friend qflt p_func::coeff_ptl(qflt&, int);
friend p_func p_deriv(p_func&);
friend void bld_nguess(pqf, pqf, p_func&, int, qflt[], qflt&);
friend int main(int, char *[]);
friend void bld_r_iguess(pqf, pqf, p_func&, int, qflt[], qflt&, qflt&);
};
/*
Construct a polynomial of given degree.
*/
inline poly::poly(int deg)
{
d = deg;
for (int i = 0; i <= d; i++)
coeff[i] = qflt_zero;
}
----------------- end poly.h ------------------
----------------- begin qflt.h ----------------
overload fabs;
#include <math.h>
double quiet_NaN(long); /* Should be in math.h but isn't */
#include <stream.h>
inline int min(int a, int b)
{
return a < b ? a : b;
}
inline int max(int a, int b)
{
return a > b ? a : b;
}
// Max with absolute value
inline double maxabs(double& a, double& b)
{
if (fabs(a) > fabs(b)) return a;
else return b;
}
// Min with absolute value
inline double minabs(double& a, double& b)
{
if (fabs(a) < fabs(b)) return a;
else return b;
}
class qflt
{
double high, low;
void split(double, double&, double&);
qflt exdblmul(double&, double&);
public:
qflt(double, double = 0.0);
inline qflt();
qflt operator + (qflt&);
qflt operator - (qflt&);
qflt operator - ();
qflt operator * (qflt&);
qflt operator / (qflt&);
inline void operator *= (qflt&);
inline int operator == (qflt&);
inline int operator != (qflt&);
inline int operator < (qflt&);
inline int operator > (qflt&);
inline int operator <= (qflt&);
inline int operator >= (qflt&);
inline double operator double();
friend ostream& operator << (ostream&, qflt&);
friend istream& operator >> (istream&, qflt&);
};
extern qflt qflt_zero;
extern qflt qflt_one;
inline qflt::qflt()
{
}
inline double qflt::operator double()
{
return(high + low);
}
inline void qflt::operator *= (qflt& a)
{
*this = *this * a;
}
inline int qflt::operator == (qflt& a)
{
return (high == a.high) && (low == a.low);
}
inline int qflt::operator != (qflt& a)
{
return (high != a.high) || (low != a.low);
}
inline int qflt::operator < (qflt& a)
{
return (high < a.high) || ((high == a.high) && (low < a.low));
}
inline int qflt::operator > (qflt& a)
{
return (high > a.high) || ((high == a.high) && (low > a.low));
}
inline int qflt::operator <= (qflt& a)
{
return (high < a.high) || ((high == a.high) && (low <= a.low));
}
inline int qflt::operator >= (qflt& a)
{
return (high > a.high) || ((high == a.high) && (low >= a.low));
}
inline qflt fabs(qflt& x)
{
return x < qflt_zero ? -x : x;
}
/* Declaration for pointer to qflt function */
typedef qflt (*pqf)(qflt&);
------------------- end qflt.h ---------------------
------------------- begin remes.h ------------------
const qflt macheps = qflt(1.0e-25, 0.0);
const int max_dstep = 5;
const int max_iter = 50;
------------------- end remes.h --------------------
I have not tried to figure out why this is happening...
-- Robert