possible glibc bug?
Hello,
I am a student taking a c programming class (I'm a fortran hacker) and
I believe I may have stumbled across a bug in libc for potato. I'm
attaching
my code. Please be kind as I am quite new at c. Anyway, if I run the
program,
tell it that I have 2 variables, t=0.0, y1 = 1.0, y2 = 10.0, RTOL =
1.0e-5,
ATOL = 1.0e-200, I have no problems, the program exits normally. If
I give it
the exact same inputs, but change RTOL to 1.0e-7, the program dies with
a
seg fault. I thought this was interesting, so I tried running it in gdb
and it said that
the program died like so:
Program received signal SIGSEGV, Segmentation fault.
0x2000011fb04 in chunk_free () from /lib/libc.so.6.1
I tried several things to clear this up, but couldn't figure it out, so
I tried using
ccc instead of gcc and got the same thing. The thing that led me to
think it
might be a bug in the alpha version of glibc is that I took the same
code and
tried compiling it and running on an x86 running woody. There were no
problems
when I tried that, I also attempted to compile it and run it using gcc
on a sparc running
solaris. Once again, I could not reproduce the same errors that I had
gotten on
my alpha.
Here is what information on gcc and the various libraries that I could
scrape together
for the different machines:
my alpha:
gcc: version 2.95.2 20000220 (Debian GNU/Linux)
libc: libc6.1_2.2-5_alpha.deb
the x86 box:
gcc version 2.95.2 20000220 (Debian GNU/Linux)
libc: glibc Version: 2.1.94-3
I can't get the c library version on the solaris box, but the gcc is:
gcc version 2.95 19990728 (release)
Like I said earlier, I'm quite new to c, so this is in all likelihood
just a programming
error of mine that for some reason only manifests itself on alpha.
Please be kind.
And don't bother sending me corrections on the homework, first, I know
this is not
the place for that and second, it was already due :-)
Luke Shulenburger
(sluke@mit.edu)
/*
* 10.001 Fall 2000
* Homework 12.0 solution
* Written by L. Shulenburger
* 12/1/00
* sluke@mit.edu
*
* This program uses the AM2 ode solver to take initial values of x[1] and x[2]
* from the user as well as the initial time and to calculate the future values
* for these variables given the functions exp_decay and exp_decayJJ.
*
*/
#include <stdio.h> /* needed for scanf, printf */
#include <stdlib.h> /* needed for malloc, calloc, and free */
#include <math.h> /* needed for exp and sqrt */
#define NR_END 1 /* needed for dmatrix */
#define FREE_ARG char* /* needed for free_dmatrix and free_dvector */
double **dmatrix(long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
double *dvector(long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void nrerror(char error_text[]);
int Newton_solve(void FJ(double X[],double F[],double **,double *,int *),int N, double X0[], double p[], int ip[], double RTOL, double ATOL, double X[]);
void gauss(double **a, double *b, double *x, int n);
void exp_decay(double Y[], double F[], double **J, double param[], int iparam[]);
void exp_decayJJ(double Y[], double FF[], double **JJ, double param[], int iparam[]);
void AM2(double y[], int neq, double h, double param[], int iparam[], double RTOL, double ATOL,
void derivs(double Y[], double F[], double **J, double param[], int iparam[]),
void minme(double Y[], double FF[], double **JJ, double param[], int iparam[]));
int main(void) {
int neq, i, j, z;
double RTOL, ATOL, *params, *save_params, *y, *yold, inith;
int *iparams, steps, initsteps, flag, nsteps;
FILE *ofp;
ofp = fopen("output", "w");
nsteps = 10;
neq = 4;
do { /* read in neq and make sure it is sane */
if (neq < 0) {
printf("A positive integer number of variables must be specified.\n");
}
printf("How many variables (excluding t) are in your model?\n");
scanf("%d", &neq);
} while (neq < 0);
/* allocate space for vectors */
params = calloc((neq+2), sizeof(double));
save_params = calloc((neq+2), sizeof(double));
y = calloc(neq, sizeof(double));
yold = calloc(neq, sizeof(double));
/* decrement pointers so y and yold will have 1 as offset */
y--;
yold--;
/* read in the current time */
printf("What is the initial value of the time in your model?\n");
scanf("%lf", ¶ms[0]);
/* read in initial values of all other variables */
for (i = 1; i <= neq; i++) {
printf("%s%d%s\n", "What is the initial value of variable ", i, " in your model?");
scanf("%lf", ¶ms[i]);
}
RTOL = 5.0;
do { /* read in ROTL and make sure it is attainable */
if (RTOL < 0) {
printf("A positive value of the relative error must be specified.\n");
} else if (RTOL < 1e-12) {
printf("That much accuracy cannot be attained with double precision\n");
printf("variables. Please choose a less restrictive error tolerance.\n");
}
printf("What is the maximum value for the relative error which is tolerable?\n");
scanf("%lf", &RTOL);
} while (RTOL < 1e-12);
ATOL = 5.0;
do { /* read in ROTL and make sure it is attainable */
if (ATOL < 0) {
printf("A positive value of the absolute error must be specified.\n");
} else if (ATOL < 1e-212) {
printf("That much accuracy cannot be attained with double precision\n");
printf("variables. Please choose a less restrictive error tolerance.\n");
}
printf("What is the maximum value for the absolute error which is tolerable?\n");
scanf("%lf", &ATOL);
} while (ATOL < 1e-212);
/* loop over all steps */
for (z = 1; z <= nsteps; z++) {
inith = params[neq+1] = 0.1;
initsteps = steps = 10;
for (i = 0; i <= neq; i++) {
save_params[i] = params[i];
}
/* set initial guess for the value of y at t = t+1 to
* the value of y at t = t */
for (i = 1; i <= neq; i++) {
yold[i] = params[i];
}
/* run AM2 over this interval with smaller and smaller h's until
* each of the y's has converged */
do {
if(steps == 1000000) {
printf("The requested accruacy cannot be achieved, sorry\n");
return 1;
}
/* reset parameters */
for (i = 0; i <=neq; i++) {
params[i] = save_params[i];
}
flag = 0;
/* get guess for value of y at t = t+1 */
for (i = 1; i <= steps; i++) {
AM2(y, neq, params[neq+1], params, iparams, RTOL, ATOL, exp_decay, exp_decayJJ);
for (j = 1; j <= neq; j++) {
params[j] = y[j];
}
}
/* check convergence criteria and then copy y into yold */
for (i = 1; i <= neq; i++) {
if (fabs((y[i] - yold[i])) >= fabs(y[i]) * RTOL + ATOL) {
flag++;
}
yold[i] = y[i];
}
params[neq+1] = params[neq+1] / 10.0;
steps = steps * 10;
} while(flag != 0);
/* write results to a file */
for (i = 0; i <= neq; i++) {
fprintf(ofp, "%.9g%s", params[i], " ");
}
fprintf(ofp, "\n");
}
/* free the space allocated for the pointers */
params++; save_params++;
y++;
yold++;
free(params);
free(save_params);
free(y);
free(yold);
fclose(ofp);
return 0;
}
/* this function uses the Adams Moulton second order method to take one step of
* size h in integrating differential equations. the new values are returned in
* y[]. The initial values are passed down in param[1] and param[2]..., The
* initial time is passed down in param[0]
*/
void AM2(double y[], int neq, double h, double param[], int iparam[], double RTOL, double ATOL,
void derivs(double Y[], double F[], double **J, double param[], int iparam[]),
void minme(double Y[], double FF[], double **JJ, double param[], int iparam[])) {
double *FE, *oldslope, **J, *BE, *newslope;
int i;
FE = calloc(neq, sizeof(double));
BE = calloc(neq, sizeof(double));
oldslope = calloc(neq, sizeof(double));
newslope = calloc(neq, sizeof(double));
FE--; BE--; oldslope--; newslope--;
J = dmatrix(1, neq, 1, neq);
/* copy initial guesses from params into y */
for (i = 1; i <= neq; i++) {
y[i] = param[i];
}
/* get the estimate of ynew from the forward euler method (put it in FE) */
derivs(y, oldslope, J, param, iparam);
for (i = 1; i <= neq; i++) {
FE[i] = y[i] + h * oldslope[i];
}
/* get slope from backward euler's method */
/* get value of y for evaluating derivative in BE method */
param[0] = param[0] + h;
Newton_solve(minme, neq, FE, param, iparam, 1.0e-10, 1.0e-200, BE);
derivs(BE, newslope, J, param, iparam);
/* calculate new value of y using AM2 */
for (i = 1; i <= neq; i++) {
y[i] = y[i] + (h / 2.0) * (newslope[i] + oldslope[i]);
}
FE++; oldslope++; BE++; newslope++;
free_dmatrix(J, 1, neq, 1, neq);
free(FE);
free(BE);
free(oldslope);
free(newslope);
return;
}
int Newton_solve(void FJ(double X[],double F[],double **,double *,int *),int N, double X0[], double p[], int ip[], double RTOL, double ATOL, double X[])
/* general multi-d Newton solver for F(X)=0.
FJ returns F(X) and its Jacobian J(X).
At output, X[] holds the solution
and Niter is the number of iterations used.
Terminates if Niter reaches 500. */
{
int Niter,i,converged;
double *F,*DX;
double **J;
J=dmatrix(1,N,1,N); /* allocate memory for Jacobian */
F=dvector(1,N); /* allocate memory for F (called "residual vector") */
DX=dvector(1,N);
for(i=1;i<=N;i++){X[i]=X0[i];} /* set X=X0 */
Niter=0;
do{
Niter=Niter+1;
converged=1;
FJ(X,F,J,p,ip); /* evaluate F(X) and J(X) */
for(i=1;i<=N;i++){F[i]=-F[i];} /* switch sign of F */
gauss(J,F,DX,N); /* solve for DX */
for(i=1;i<=N;i++){
X[i]=X[i]+DX[i]; /* update X */
/* test for failure to converge: */
if((fabs(DX[i]))>(RTOL*fabs(X[i])+ATOL))converged=-1;
}
}while((converged<0)&&(Niter<500));
/* if not converged yet and not too many iterations, try again */
free_dmatrix(J,1,N,1,N); /* free memory used for J */
free_dvector(F,1,N); /* free memory used for F */
free_dvector(DX,1,N); /* free memory used for DX */
return Niter;
}
void exp_decay(double Y[], double F[], double **J, double param[], int iparam[]) {
/* returns F(Y) and J(Y).With this F, both Y[1] and Y[2] decay exponentially from their initial values */
F[1] = -Y[1];
F[2] = -Y[2];
J[1][1] = -1;
J[2][2] = -1;
J[1][2] = J[2][1] = 0;
return;
}
void exp_decayJJ(double Y[], double FF[], double **JJ, double param[], int iparam[]) {
/* returns FF(Y) and JJ(Y) suitable for use with AM2 solver.
* param[1] = Yold[1]
* param[2] = Yold[2]
* param[3] = the time step h */
double h;
h=param[3];
FF[1] = param[1] - Y[1] - (0.5*h*(param[1]+Y[1]));
FF[2] = param[2] - Y[2] - (0.5*h*(param[2]+Y[2]));
JJ[1][1] = -1 - 0.5*h;
JJ[2][2] = -1 - 0.5*h;
JJ[1][2] = JJ[2][1]= 0;
return;
}
void gauss(double **a, double *b, double *x, int n) {
/* This program uses Gaussian Elimination with
* pivoting to solve the problem A x =b.
* Use is made of array offsets so that the indices go from 1 to n.
*
* Author(s): R. Sureshkumar (10 January 1997)
* Modified by: Gregory J. McRae (22 October 1997)
*
* Usage: gauss(a,b,x,n)
*
* a - Matrix a[1..n][1..n]
* b - Right hand side vector b[1..n]
* x - Desired solution vector
* n - Matrix dimensions
*/
int i,j,k,m,rowx;
double xfac,temp,temp1,amax;
/*_______________________________________
Do the forward reduction step.
_______________________________________*/
rowx = 0; /* Keep count of the row interchanges */
for (k=1; k<=n-1; ++k) {
amax = (double) fabs(a[k][k]) ;
m = k;
for (i=k+1; i<=n; i++){ /* Find the row with largest pivot */
xfac = (double) fabs(a[i][k]);
if(xfac > amax) {amax = xfac; m=i;}
}
if(m != k) { /* Row interchanges */
rowx = rowx+1;
temp1 = b[k];
b[k] = b[m];
b[m] = temp1;
for(j=k; j<=n; j++) {
temp = a[k][j];
a[k][j] = a[m][j];
a[m][j] = temp;
}
}
for (i=k+1; i<=n; ++i) {
xfac = a[i][k]/a[k][k];
for (j=k+1; j<=n; ++j) {
a[i][j] = a[i][j]-xfac*a[k][j];
}
b[i] = b[i]-xfac*b[k];
}
}
/*_______________________________________
Do the back substitution step
_______________________________________*/
for (j=1; j<=n; ++j) {
k=n-j+1;
x[k] = b[k];
for(i=k+1; i<=n; ++i) {
x[k] = x[k]-a[k][i]*x[i];
}
x[k] = x[k]/a[k][k];
}
return;
}
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
/* sample output:
[sluke@bartender] ~/10.001/ps12/> hw11_0.x
How many variables (excluding t) are in your model?
2
What is the initial value of the time in your model?
0.0
What is the initial value of variable 1 in your model?
1.0
What is the initial value of variable 2 in your model?
10.0
What is the maximum value for the relative error which is tolerable?
1.0e-5
What is the maximum value for the absolute error which is tolerable?
1.0e-200
This generates the file output which looks like this:
1 0.367879411 3.67879411
2 0.135335261 1.35335261
3 0.0497870559 0.497870559
4 0.0183156328 0.183156328
5 0.00673794419 0.0673794419
6 0.00247875094 0.0247875094
7 0.000911881434 0.00911881434
8 0.000335462404 0.00335462404
9 0.000123409712 0.00123409712
10 4.53998919e-05 0.000453998919
*/
Reply to: