Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Error freeing memory....

6 views
Skip to first unread message

Sheldon

unread,
Mar 14, 2010, 1:57:03 PM3/14/10
to
Hi everyone,

I have a fairly large program that has a problem that I cannot seem to
solve.
The program runs smoothly but when it tries to free allocated memory
before exiting I get a segmentation fault.
I have checked many times the following:

1. That I have not freed the memory before hand and call free this
way: if(array) free(array);

2. The arrays are vectors that was allocated with error checking:
if(mem == NULL)...

3. There are no conflicts like a double definition of a variable.

4. Variables are assigned in structs

But still there I get his error. Can anyone offer some advice on what
may be the problem?
The code is nearly 1000 lines long. I can post it if anyone wants to
look at it or I can send it via email.
I would greatly appreciate any help with this problem.


/S

Seebs

unread,
Mar 14, 2010, 2:09:33 PM3/14/10
to

A thousand lines isn't all that much. A couple thoughts:

1. Go through removing functional code to strip it down. If you find
that, after removing something, you suddenly don't have a problem, then
that code is what did it.
2. Once you get it stripped down, it should be a lot shorter.

My guess would be that you're overrunning something -- allocating the
wrong size of space. The last time I had this, I'd allocated enough space
for a pointer, rather than enough space for the thing pointed to...

-s
--
Copyright 2010, all wrongs reversed. Peter Seebach / usenet...@seebs.net
http://www.seebs.net/log/ <-- lawsuits, religion, and funny pictures
http://en.wikipedia.org/wiki/Fair_Game_(Scientology) <-- get educated!

Eric Sosman

unread,
Mar 14, 2010, 2:53:21 PM3/14/10
to
On 3/14/2010 1:57 PM, Sheldon wrote:
> Hi everyone,
>
> I have a fairly large program that has a problem that I cannot seem to
> solve.
> The program runs smoothly but when it tries to free allocated memory
> before exiting I get a segmentation fault.
> I have checked many times the following:
>
> 1. That I have not freed the memory before hand and call free this
> way: if(array) free(array);

A note in passing: This makes no difference, because
free(NULL) is legal (and does nothing).

> 2. The arrays are vectors that was allocated with error checking:
> if(mem == NULL)...
>
> 3. There are no conflicts like a double definition of a variable.
>
> 4. Variables are assigned in structs
>
> But still there I get his error. Can anyone offer some advice on what
> may be the problem?
> The code is nearly 1000 lines long. I can post it if anyone wants to
> look at it or I can send it via email.
> I would greatly appreciate any help with this problem.

Almost certainly, you have somehow corrupted the "private"
data with which malloc() and free() and so on keep track of what
is allocated and what isn't. Unfortunately, there is a very
large number of ways this corruption could have occurred, and
a lot of them will leave no visible traces until long after the
damage has been done (so, for example, the faulty function may
not appear in a stack trace when the crash occurs). These can
be very difficult to track down, but there are a few patterns
that occur often enough that they're worth checking for:

- Calling free() on a pointer to memory that wasn't obtained
from malloc() (include realloc() and calloc() as appropriate, here
and throughout):

char array[] = "Hello, world!";
...
free (array); /* Bad: not dynamic memory */

- Calling free() with a pointer to a block of dynamic memory,
but not to its start:

char *ptr = malloc(42);
...
++ptr;
...
free (ptr); /* Bad: not the original value */

- Calling free() twice on the same block:

char *ptr = malloc(42);
...
free (ptr); /* No problem */
...
free (ptr); /* Bad: not "your" memory any more */

- Storing outside the bounds of an allocated area, and thus
stomping on whatever is next to it. The classic example in C is:

char *ptr = malloc(strlen(string));
strcpy (ptr, string); /* Bad: no space for the '\0' */

- Using an uninitialized or otherwise invalid pointer or
array index. Several examples:

char *ptr;
...
strcpy (ptr, "Hello"); /* Bad: ptr "points to garbage" */

char *ptr;
...
free (ptr); /* Bad: ptr "points to garbage" */

char array[42];
int idx;
...
array[idx] = 'X'; /* Bad: idx "is garbage" */

If inspection doesn't reveal any of these (or similar) bugs
in your code, you'll have to resort to run-time checking. Tools
to assist with this are platform-specific: Linux has valgrind,
Solaris has libumem, things like dbmalloc and Electric Fence are
fairly widely available and you may be able to get them running
on your system, and there are commercial products like Purify.
The very existence of all these debugging tools shows that the
problem is (1) not uncommon and (2) uncommonly difficult to solve
unaided. Good luck!

--
Eric Sosman
eso...@ieee-dot-org.invalid

Sheldon

unread,
Mar 14, 2010, 3:46:50 PM3/14/10
to
> esos...@ieee-dot-org.invalid

Thanks guys! Both your answers are very good and insightful.
I'm on my way home and will look at this deeper. But I'll post the
code just in case your trained eyes sees something I've missed.
Please remember that I'm not a programmer. I learned C on my own and
mostly via trail and error. I need this program to aid in my studies.
To do this in MATLAB, PYTHON, etc would take too long time and too
much memory. I would love to learn more and to be a better programmer.
Your help is truly appreciated. :-))


Header file:
*********
/*********************************************************
* This program reads a directory for GRIB files. *
* *
* *
* Functions: Uses buildin functions from stdlib and *
* dirent library as well NETCDF, and *
* GRIBEX library. *
* *
* See gribread.h for specifics *
* *
*********************************************************/

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <dirent.h>
#include <errno.h>
#include <ctype.h>
#include <stdarg.h>
#include <assert.h>
#include <unistd.h>
#include <math.h>
#include <netcdf.h>
#include <cdi.h>
#include "utilities.h"

// Macros
#define ERRCODE 2
#define ERR(e) {printf("Error: %s\n", nc_strerror(e)); exit(ERRCODE);}

#define STR_LENGTH 500
#define STR_NAME 128
#define G0 9.80665 // m s-2
#define EQRAD 6378.14e3 // meters
#define PI 3.14159265
#define FLAT 1.0/298.257
#define Boltzmann 1.380658e-23
#define Avogadro 6.0221367e23
#define UnGascon Boltzmann * Avogadro
#define Md 28.964
#define Mv 18.016
#define RD 1000.0 * UnGascon/Md
#define RV 1000.0 * UnGascon/Mv
#define CPD 3.5 * RD
#define CVD CPD - RD
#define CPV 4.0 * RV
#define TS 3600 // time step for the model run
#define FOD 10000
#define TIMESTEP 3*FOD

#define PACKED_TYPE 16
#define SCALE "scale_factor"
#define OFFSET "add_offset"
#define TITLE "EC EARTH Comparison Data Radio and Space (GEM)"
#define TTLE "TITLE"
#define LOCATION "Institution"
#define INST "Radio and Space Science Chalmers Technical University
Gothenburg"
#define OPATH "/home/marston/ecearth_data/v21/t159l62/mmh_"

#define DEGREES_EAST "degrees_east"
#define DEGREES_NORTH "degrees_north"
#define LAT_NAME "latitude"
#define LON_NAME "longitude"
#define NREC 8
#define REC_NAME "time"
#define LVL_NAME "level"

// Names of things
#define PRES_NAME "pressure"
#define TEMP_NAME "temperature"
#define TTR_NAME "Top of the atmosphere thermal radiation"
#define UNITS "units"

#define PRES_UNITS "Pa"
#define TEMP_UNITS "K"
#define WP_UNITS "g m-2"
#define IWC_UNITS "g m-3"
#define MAX_ATT_LEN 80

#define TEMP "T"
#define SFCP "SP"
#define SPH "Q"
#define OLR "TTR"
#define ELEV "Z"
#define SSN "var100"
#define CSN "var101"
#define CIW "CIWC"

struct V {
char date[NC_MAX_NAME+1];
char pdate[NC_MAX_NAME+1];;
size_t ilev;
size_t mlev;
size_t nlon;
size_t nlat;
size_t nrec;
size_t pnrec;
size_t nlev;
size_t ndims;
char vname[NC_MAX_NAME+1];
char pvname[NC_MAX_NAME+1];
char lname[NC_MAX_NAME+1];
char units[NC_MAX_NAME+1];
int varid;
int pvarid;
size_t size;
size_t gridsize;
char infile[NC_MAX_NAME+1];
char pinfile[NC_MAX_NAME+1];
float *mean;
};

struct AR {
float *gph;
double *lats;
float *ph;
float *hgph;
double *lon;
double *lat;
double *a;
double *b;
double *time;
float *FP;
float *GMH;
float *RH;
float *CIWP;
float *TCIWP;
};

struct W {
int rec;
int ncid;
int lon_dimid;
int lat_dimid;
int lvl_dimid;
int rec_dimid;
int latid;
int lonid;
int lvlid;
int fpid;
int tid;
int ciid;
int csid;
int ssid;
int ttrid;
int rhid;
int iwpid;
int tiwpid;
int gmhid;
int dimids4[4];
int dimids3[3];
size_t start4[4];
size_t count4[4];
size_t start3[3];
size_t count3[3];
};

double mlvls[62] =
{62,61,60,59,58,57,56,55,54,53,52,51,50,49,48,47,46,45,44,\
43,42,41,40,39,38,37,36,35,34,33,32,31,30,29,28,27,26,25,\
24,23,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1};

struct dirent **namelist;
char FILESTR[NC_MAX_NAME+1];
char filter1[NC_MAX_NAME+1];
char filter2[NC_MAX_NAME+1];
char filter3[NC_MAX_NAME+1];
/*
Function Prototypes
*/
int read_dir(char *path);
int file_select(const struct dirent *entry);
void USAGE(int nargs, char **args);
void handle_error(int status);
float GeoidRadius(float lat);
void FullPresGph(struct V *T,
struct V *Q,
struct V *Z,
struct V *SP,
struct AR *A);
float EI(float T);
float EW(float T);
void RHRho(struct V *Q,
struct V *T,
struct V *SS,
struct V *CS,
struct V *CI,
struct AR *A);
void IW_PATH(struct V *CS,
struct V *SS,
struct V *CI,
struct AR *A);
void GetInfo(struct V *info, char *name, char *date, char *pdate, char
*path, int nlev);
void GetXtra(struct AR *A, char *infile);
void A1DTo2D(struct V *T, float *a1d, float **a2d);
void A1DTo3D(struct V *T, float *a1d, float ***a3d);
void GetArDims(int ncid, int ndims, struct V *p, int *dimids, size_t
*array_lenp);
void GetDims(int ncid, int ndims, char *vname, struct V *p);
void GetAtt(int ncid, int natts, struct V *p);
void GetMean(struct AR *A, struct V *ret, int time, double *array);
void SetupWrite(struct V *T, struct AR *A, struct W *w);
void WRITE(struct W *w, float ***ar3D, float **ar2D, int id);

********

Main program:
*****
#include "gribread.h"

int main(int argc, char *argv[]) {

int i, ONCE = 0, FIRST = atoi(argv[4]);
char path[NC_MAX_NAME+1], buf[NC_MAX_NAME+1];
char date[NC_MAX_NAME+1], pdate[NC_MAX_NAME+1];
char FILE_NAME[NC_MAX_NAME+1];
struct V *T = NULL, *SP = NULL, *Q = NULL, *TTR = NULL, *Z = NULL,
*SS = NULL, *CS = NULL, *CI = NULL;
struct AR *A = NULL;
struct W *w;

double *array2d;
double *array3d;
float ***ar3D;
float **ar2D;

A = (struct AR *)smalloc(sizeof(struct AR));
T = (struct V *)smalloc(sizeof(struct V));
SP = (struct V *)smalloc(sizeof(struct V));
Q = (struct V *)smalloc(sizeof(struct V));
TTR = (struct V *)smalloc(sizeof(struct V));
Z = (struct V *)smalloc(sizeof(struct V));
SS = (struct V *)smalloc(sizeof(struct V));
CS = (struct V *)smalloc(sizeof(struct V));
CI = (struct V *)smalloc(sizeof(struct V));
w = (struct W *)smalloc(sizeof(struct W));

// Loop indexes.
int rec;
// Error handling.
int retval;

if(argc != 5) {
USAGE(argc,argv);
printf("argc = %d\n",argc);
}

printf("INFO: Calling program %s for files %s in dir:
\n",argv[0],argv[1]);

memset(path,'\0',sizeof(char)*(NC_MAX_NAME+1));
strcpy(FILESTR,argv[2]); // Experiment name
strcpy(path,argv[1]); // directory
memset(filter1,'\0',sizeof(char)*(NC_MAX_NAME+1));
strcpy(filter1,".NC"); // file type

if(path[strlen(path)-1] != '/') strcat(path,"/");
printf("INFO: path= %s\n",path);

// Read directory for only NETCDF files and sorting them

for(i = 0; i < 12; i++) {

memset(date,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(buf,'\0',sizeof(char)*(NC_MAX_NAME+1));
if(i == 0) {
strcpy(date,argv[3]);
sprintf(buf,"%02d",i+1);
strcat(date,buf);
if(FIRST == 0) {
strcpy(pdate,"000000");
} else {
sprintf(buf,"%04d",(atoi(argv[3])-1));
strcpy(pdate,buf);
sprintf(buf,"%02d",12);
strcat(pdate,buf);
}
} else {
strcpy(date,argv[3]);
strcpy(pdate,argv[3]);
sprintf(buf,"%02d",i+1);
strcat(date,buf);
sprintf(buf,"%02d",i);
strcat(pdate,buf);
}

memset(filter2,'\0',sizeof(char)*(NC_MAX_NAME+1));
strcpy(filter2,"SH");
GetInfo(T,TEMP,date,pdate,path,1);
GetInfo(SP,SFCP,date,pdate,path,0);
strcpy(filter2,"GG");
GetInfo(Q,SPH,date,pdate,path,1);
GetInfo(TTR,OLR,date,pdate,path,0);
strcpy(filter2,"SH");
GetInfo(Z,ELEV,date,pdate,path,0);
strcpy(filter2,"GG");
GetInfo(SS,SSN,date,pdate,path,1);
GetInfo(CS,CSN,date,pdate,path,1);
GetInfo(CI,CIW,date,pdate,path,1);

if(ONCE == 0) {
ar2D = ConMem2DF(T->nlat,T->nlon);
ar3D = ConMem3DF(T->nlev,T->nlat,T->nlon);
A->gph = smalloc(sizeof(float)*T->size);
A->hgph = smalloc(sizeof(float)*T->ilev*T->gridsize);
A->ph = smalloc(sizeof(float)*T->ilev*T->gridsize);
A->lats = smalloc(sizeof(double)*T->gridsize);
A->lon = smalloc(sizeof(double)*T->nlon);
A->lat = smalloc(sizeof(double)*T->nlat);
A->a = smalloc(sizeof(double)*T->ilev);
A->time = smalloc(sizeof(double)*T->nrec);
A->b = smalloc(sizeof(double)*(T->ilev));
array2d = smalloc(sizeof(double)*SP->gridsize);
array3d = smalloc(sizeof(double)*T->size);
T->mean = smalloc(sizeof(float)*T->size);
Q->mean = smalloc(sizeof(float)*Q->size);
SP->mean = smalloc(sizeof(float)*SP->size);
Z->mean = smalloc(sizeof(float)*Z->size);
SS->mean = smalloc(sizeof(float)*SS->size);
TTR->mean = smalloc(sizeof(float)*TTR->size);
CS->mean = smalloc(sizeof(float)*CS->size);
CI->mean = smalloc(sizeof(float)*CI->size);
A->FP = smalloc(sizeof(float)*T->size);
A->GMH = smalloc(sizeof(float)*T->size);
A->RH = smalloc(sizeof(float)*T->size);
A->CIWP = smalloc(sizeof(float)*T->gridsize);
A->TCIWP = smalloc(sizeof(float)*T->gridsize);
ONCE = 1;
GetXtra(A,T->infile);
}
memset(FILE_NAME,'\0',sizeof(char)*(NC_MAX_NAME+1));
strcpy(FILE_NAME,OPATH);
strcat(FILE_NAME,date);
strcat(FILE_NAME,"_");
strcat(FILE_NAME,FILESTR);
strcat(FILE_NAME,"_");
strcat(FILE_NAME,"monmean.nc");
printf("INFO: OUTFILE: %s\n",FILE_NAME);

// Create the file.
if((retval = nc_create(FILE_NAME,NC_64BIT_OFFSET | NC_CLOBBER,&w-
>ncid))) ERR(retval);

/* These settings tell netcdf to write one timestep of data. (The
setting of start[0] inside the loop below tells netCDF which
timestep to write.) */
w->count4[0] = 1;
w->count4[1] = T->nlev;
w->count4[2] = T->nlat;
w->count4[3] = T->nlon;
w->start4[1] = 0;
w->start4[2] = 0;
w->start4[3] = 0;

w->count3[0] = 1;
w->count3[1] = T->nlat;
w->count3[2] = T->nlon;
w->start3[1] = 0;
w->start3[2] = 0;

SetupWrite(T,A,w);

// This will write the data. The arrays only hold one timestep
worth
// of data. The data would change between timesteps.
for(rec = 0; rec < NREC; rec++) {

printf("INFO: Fetching timestep %d\n",rec*TIMESTEP);

printf("INFO: Temp\n");
GetMean(A,T,rec*TIMESTEP,array3d);
MMF(T->mean,T->size);
printf("INFO: Q\n");
GetMean(A,Q,rec*TIMESTEP,array3d);
MMF(Q->mean,Q->size);
printf("INFO: SP\n");
GetMean(A,SP,rec*TIMESTEP,array2d);
MMF(SP->mean,SP->size);
printf("INFO: Z\n");
GetMean(A,Z,rec*TIMESTEP,array2d);
MMF(Z->mean,Z->size);
printf("INFO: Calculate FP and GPH\n");
FullPresGph(T,Q,Z,SP,A);
MMF(A->FP,T->size);
MMF(A->GMH,T->size);
printf("INFO: SS\n");
GetMean(A,SS,rec*TIMESTEP,array3d);
MMF(SS->mean,SS->size);
printf("INFO: CS\n");
GetMean(A,CS,rec*TIMESTEP,array3d);
MMF(CS->mean,CS->size);
printf("INFO: CI\n");
GetMean(A,CI,rec*TIMESTEP,array3d);
MMF(CI->mean,CI->size);
printf("INFO: Calculate RH and convert to g m3\n");
RHRho(Q,T,SS,CS,CI,A);
MMF(A->RH,T->size);
printf("INFO: Calculate IWP\n");
IW_PATH(CS,SS,CI,A);
printf("INFO: TTR\n");
GetMean(A,TTR,rec*TIMESTEP,array2d);
MMF(TTR->mean,TTR->size);

printf("INFO: Writing variables...\n");

w->rec = rec;

A1DTo3D(T,A->FP,ar3D);
WRITE(w,ar3D,0,w->fpid);

A1DTo3D(T,A->GMH,ar3D);
WRITE(w,ar3D,0,w->gmhid);

A1DTo3D(T,A->RH,ar3D);
WRITE(w,ar3D,0,w->rhid);

A1DTo3D(T,T->mean,ar3D);
WRITE(w,ar3D,0,w->tid);

A1DTo3D(T,CS->mean,ar3D);
WRITE(w,ar3D,0,w->csid);

A1DTo3D(T,SS->mean,ar3D);
WRITE(w,ar3D,0,w->ssid);

A1DTo3D(T,CI->mean,ar3D);
WRITE(w,ar3D,0,w->ciid);

A1DTo2D(T,TTR->mean,ar2D);
WRITE(w,0,ar2D,w->ttrid);

A1DTo2D(T,A->CIWP,ar2D);
WRITE(w,0,ar2D,w->iwpid);

A1DTo2D(T,A->TCIWP,ar2D);
WRITE(w,0,ar2D,w->tiwpid);

printf("INFO: Completed time averaging of time steps...\n");
} // end loop over time steps
Print("Closing output file...");
if((retval = nc_close(w->ncid))) ERR(retval);
//if(retval != NC_NOERR) handle_error(retval);
printf("INFO: *** SUCCESS writing file %s!\n", FILE_NAME);
}
printf("INFO: Freeing arrays");
if(ar2D) FreeCon2DF(ar2D);
if(ar3D) FreeCon3DF(ar3D);
if(A->a) free(A->a);
if(A->b) free(A->b);
if(A->lat) free(A->lat);
if(A->lon) free(A->lon);
if(A->gph) free(A->gph);
if(A->hgph) free(A->hgph);
if(A->ph) free(A->ph);
if(A->lats) free(A->lats);

if(SS->mean) free(SS->mean);
if(CS->mean) free(CS->mean);
if(CI->mean) free(CI->mean);
if(TTR->mean) free(TTR->mean);
if(T->mean) free(T->mean);
if(Z->mean) free(Z->mean);
if(Q->mean) free(Q->mean);
if(SP->mean) free(SP->mean);
if(A->FP) free(A->FP);
if(A->GMH) free(A->GMH);
if(A->CIWP) free(A->CIWP);
if(A->TCIWP) free(A->TCIWP);
if(A->RH) free(A->RH);
if(A->gph) free(A->gph);
if(A->ph) free(A->ph);
if(A->hgph) free(A->hgph);
if(A->lats) free(A->lats);

if(T) free(T);
if(TTR) free(TTR);
if(SS) free(SS);
if(CS) free(CS);
if(CI) free(CI);
if(Z) free(Z);
if(Q) free(Q);
if(SP) free(SP);
if(A) free(A);
if(array2d) free(array2d);
if(array3d) free(array3d);
if(namelist) free(namelist);
return EXIT_SUCCESS;
}

void WRITE(struct W *w, float ***ar3D, float **ar2D, int id) {

int retval;

w->start4[0] = w->rec;
w->start3[0] = w->rec;

if(ar3D) {
if((retval = nc_put_vara_float(w->ncid,id,w->start4,w-
>count4,&ar3D[0][0][0]))) ERR(retval);
}
if(ar2D) {
if((retval = nc_put_vara_float(w->ncid,id,w->start3,w-
>count3,&ar2D[0][0]))) ERR(retval);
}
return;
}

void SetupWrite(struct V *T, struct AR *A, struct W *w) {

int retval;

/* Define the dimensions. The record dimension is defined to have
* unlimited length - it can grow as needed. In this example it is
* the time dimension.*/
if((retval = nc_def_dim(w->ncid,LVL_NAME,T->mlev,&w->lvl_dimid)))
ERR(retval);
//if(retval != NC_NOERR) handle_error(retval);

if((retval = nc_def_dim(w->ncid,LAT_NAME,T->nlat,&w->lat_dimid)))
ERR(retval);

if((retval = nc_def_dim(w->ncid,LON_NAME,T->nlon,&w->lon_dimid)))
ERR(retval);

if((retval = nc_def_dim(w->ncid,REC_NAME,NC_UNLIMITED,&w-
>rec_dimid))) ERR(retval);

if((retval = nc_def_var(w->ncid,LAT_NAME,NC_DOUBLE,1,&w-
>lat_dimid,&w->latid))) ERR(retval);

if((retval = nc_def_var(w->ncid,LON_NAME,NC_DOUBLE,1,&w-
>lon_dimid,&w->lonid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"Levels",NC_DOUBLE,1,&w-
>lvl_dimid,&w->lvlid))) ERR(retval);

/* Assign units attributes to coordinate variables. */
if((retval = nc_put_att_text(w->ncid,w-
>latid,UNITS,strlen(DEGREES_NORTH),DEGREES_NORTH))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w-
>lonid,UNITS,strlen(DEGREES_EAST),DEGREES_EAST))) ERR(retval);

w->dimids4[0] = w->rec_dimid;
w->dimids4[1] = w->lvl_dimid;
w->dimids4[2] = w->lat_dimid;
w->dimids4[3] = w->lon_dimid;

w->dimids3[0] = w->rec_dimid;
w->dimids3[1] = w->lat_dimid;
w->dimids3[2] = w->lon_dimid;

// Define the netCDF variables

if((retval = nc_def_var(w->ncid,"P",NC_FLOAT,4,w->dimids4,&w-
>fpid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"T",NC_FLOAT,4,w->dimids4,&w->tid)))
ERR(retval);

if((retval = nc_def_var(w->ncid,"RH",NC_FLOAT,4,w->dimids4,&w-
>rhid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"GMH",NC_FLOAT,4,w->dimids4,&w-
>gmhid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"CS",NC_FLOAT,4,w->dimids4,&w-
>csid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"CI",NC_FLOAT,4,w->dimids4,&w-
>ciid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"SS",NC_FLOAT,4,w->dimids4,&w-
>ssid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"TTR",NC_FLOAT,3,w->dimids3,&w-
>ttrid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"CIWP",NC_FLOAT,3,w->dimids3,&w-
>iwpid))) ERR(retval);

if((retval = nc_def_var(w->ncid,"TCIWP",NC_FLOAT,3,w->dimids3,&w-
>tiwpid))) ERR(retval);

// Assign units attributes to the netCDF variables.
if((retval = nc_put_att_text(w->ncid,w-
>fpid,UNITS,strlen(PRES_UNITS),PRES_UNITS))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w-
>tid,UNITS,strlen(TEMP_UNITS),TEMP_UNITS))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w->csid,UNITS,strlen("g
m-3"),"g m-3"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w->ssid,UNITS,strlen("g
m-3"),"g m-3"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w->ciid,UNITS,strlen("g
m-3"),"g m-3"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w-
>gmhid,UNITS,strlen("m"),"m"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w-
>rhid,UNITS,strlen("%"),"%"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w->iwpid,UNITS,strlen("g
m-2"),"g m-2"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w->tiwpid,UNITS,strlen("g
m-2"),"g m-2"))) ERR(retval);

if((retval = nc_put_att_text(w->ncid,w->ttrid,UNITS,strlen("W
m-2"),"W m-2"))) ERR(retval);

// End define mode.
if((retval = nc_enddef(w->ncid))) ERR(retval);


// Write the coordinate variable data. This will put the latitudes
// and longitudes of our data grid into the netCDF file.
if((retval = nc_put_var_double(w->ncid,w->latid,&A->lat[0])))
ERR(retval);

if((retval = nc_put_var_double(w->ncid,w->lonid,&A->lon[0])))
ERR(retval);

if((retval = nc_put_var_double(w->ncid,w->lvlid,&mlvls[0])))
ERR(retval);

return;
}


void GetInfo(struct V *info, char *name, char *date, char *pdate, char
*path, int nlev) {

int status = 0, j = 0, k = 0;
int ncid = 0, ndims = 0,dimids[NC_MAX_VAR_DIMS] = {0},natts = 0;
int ndimsp = 0, nvarsp = 0, ngattsp = 0, unlimdimidp = 0;
char getname[NC_MAX_NAME+1];
int n = 0;
size_t len = 0;

memset(getname,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->date,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->pdate,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->vname,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->lname,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->units,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->infile,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(info->pinfile,'\0',sizeof(char)*(NC_MAX_NAME+1));
memset(filter3,'\0',sizeof(char)*(NC_MAX_NAME+1));

info->ilev = -1;
info->mlev = -1;
info->nlon = -1;
info->nlat = -1;
info->nrec = -1;
info->ndims = -1;
info->varid = -1;
info->size = -1;
info->gridsize = -1;

strcpy(info->date,date);
strcpy(filter3,info->date);
n = read_dir(path);
if(n != 1) {
fprintf(stderr,"Wrong number of files: (%d) found!",n);
nrerror("Should only find 1 file!");
}
strcpy(info->infile,path);
strcat(info->infile,namelist[0]->d_name);

// finding the previous file
strcpy(info->pdate,pdate);
strcpy(filter3,info->pdate);
n = read_dir(path);
if(n != 1) {
fprintf(stderr,"Wrong number of files: (%d) found!",n);
nrerror("Should only find 1 file!");
}
strcpy(info->pinfile,path);
strcat(info->pinfile,namelist[0]->d_name);

printf("INFO: Opening file ->%s\n",info->infile);

status = nc_open(info->infile, NC_NOWRITE, &ncid);
if (status != NC_NOERR) handle_error(status);

// Inquire what's in one file and extract meta data
status = nc_inq(ncid,&ndimsp,&nvarsp,&ngattsp,&unlimdimidp);
if(status != NC_NOERR) handle_error(status);
for(j = 0; j < nvarsp; j++) {
status = nc_inq_var(ncid,j,getname,0,&ndims,dimids,&natts);
if(status != NC_NOERR) handle_error(status);
//printf("Trying to match %s -> %s\n",name,getname);
if(strcmp(getname,name) == 0) {
printf("Matched %s -> %s and ID %d\n",name,getname,j);
strcpy(info->vname,getname);
info->varid = j;
break;
} else {
continue;
}
}
info->ndims = ndims-1;
// Loop over the the number of dimensions
GetDims(ncid,ndimsp,getname,info);
//GetArDims(ncid,ndims,info,dimids,array_lenp);
if(nlev == 1) info->nlev = info->mlev; else info->nlev = 1;
info->size = info->nlat*info->nlon*info->nlev;
info->gridsize = info->nlat*info->nlon;
printf("Nlev = %ld size = %ld\n",info->nlev,info->size);
GetAtt(ncid,natts,info);
status = nc_close(ncid);
if(status != NC_NOERR) handle_error(status);


// Open previous file
printf("Opening previous file %s\n",info->pinfile);
status = nc_open(info->pinfile, NC_NOWRITE, &ncid);
if (status != NC_NOERR) handle_error(status);

status = nc_inq(ncid,&ndimsp,&nvarsp,&ngattsp,&unlimdimidp);
if(status != NC_NOERR) handle_error(status);
for(j = 0; j < nvarsp; j++) {
status = nc_inq_var(ncid,j,getname,0,&ndims,dimids,&natts);
if(status != NC_NOERR) handle_error(status);
//printf("Trying to match %s -> %s\n",name,getname);
if(strcmp(getname,name) == 0) {
printf("Matched %s -> %s and ID %d\n",name,getname, j);
strcpy(info->pvname,getname);
info->pvarid = j;
break;
} else {
continue;
}
}

// Loop over the the number of dimensions
for(k = 0; k < ndimsp; k++) {
status = nc_inq_dim(ncid,k,getname,&len);
if(status != NC_NOERR) handle_error(status);
if(strcmp(getname,"time") == 0) info->pnrec = len;
}
printf("Nrec = %ld\n",info->pnrec);
status = nc_close(ncid);
if(status != NC_NOERR) handle_error(status);
Print("Complete..");
return;
}

void GetAtt(int ncid, int natts, struct V *p) {

int k = 0, status = 0;
char name[NC_MAX_NAME+1];
memset(name,'\0',sizeof(char)*(NC_MAX_NAME+1));

if(strcmp(p->vname,"var100") == 0) {
strcpy(p->lname,"Stratiform snow on level");
strcpy(p->units,"kg kg**-1");
printf("INFO: Lname = %s Units = %s\n",p->lname,p->units);
return;
} else if(strcmp(p->vname,"var101") == 0) {
strcpy(p->lname,"Convective snow on level");
strcpy(p->units,"kg kg**-1");
printf("INFO: Lname = %s Units = %s\n",p->lname,p->units);
return;
}

for(k = 0; k < natts; k++) {
//Print("Get the attributes of the variables");
//Print(p->vname);
status = nc_inq_attname(ncid,p->varid,k,name);
if(status != NC_NOERR) handle_error(status);
if(strcmp(name,"long_name") == 0) {
status = nc_get_att_text(ncid,p->varid,name,p->lname);
if(status != NC_NOERR) handle_error(status);
printf("INFO: Lname = %s\n",p->lname);
}
if(strcmp(name,"units") == 0) {
status = nc_get_att_text(ncid,p->varid,name,p->units);
if(status != NC_NOERR) handle_error(status);
printf("INFO: Units = %s\n",p->units);
}
}
}

void GetDims(int ncid, int ndims, char *vname, struct V *p) {

int k = 0, status = 0;
size_t len = 0;

for(k = 0; k < ndims; k++) {
status = nc_inq_dim(ncid,k,vname,&len);
if(status != NC_NOERR) handle_error(status);
if(strcmp(vname,"lon") == 0) p->nlon = len;
if(strcmp(vname,"lat") == 0) p->nlat = len;
if(strcmp(vname,"mlev") == 0) p->mlev = len;
if(strcmp(vname,"ilev") == 0) p->ilev = len;
if(strcmp(vname,"time") == 0) p->nrec = len;
}
printf("INFO: Nlon = %ld Nlat = %ld Nrec = %ld ilev = %ld mlev = %ld
\n",p->nlon,p->nlat,p->nrec,p->ilev,p->mlev);
return;
}

/*
void GetArDims(int ncid, int ndims, struct V *p, int *dimids, size_t
*array_lenp) {

int k, status;
char name[NC_MAX_NAME+1];

// Get the dimensions for the arrays
for(k = 0; k < ndims; k++) {
status = nc_inq_dimlen(ncid,dimids[k],&array_lenp[k]);
if(status != NC_NOERR) handle_error(status);
status = nc_inq_dimname(ncid,dimids[k],name);
if(status != NC_NOERR) handle_error(status);
printf("%s = %ld\n",name,(int)array_lenp[k]);
if(ndims == 4) {
p->nlev = p->mlev;
p->size = p->nlat*p->nlon*p->nlev;
} else if(ndims == 3) {
p->nlev = 1;
p->size = p->nlat*p->nlon*p->nlev;
}
}
printf("Nlev = %ld size = %ld\n",p->nlev,p->size);
}
*/
void GetXtra(struct AR *A, char *infile) {

int ncid = -1, i;
int status = -1;
char vname[NC_MAX_NAME+1];
int nvarsp = -1;

memset(vname,'\0',sizeof(char)*(NC_MAX_NAME+1));

status = nc_open(infile, NC_NOWRITE, &ncid);
if (status != NC_NOERR) handle_error(status);

status = nc_inq(ncid,0,&nvarsp,0,0);
for(i = 0; i < nvarsp; i++) {
status = nc_inq_var(ncid,i,vname,0,0,0,0);
//printf("INFO: vname = %s\n",vname);
if(strcmp(vname,"lat") == 0) {
status = nc_get_var_double(ncid,i,&A->lat[0]);
if(status != NC_NOERR) handle_error(status);
}
if(strcmp(vname,"lon") == 0) {
status = nc_get_var_double(ncid,i,&A->lon[0]);
if(status != NC_NOERR) handle_error(status);
}
if(strcmp(vname,"hyai") == 0) {
status = nc_get_var_double(ncid,i,&A->a[0]);
if(status != NC_NOERR) handle_error(status);
}
if(strcmp(vname,"hybi") == 0) {
status = nc_get_var_double(ncid,i,&A->b[0]);
if(status != NC_NOERR) handle_error(status);
}
}
status = nc_close(ncid);
if(status != NC_NOERR) handle_error(status);
}

// print files in current directory in reverse order
int read_dir(char *path) {

int n = -1;
/* Read directory and find desired files */
n = scandir(path, &namelist, file_select, alphasort);
if (n < 1) {
perror(path); // This returns the system error to C
printf("ERROR: Number of files found = %d\n",n);
if (namelist) free(namelist);
nrerror("ERROR: Either the directory is empty and n = 0, Or there
is an error in the search pattern and n < 0");
}
return n;
}

/* Function that filters the input files */
//int file_select(const struct dirent *entry) {
int file_select(const struct dirent *entry) {
//printf("INFO: entry->d_name = %s\n",entry->d_name);
if ((strstr(entry->d_name,FILESTR) != NULL) &&
(strstr(entry->d_name,filter1) != NULL) &&
(strstr(entry->d_name,filter2) != NULL) &&
(strstr(entry->d_name,filter3) != NULL)) {
printf("INFO: entry->d_name = %s\n",entry->d_name);
return 1; // files found
}
else {
//printf("INFO: entry->d_name = %s\n",entry->d_name);
return 0; // files excluded
}
}

/*
Program usage
*/
void USAGE(int nargs,char **args) {

int i = -1;

fprintf(stderr,"ERROR: 4 arguments: a full path file(s), experiment,
year and if initial year\n");
if (nargs > 0) {
fprintf(stderr,"ERROR: Number of arguments -> %d\n",nargs);
for(i = 0; i < nargs; i++) {
fprintf(stderr,"ERROR: ARGV[%d] -> %s\n",i,args[i]);
}
}
nrerror("Incorrect number of arguments given");
}

/*
Handle errors by printing an error message and exiting with
a
non-zero status.
*/
void handle_error(int status) {
if(status != NC_NOERR) {
fprintf(stderr,"Error: %s\n",nc_strerror(status));
exit(EXIT_FAILURE);
}
}

float GeoidRadius(float lat) {

//[Re] = geoid radius as a function of latitude calculates the
radius of the geoid (m) at the
// given latitude (degrees). If lat is given in radians then
lat[i]*DEG2RAD is just lat[i]

float Rmax = EQRAD;
float Rmin = Rmax*(1.0-FLAT);
float D2R = PI/180;

return sqrt(1/(pow(cos(lat*D2R)/Rmax,2) + pow(sin(lat*D2R)/Rmin,
2)));
}

/*
Calculates the half pressures using the formula: P_level =
A_half_level + B_half_level * P_sfc
You must first calculate the half level before you can calculate the
full pressure levels using:
Pk = 1/2(Pk_+1/2 + Pk_1/2).
*/
void FullPresGph(struct V *T,
struct V *Q,
struct V *Z,
struct V *SP,
struct AR *A) {

float alpha = 0, g = 0, Re = 0, G = 0;
int i = 0, k = 0, lev = 0, lev2 = 0, j = 0, size = T->gridsize;
float D2R = PI/180;
float fpmax = 0.0, fpmin = 0.0, gphmax = 0.0, gphmin = 0.0;
float gm_max = 0.0, gm_min = 0.0;
float lvl_fpmax = 0.0, lvl_fpmin = 0.0;
float lvl_gphmax = 0.0, lvl_gphmin = 0.0;
float lvlgm_max = 0.0, lvlgm_min = 0.0;
float logph = 0;
float dp = 0, Tv = 0;

for (i = 0; i < T->nlat; i++) {
for(j = 0; j < T->nlon; j++) {
A->lats[(i*T->nlon)+j] = A->lat[i];
}
}
printf("INFO: Calculating the half pressure PH and half GPH\n");
for (k = T->ilev-1; k > -1; k--) {// Integrate from the ground up.
Loop over levels
lev = k*size;
for(i = 0; i < size; i++) {
// A and B values at pos 0 are 0. Ergo A/B levels length 0-90
A->ph[lev+i] = A->a[k] + A->b[k]*SP->mean[i]; // first half
level is pos 91 = SFC pressure: A= 0 & B = 1
if(k == T->ilev-1) A->hgph[lev+i] = Z->mean[i]/G0;
}
}
printf("INFO: Calculating the Full pressure, and GPH\n");

for (k = T->nlev-1; k > -1; k--) { // Integrate from the ground up.
Loop over levels
lev = k*size;
lev2 = (k+1)*size;
for(i = 0; i < size; i++) { // loop i
Tv = T->mean[lev+i]*(1 + (RV/RD -1)*Q->mean[lev+i]);
// The GPH at half levels (hgph) includes the surface elevation
at k == 90
if(k == 0) { // At the top of the model level A & B are both 0
alpha = log(2);
A->FP[lev+i] = (A->ph[lev+i] + A->ph[lev2+i])*0.5; // Full pressure
level
// integrate from previous (lower) half-level Fi to the full level.
A->gph[lev+i] = A->hgph[lev2+i] + (RD*Tv*alpha); // Convert to
decimeters
} else {
logph = log(A->ph[lev2+i]/A->ph[lev+i]);
dp = A->ph[lev2+i] - A->ph[lev+i];
alpha = 1.0 - ((A->ph[lev+i]/dp)*logph);
A->FP[lev+i] = (A->ph[lev+i] + A->ph[lev2+i])*0.5; // Full pressure
level
// integrate from previous (lower) half-level hgph to the full level.
A->gph[lev+i] = A->hgph[lev2+i] + (RD*Tv*alpha); // Convert to
decimeters
}
A->hgph[lev+i] = A->hgph[lev2+i] + RD*Tv*log(A->ph[lev2+i]/A-
>ph[lev+i]);
// Calculate the geometric height
g = 9.80616*(1-(0.002637*cos(2*A->lats[i]*D2R)) +
0.0000059*pow(cos(2*A->lats[i]*D2R),2));
G = g/G0;
Re = GeoidRadius(A->lats[i]);
A->GMH[lev+i] = (A->gph[lev+i]*Re/(G*Re - A->gph[lev+i]))/G0;
if(i == 0 && k == T->nlev-1) {
lvl_gphmax = A->gph[lev+i];
lvl_gphmin = A->gph[lev+i];
lvl_fpmax = A->FP[lev+i];
lvl_fpmin = A->FP[lev+i];
lvlgm_max = A->GMH[lev+i];
lvlgm_min = A->GMH[lev+i];
} else {
if(A->GMH[lev+i] > lvlgm_max) lvlgm_max = A->GMH[lev+i];
if(A->GMH[lev+i] < lvlgm_min) lvlgm_min = A->GMH[lev+i];
if(A->gph[lev+i] > lvl_gphmax) lvl_gphmax = A->gph[lev+i];
if(A->gph[lev+i] < lvl_gphmin) lvl_gphmin = A->gph[lev+i];
if(A->FP[lev+i] > lvl_fpmax) lvl_fpmax = A->FP[lev+i];
if(A->FP[lev+i] < lvl_fpmin) lvl_fpmin = A->FP[lev+i];
}
} // end for loop for i
if(k == T->nlev-1) {
fpmax = lvl_fpmax;
fpmin = lvl_fpmin;
gphmax = lvl_gphmax;
gphmin = lvl_gphmin;
gm_max = lvlgm_max;
gm_min = lvlgm_min;
} else {
if(lvlgm_max > gm_max) gm_max = lvlgm_max;
if(lvlgm_min < gm_min) gm_min = lvlgm_min;
if(lvl_gphmax > gphmax) gphmax = lvl_gphmax;
if(lvl_gphmin < gphmin) gphmin = lvl_gphmin;
if(lvl_fpmax > fpmax) fpmax = lvl_fpmax;
if(lvl_fpmin < fpmin) fpmin = lvl_fpmin;
}
} // end loop for k

printf("INFO: Completed GMH and FP calculations.\n");

return;
}

void RHRho(struct V *Q,
struct V *T,
struct V *SS,
struct V *CS,
struct V *CI,
struct AR *A) {

float RS = 8.314;
float Rho_d = 0.0;
float Rho_v = 0.0;
float Rho = 0.0;
float Rho2 = 0.0;
float e = 0.0;
int k;

for (k = 0; k < T->size; k++) {
// Density of h2o
Rho_v = (A->FP[k]/(RS*T->mean[k])) * 1/( (1 - Q->mean[k])/(Q-
>mean[k]*Md*1e-3) + 1/Mv*1e-3);
// Density for dry air
Rho_d = Rho_v*(1-Q->mean[k])/Q->mean[k];
// Total density
Rho = Rho_d*(1 + Q->mean[k])/(1 + 1.609*Q->mean[k]);
Rho2 = Rho_v + Rho_d;
e = Rho_v*RS/Mv*T->mean[k]; // Partial pressure
//EI(T[k]); // saturation pressure
A->RH[k] = e/(EI(T->mean[k]))*100; // relative humidty

// Convert to g/m2
CS->mean[k] = CS->mean[k]*Rho*1000;
SS->mean[k] = SS->mean[k]*Rho*1000;
CI->mean[k] = CI->mean[k]*Rho*1000;
}
return;
}

void IW_PATH(struct V *CS,
struct V *SS,
struct V *CI,
struct AR *A) {
int k, i;
int size = CS->gridsize;

for(i = 0; i < CS->gridsize; i++) A->CIWP[i] = 0.0;
for(i = 0; i < CS->gridsize; i++) A->TCIWP[i] = 0.0;

for (k = CS->nlev-2; k > -1; k--) {
int lev = k*size;
int lev2 = (k+1)*size;
for(i = 0; i < size; i++) {
A->CIWP[i] += (A->GMH[lev+i]-A->GMH[lev2+i])*((CI->mean[lev+i]
+CI->mean[lev2+i])/2.0);
A->TCIWP[i] += (A->GMH[lev+i]-A->GMH[lev2+i])*(((CS->mean[lev+i]
+CI->mean[lev+i]+SS->mean[lev+i])+(SS->mean[lev2+i]+CS->mean[lev2+i]
+CI->mean[lev2+i]))/2.0);
}
}
printf("INFO: Printing the max, min and mean of CIWP and TCIWP\n");
MMF(A->CIWP,CS->gridsize);
MMF(A->TCIWP,CS->gridsize);
MEAN(A->CIWP,CS->gridsize);
MEAN(A->TCIWP,CS->gridsize);

return;
}

float EW(float T) {

// Coefficients for Ew:
float a = -6096.9385;
float b = 21.2409642;
float c = -2.711193e-2;
float d = 1.673952e-5;
float e = 2.433502;

return exp(a/T + b + c*T + d*pow(T,2) + e*log(T));
}

float EI(float T) {

// Coefficients for Ei:
float a = -6024.5282;
float b = 29.32707;
float c = 1.0613868e-2;
float d = -1.3198825e-5;
float e = -0.49382577;

return exp(a/T + b + c*T + d*pow(T,2) + e*log(T));
}


void A1DTo2D(struct V *T, float *a1d, float **a2d) {

int lat,lon,count = 0;

for(lat = 0; lat < T->nlat; lat++) {
for (lon = 0; lon < T->nlon; lon++) {
a2d[lat][lon] = (float)a1d[count++];
}
}
return;
}

void A1DTo3D(struct V *T, float *a1d, float ***a3d) {

int lvl, lat, lon, count = 0;

for (lvl = 0; lvl < T->nlev; lvl++) {
for (lat = 0; lat < T->nlat; lat++) {
for (lon = 0; lon < T->nlon; lon++) {
a3d[lvl][lat][lon] = (float)a1d[count++];
}
}
}
return;
}

void GetMean(struct AR *A, struct V *ret, int timestep, double *array)
{
int n;
int tsID = 0, vdate = 0, vtime = 0;
int nmiss = 0, ts = 0, cdiID = -1;
int streamID, taxisID, vlistID;
char name[NC_MAX_NAME+1];

memset(name,'\0',sizeof(char)*(NC_MAX_NAME+1));
for(n = 0; n < ret->size; n++) ret->mean[n] = 0.0;

streamID = streamOpenRead(ret->infile);
if(streamID < 0) {
fprintf(stderr,"ERROR: %s\n", cdiStringError(streamID));
exit(EXIT_FAILURE);
}
vlistID = streamInqVlist(streamID);
taxisID = vlistInqTaxis(vlistID);
for(n = 0; n < vlistNvars(vlistID); n++) {
vlistInqVarName(vlistID,n,name);
//printf("CDI %d -> name: %s\n",n,name);
if(strcmp(name,ret->vname) == 0) {
printf("INFO: Assigned CDI %d -> name: %s\n",n,name);
cdiID = n;
}
}
while((streamInqTimestep(streamID,tsID++)) != 0) {
vdate = taxisInqVdate(taxisID)/100;
vtime = taxisInqVtime(taxisID);
//printf("INFO: TS[%03d]-> VDATE = %d VTIME = %d
\n",tsID,vdate,vtime);
if(vdate == atoi(ret->date)) {
if(vtime == timestep) {
printf("INFO: TS[%03d]-> VDATE = %d VTIME = %d date = %s TS = %d
\n",tsID,vdate,vtime,ret->date,timestep);
streamReadVar(streamID,cdiID,array,&nmiss);
if(strcmp(ret->vname,"Z") == 0) {
for(n = 0; n < ret->size; n++) ret->mean[n] = (float)array[n];
return;
}
printf("INFO: Adding the values to the mean array\n");
for(n = 0; n < ret->size; n++) ret->mean[n] += (float)array[n];
ts++;
printf("INFO: Extracted time step: %d\n",ts);
break;
}
}
} // End while loop over time steps
// Extract the time step from the previous date
streamClose(streamID);

tsID = 0;

streamID = streamOpenRead(ret->pinfile);
if(streamID < 0) {
fprintf(stderr,"ERROR: %s\n", cdiStringError(streamID));
exit(EXIT_FAILURE);
}
vlistID = streamInqVlist(streamID);
taxisID = vlistInqTaxis(vlistID);
//printf("Number of variables in file %d\n",vlistNvars(vlistID));
for(n = 0; n < vlistNvars(vlistID); n++) {
vlistInqVarName(vlistID,n,name);
//printf("CDI %d -> name: %s\n",n,name);
if(strcmp(name,ret->vname) == 0) {
printf("INFO: Assigned CDI %d -> name: %s\n",n,name);
cdiID = n;
}
}
while((streamInqTimestep(streamID,tsID++)) != 0) {
vdate = taxisInqVdate(taxisID)/100;
vtime = taxisInqVtime(taxisID);
// Read variable
//printf("INFO: TS[%03d]-> VDATE = %d VTIME = %d
\n",tsID,vdate,vtime);
if(vdate == atoi(ret->date)) {
if(vtime == timestep) {
printf("INFO: TS[%03d]-> VDATE = %d VTIME = %d date = %s TS = %d
\n",tsID,vdate,vtime,ret->date,timestep);
memset(array,0.0,sizeof(array));
streamReadVar(streamID,cdiID,array,&nmiss);
Print("INFO: Adding the values to the mean array");
for(n = 0; n < ret->size; n++) ret->mean[n] += (float)array[n];
ts++;
printf("INFO: Extracted time step: %d\n",ts);
}
}
} // End while loop over time steps
printf("INFO: Returning the mean of %d days!\n",ts);

if(strcmp(ret->vname,"TTR") == 0) {
if(timestep == 0) {
for(n = 0; n < ret->size; n++) ret->mean[n] = ret->mean[n]/
(ts*3600);
} else {
for(n = 0; n < ret->size; n++) ret->mean[n] = ret->mean[n]/
(ts*timestep);
}
} else {
for(n = 0; n < ret->size; n++) ret->mean[n] = ret->mean[n]/ts;
}
streamClose(streamID);
return;
}
*****

Some extra functions:

**********
#include "utilities.h"

void MMD(double *array, size_t size) {

int i;
double min = array[0], max = array[0];

for(i = 0; i < size; i++) {
if(array[i] < min) min = array[i];
if(array[i] > max) max = array[i];
}
printf("INFO: Min = %12.6lf\tMax = %12.6lf\n",min,max);
}

void MMF(float *array, size_t size) {

int i;
float min = array[0], max = array[0];

for(i = 0; i < size; i++) {
if(array[i] < min) min = array[i];
if(array[i] > max) max = array[i];
}
printf("INFO: Min = %12.6f\tMax = %12.6f\n",min,max);
}

float MaxF(double *array, size_t size) {

int i;
float max = array[0];

for(i = 0; i < size; i++) if(array[i] > max) max = array[i];
return max;
}

float MinF(double *array, size_t size) {

int i;
float min = array[0];

for(i = 0; i < size; i++) if(array[i] < min) min = array[i];

return min;
}


void MMI(int *array, size_t size) {

int i;
int min = array[0], max = array[0];

for(i = 0; i < size; i++) {
if(array[i] < min) min = array[i];
if(array[i] > max) max = array[i];
}
printf("INFO: Min = %-10d\tMax = %-10d\n",min,max);
}

void MinMaxD(double **array, size_t nrows, size_t ncols) {

int i,j;
double min = array[0][0], max = array[0][0];

for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
if(array[i][j] < min) min = array[i][j];
if(array[i][j] > max) max = array[i][j];
}
}
printf("INFO: Min = %12.6lf\tMax = %12.6lf\n",min,max);
}

void MinMaxF(float **array, size_t nrows, size_t ncols) {

int i,j;
float min = array[0][0], max = array[0][0];

for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
if(array[i][j] < min) min = array[i][j];
if(array[i][j] > max) max = array[i][j];
}
}
printf("INFO: Min = %-12.6lf\tMax = %-12.6lf\n",min,max);
}

void nrerror(const char *error_text) {
/* standard error handler */
fprintf(stderr,"ERROR: %s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(EXIT_FAILURE);
}

void Print(const char *verbose_text) {
/* standard info handler */
printf("INFO: %s\n",verbose_text);
}

void PrintI(int num) {

printf("INFO: Integer number is %d\n",num);

}

void PrintF(float num) {

printf("INFO: Float number is %12.6f\n",num);

}

void PrintL(long num) {

printf("INFO: Long integer number is %ld\n",num);

}

void PrintD(double num) {

printf("INFO: Double number is %-12.6lf\n",num);

}

char **cmatrix(size_t nrows, size_t ncols) {

/* allocate an unsigned char matrix with subscript range
m[0..nrows-1][0..ncols-1] */
int i;
char **m;

/* allocate pointers to rows */
m = (char **)smalloc(nrows*sizeof(char *));
if (!m) nrerror("allocation failure 1 in imatrix()");
/* allocate rows and set pointers to them */
for(i = 0; i < nrows; i++) {
m[i] = (char *)smalloc(ncols*sizeof(char));
if (!m[i]) nrerror("allocation failure 2 in imatrix()");
}

/* return pointer to array of pointers to rows */
return m;
}

void FreeCon2DF(float **array) {
free((void *)array[0]);
free((void *)array);
//Print("Array deallocated!");
}

void FreeCon2DD(double **array) {
free((void *)array[0]);
free((void *)array);
//Print("Array deallocated!");
}

void FreeCon3DF(float ***array) {

free((void *)array[0]);
free((void *)array);
}

void FreeCon3DD(double ***array) {

free((void *)array[0]);
free((void *)array);
}

void FreeCon3DSI(short int ***array) {

free((void *)array[0]);
free((void *)array);
}


/*
In this method we allocate a block of memory to hold the whole array
first.
We then create an array of pointers to point to each row. Thus even
though
the array of pointers is being used, the actual array in memory is
contiguous.
The code looks like this:
*/
short int ***ConMem3DSI(size_t nlevs, size_t nrows, size_t ncols) {

short int *space;
short int ***Arr3D;
int y, z, i, j, k;

/* first we set aside space for the array itself */
space = smalloc(nlevs*nrows*ncols*sizeof(*space));

/* next we allocate space of an array of pointers, each
to eventually point to the first element of a
2 dimensional array of pointers to pointers
*/
Arr3D = smalloc(nlevs*sizeof(**Arr3D));

/* and for each of these we assign a pointer to a newly
allocated array of pointers to a row
*/
for (z = 0; z < nlevs; z++) {
Arr3D[z] = smalloc(nrows*sizeof(*Arr3D));

/* and for each space in this array we put a pointer to
the first element of each row in the array space
originally allocated
*/
for (y = 0; y < nrows; y++) {
Arr3D[z][y] = space + (z*(nrows*ncols) + y*ncols);
}
}
for(k = 0; k < nlevs; k++) {
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
Arr3D[k][i][j] = 0.0;
}
}
}
return Arr3D;
}

float ***ConMem3DF(size_t nlevs, size_t nrows, size_t ncols) {

float *space;
float ***Arr3D;
int y, z, i, j, k;

/* first we set aside space for the array itself */
space = smalloc(nlevs*nrows*ncols*sizeof(*space));

/* next we allocate space of an array of pointers, each
to eventually point to the first element of a
2 dimensional array of pointers to pointers
*/

Arr3D = smalloc(nlevs*sizeof(**Arr3D));

/* and for each of these we assign a pointer to a newly
allocated array of pointers to a row
*/
for (z = 0; z < nlevs; z++) {
Arr3D[z] = smalloc(nrows*sizeof(*Arr3D));

/* and for each space in this array we put a pointer to
the first element of each row in the array space
originally allocated
*/
for (y = 0; y < nrows; y++) {
Arr3D[z][y] = space + (z*(nrows*ncols) + y*ncols);
}
}
for(k = 0; k < nlevs; k++) {
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
Arr3D[k][i][j] = 0.0;
}
}
}
return Arr3D;
}

double ***ConMem3DD(size_t nlevs, size_t nrows, size_t ncols) {

double *space;
double ***Arr3D;
int y, z, k, i, j;

/* first we set aside space for the array itself */
space = smalloc(nlevs*nrows*ncols*sizeof(*space));

/* next we allocate space of an array of pointers, each
to eventually point to the first element of a
2 dimensional array of pointers to pointers
*/

Arr3D = smalloc(nlevs*sizeof(**Arr3D));

/* and for each of these we assign a pointer to a newly
allocated array of pointers to a row
*/
for (z = 0; z < nlevs; z++) {
Arr3D[z] = smalloc(nrows*sizeof(*Arr3D));

/* and for each space in this array we put a pointer to
the first element of each row in the array space
originally allocated
*/
for (y = 0; y < nrows; y++) {
Arr3D[z][y] = space + (z*(nrows*ncols) + y*ncols);
}
}
for(k = 0; k < nlevs; k++) {
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
Arr3D[k][i][j] = 0.0;
}
}
}
return Arr3D;
}

float **ConMem2DF(size_t nrows, size_t ncols) {

int i,j;

float **array = smalloc(nrows * sizeof(*array));
array[0] = smalloc(nrows * ncols*sizeof(**array));
if(array[0] == NULL) nrerror("Failed to allocate memory!");
for(i = 0; i < nrows; i++) {
array[i] = array[0] + i * ncols;
}
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
array[i][j] = 0;
}
}
return array;
}

double **ConMem2DD(size_t nrows, size_t ncols) {

int i,j;

double **array = smalloc(nrows * sizeof(*array));
array[0] = smalloc(nrows * ncols*sizeof(**array));
if(array[0] == NULL) nrerror("Failed to allocate memory!");
for(i = 0; i < nrows; i++) {
array[i] = array[0] + i * ncols;
}
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
array[i][j] = 0;
}
}
return array;
}

/*
Calculates the standard deviation
*/
float SD_D(double **array, size_t nrows, size_t ncols) {

double average = 0;
double sum1 = 0, sum2 = 0;
int i, j;

for( i = 0; i < nrows; i++) {
for( j = 0; j < ncols; j++) {
sum1 += array[i][j];
}
}
average = sum1/(nrows*ncols);

for( i = 0; i < nrows; i++) {
for( j = 0; j < ncols; j++) {
sum2 += ((array[i][j] - average) * (array[i][j] - average));
}
}
return (float)sqrt((1.0/(nrows*ncols - 1)) * sum2);
}

float SD_F(float ***array, int k, int nrows, int ncols) {

double average = 0;
double sum1 = 0, sum2 = 0;
int i, j;

for( i = 0; i < nrows; i++) {
for( j = 0; j < ncols; j++) {
sum1 += array[k][i][j];
}
}
average = sum1/(nrows*ncols);

for( i = 0; i < nrows; i++) {
for( j = 0; j < ncols; j++) {
sum2 += ((array[k][i][j] - average) * (array[k][i][j] -
average));
}
}
return (float)sqrt((1.0/(nrows*ncols - 1)) * sum2);
}

/*
Returns the min value
*/
void MMValueI(int ***array, int nlevs, int nrows, int ncols) {

int i,j,k;
int min = array[0][0][0], max = array[0][0][0];

for(k = 0; k < nlevs; k++) {
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
if(array[k][i][j] < min) min = array[k][i][j];
if(array[k][i][j] > max) max = array[k][i][j];
}
}
}
printf("INFO: Min = %-10d\tMax = %-10d\n",min,max);
return;
}

void MMValueF(float ***array, int nlevs, int nrows, int ncols) {

int i,j,k;
float min = array[0][0][0], max = array[0][0][0];


for(k = 0; k < nlevs; k++) {
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
if(array[k][i][j] < min) min = array[k][i][j];
if(array[k][i][j] > max) max = array[k][i][j];
}
}
}
printf("INFO: Min = %12.4f\tMax = %12.4f\n",min,max);
return;
}

/*
Returns the max vaule
*/
void MMValueD(double ***array, int nlevs, int nrows, int ncols) {

int i,j,k;
double min = array[0][0][0], max = array[0][0][0];

for(k = 0; k < nlevs; k++) {
for(i = 0; i < nrows; i++) {
for(j = 0; j < ncols; j++) {
if(array[k][i][j] < min) min = array[k][i][j];
if(array[k][i][j] > max) max = array[k][i][j];
}
}
}
printf("INFO: Min = %12.6lf\tMax = %12.6lf\n",min,max);
return;
}

void MEAN(float *array, size_t size) {

int k;
float sum = 0.0;

for(k = 0; k < size; k++) sum += array[k];

printf("INFO: Mean = %12.6lf\n",sum/size);
return;
}

void *smalloc(size_t size) {

void *mem = malloc(size);
if(!(mem)) nrerror("Not enough memory!");
return mem;
}

**********

Eric Sosman

unread,
Mar 14, 2010, 4:18:25 PM3/14/10
to
On 3/14/2010 3:46 PM, Sheldon wrote:
> On 14 mar, 19:53, Eric Sosman<esos...@ieee-dot-org.invalid> wrote:
>> On 3/14/2010 1:57 PM, Sheldon wrote:
>>
>>> Hi everyone,
>>
>>> I have a fairly large program that has a problem that I cannot seem to
>>> solve.
>>> [...]

It's one of the possibilities I suggested you look for:
calling free() on memory not obtained from malloc() and friends.
In the case at hand, you're obtaining memory from something called
smalloc(), which is not part of the Standard C library. We're
not sure exactly what it is, but there's at least one package of
that name floating around the Net -- and that package, according
to its description, allocates from static memory, not from dynamic
memory. So when you pass it to free() (or to Free() -- Lord only
knows that *that* is), you have sinned and you get punished.

Either find and use the smalloc() package's "moral equivalent"
to free() -- maybe sfree()? -- or switch to using true malloc()
throughout.

I also notice that you often fail to check whether an
allocation request succeeded or failed. Yes, I know you *said*
you checked, but in the actual code you frequently don't. This
slipshod practice may eventually bite you in the behind -- and
perhaps it already has.

Go, and sin no more.

--
Eric Sosman
eso...@ieee-dot-org.invalid

Phil Carmody

unread,
Mar 14, 2010, 5:19:27 PM3/14/10
to
Sheldon <shej...@gmail.com> writes:
> Hi everyone,
>
> I have a fairly large program that has a problem that I cannot seem to
> solve.
> The program runs smoothly but when it tries to free allocated memory
> before exiting I get a segmentation fault.
> I have checked many times the following:
>
> 1. That I have not freed the memory before hand and call free this
> way: if(array) free(array);

I know you say you've checked many times, but if you want to
decorate your use of free(), then do
free(array); array=NULL;
instead. This might help your error show itself more quickly.
But setting array to something like 0x00000001 might make it
show itself even quicker still, as a double-free will almost
certainly make itself known immediately. This is called 'poison',
and is a useful technique for debugging.

> The code is nearly 1000 lines long. I can post it if anyone wants to
> look at it or I can send it via email.
> I would greatly appreciate any help with this problem.

Quicker would be for you to link with something like 'electric fence',
and then run it within a debugging environment.

Phil
--
I find the easiest thing to do is to k/f myself and just troll away
-- David Melville on r.a.s.f1

Sheldon

unread,
Mar 14, 2010, 5:54:43 PM3/14/10
to

Hi Eric,

smalloc is just a function that I use. Nothing from the net but an
idea I got from a friend:

void *smalloc(size_t size) {
void *mem = malloc(size);
if(!(mem)) nrerror("Not enough memory!");
return mem;
}

As you can see the error checking is embedded. And I'm sure it works
since all the memory was allocated and the program actually does very
detailed calculations without any errors. You are right though, it's
one of the cases you mentioned - I think. I know that I loop along the
arrays many times in the program. So I'm currently checking this. I've
used Valgrind and I run the program in GDB but neither showed me the
error.

I don't know how to use electric fence but I guess I better learn.

The other "Free..." is something I got from this site when I wanted to
how to create 2D and 3D arrays with contiguous memory. All of these
functions I've included above - I know, it's a lot to sift through so
here they are:

void FreeCon2DF(float **array) {
free((void *)array[0]);
free((void *)array);
//Print("Array deallocated!");
}

void FreeCon2DD(double **array) {
free((void *)array[0]);
free((void *)array);
//Print("Array deallocated!");
}

void FreeCon3DF(float ***array) {
free((void *)array[0]);
free((void *)array);
}

void FreeCon3DD(double ***array) {
free((void *)array[0]);
free((void *)array);
}

void FreeCon3DSI(short int ***array) {
free((void *)array[0]);
free((void *)array);
}

I'm beginning to think this is a better way to free memory whether
it's 1D, or greater. It's late here so I'll test these things tomorrow
morning and include your good suggestions.

/M

Ian Collins

unread,
Mar 14, 2010, 6:40:56 PM3/14/10
to
On 03/15/10 10:54 AM, Sheldon wrote:

Did you really have to quote all that again?

> smalloc is just a function that I use. Nothing from the net but an
> idea I got from a friend:

<snip>

> The other "Free..." is something I got from this site when I wanted to
> how to create 2D and 3D arrays with contiguous memory. All of these
> functions I've included above - I know, it's a lot to sift through so
> here they are:

<snip>

> void FreeCon3DF(float ***array) {
> free((void *)array[0]);
> free((void *)array);
> }
>
> void FreeCon3DD(double ***array) {
> free((void *)array[0]);
> free((void *)array);
> }
>
> void FreeCon3DSI(short int ***array) {
> free((void *)array[0]);
> free((void *)array);
> }

These look well dodgy. The corresponding allocation functions call
smalloc multiple times, so there will be more allocations than frees.

--
Ian Collins

pete

unread,
Mar 14, 2010, 10:33:12 PM3/14/10
to
Eric Sosman wrote:
> On 3/14/2010 1:57 PM, Sheldon wrote:
>
>> Hi everyone,
>>
>> I have a fairly large program that has a problem that I cannot seem to
>> solve.
>> The program runs smoothly but when it tries to free allocated memory
>> before exiting I get a segmentation fault.
>> I have checked many times the following:
>>
>> 1. That I have not freed the memory before hand and call free this
>> way: if(array) free(array);
>
>
> A note in passing: This makes no difference, because
> free(NULL) is legal (and does nothing).

And it also makes no difference because

if(freed_pointer) free(freed_pointer);

is still undefined.


--
pete

Nick Keighley

unread,
Mar 15, 2010, 4:18:36 AM3/15/10
to
On 14 Mar, 17:57, Sheldon <shejo...@gmail.com> wrote:

> I have a fairly large program that has a problem that I cannot seem to
> solve.
> The program runs smoothly but when it tries to free allocated memory
> before exiting I get a segmentation fault.
> I have checked many times the following:
>
> 1. That I have not freed the memory before hand and call free this
> way: if(array) free(array);

this has no affect if array is not set to NULL (or someother bad
value). You are aware that free() doesn't set array to NULL?

This is erroneous code:-

if (array) free (array);
if (array) free (array);

ass array gets freed twice.

I really don't like your FreeConxxx() functions and similar. Why the
cast to void*? Why do you think freeing array[0] id going to do
anything magic? What about the rest of the array?
You could try adding some logging (this is assuming you aren't just
going to run electric fence or whatever) to you program. Log every
malloc and free call to a file. I suspect the mallocs and fgrees won't
match.

malloc 100 bytes at BFEE0010
malloc 200 bytes at BFEE0114
free at BFEE0010

I bet they won't all match up correctly.

BTW: Windows has stuff to help debug these sorts of problems. Not as
good as valgrind etc. but readily available.

Nick Keighley

unread,
Mar 15, 2010, 5:30:36 AM3/15/10
to

why?

Sheldon

unread,
Mar 15, 2010, 6:49:30 AM3/15/10
to
On 15 mar, 09:18, Nick Keighley <nick_keighley_nos...@hotmail.com>
wrote:

Thanks guys for your help. The problem is now solved. I switched from
malloc to calloc and that did it.
I also added some of the suggested improvements but still the error
remained.
Any ideas why calloc worked and malloc didn't?

Thanks for your help and ideas/critiques :-)

/S

Ben Bacarisse

unread,
Mar 15, 2010, 9:07:42 AM3/15/10
to
Sheldon <shej...@gmail.com> writes:
<snip>

> Thanks guys for your help. The problem is now solved. I switched from
> malloc to calloc and that did it.
> I also added some of the suggested improvements but still the error
> remained.
> Any ideas why calloc worked and malloc didn't?

I'd be worried! You may have fixed it, but in my experience that is
very rare. Far more often what you will have done is to mask the
error so that there is no longer any visible sign of it. I'd switch
right back because I don't like to loose track of errors before I've
fixed them.

BTW, since you like to zero out your data, I would switch back to
calloc for the "main" part of the arrays and remove the zeroing loops
once you find the bug.

It's hard for people here to do more than inspect the code because
it is not complete and there is no data to test with. I am not
certain, but might get more help if you posted in a form that was
complete along with a test data set. I'd use a website for this since
the program is just large enough that cutting out the various files
from a linear posting is tedious.

> Thanks for your help and ideas/critiques :-)

Just a few thoughts...

Someone else has brought this up, but you certainly have a memory
leak. The FreeCon3Dxx functions don't free up the patterns of
allocation made by the Con3Dxx functions. The 2D ones look OK.

There is much duplication of code. I'd consider switching to C++ just
so you can write all you allocation and stats functions as templates.
This will save you maintaining three versions of each one. This is a
simple use of C++ and there is no need to learn the whole language to
do this.

If you want to stick with C, you could consider using macros to make
the various differently types versions of your utility functions. All
lo priority since the work is done. It's just a tidy up.

Because of the way you allocate 1, 2 and 3D allays, you could write
the various min and max function all in terms of a single helper
function (at least for each base type) that treats the array a 1D.

You will get some simplification if you believe (as you should) that
sizeof(char) is 1. If you worry that you might change the char type
to one that is wider (wchar_t) then the usual:

p = malloc(N * sizeof *p);

idiom that you use elsewhere will do the trick.

--
Ben.

pete

unread,
Mar 15, 2010, 9:09:15 AM3/15/10
to

I intended for (freed_pointer) to represent a previously freed pointer.

--
pete

Nick Keighley

unread,
Mar 15, 2010, 10:13:20 AM3/15/10
to


ah, sorry I get you now. I was treating at just an arbitary name

if (zczc_floop) free(zczc_floop);

Keith Thompson

unread,
Mar 15, 2010, 12:54:36 PM3/15/10
to
Nick Keighley <nick_keigh...@hotmail.com> writes:
[...]

> This is erroneous code:-
>
> if (array) free (array);
> if (array) free (array);
>
> ass array gets freed twice.
[...]

In practice, yes, that will probably result in a second call to free()
on the same pointer, which invokes undefined behavior.

In principle, it's the second test ("if (array)") that invokes
undefined behavior. Once a pointer value has been passed to free(),
that pointer value becomes indeterminate, and any reference to it is
undefined.

(Which means that, given some smarts from the compiler, the value of
"array" *could* be set to NULL after the call to free() -- though
I don't think the compiler would be doing you any favors if it
performed this kind of magic.)

I'm assuming, of course, that "array" is the name of a pointer
variable. If so, it should probably have a better name.

--
Keith Thompson (The_Other_Keith) ks...@mib.org <http://www.ghoti.net/~kst>
Nokia
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

pete

unread,
Mar 15, 2010, 3:32:14 PM3/15/10
to


The fault was mine. I was overly terse.
I thank you for getting me to clarify what I wrote.

--
pete

Richard Bos

unread,
Mar 16, 2010, 6:13:11 AM3/16/10
to
Keith Thompson <ks...@mib.org> wrote:

> Nick Keighley <nick_keigh...@hotmail.com> writes:
> [...]
> > This is erroneous code:-
> >
> > if (array) free (array);
> > if (array) free (array);
> >
> > ass array gets freed twice.
> [...]
>
> In practice, yes, that will probably result in a second call to free()
> on the same pointer, which invokes undefined behavior.
>
> In principle, it's the second test ("if (array)") that invokes
> undefined behavior. Once a pointer value has been passed to free(),
> that pointer value becomes indeterminate, and any reference to it is
> undefined.
>
> (Which means that, given some smarts from the compiler, the value of
> "array" *could* be set to NULL after the call to free() -- though
> I don't think the compiler would be doing you any favors if it
> performed this kind of magic.)

Not within reason, it couldn't, since the bit pattern, as seen through
an unsigned char pointer, has to remain the same. You could possibly
argue that, after free()ing, that unaltered bit pattern could now
_compare_ equal to a null pointer when used as a pointer, but that is
not the same thing, IMO, as _setting_ it to NULL. Specifically, after

int *pointer;
unsigned memory_null[sizeof pointer], memory_ptr[sizeof pointer];

pointer=NULL;
memcpy(memory_null, &pointer, sizeof pointer);
pointer=malloc(42);
if (pointer==NULL)
puts("There's bugger all down here on Earth.");
else {
memcpy(memory_ptr, &pointer, sizeof pointer);
free(pointer);
if (memcmp(memory_null, &pointer, sizeof pointer)==0)
puts("Run away! Run away!");
if (memcmp(memory_ptr, &pointer, sizeof pointer)!=0)
puts("I soiled my armour...");
}

your implementation must _not_ have quoted from Monty Python and the
Holy Grail (though it may have quoted from The Meaning of Life).

Richard

Keith Thompson

unread,
Mar 16, 2010, 12:27:59 PM3/16/10
to
ral...@xs4all.nl (Richard Bos) writes:
> Keith Thompson <ks...@mib.org> wrote:
>> Nick Keighley <nick_keigh...@hotmail.com> writes:
>> [...]
>> > This is erroneous code:-
>> >
>> > if (array) free (array);
>> > if (array) free (array);
>> >
>> > ass array gets freed twice.
>> [...]
>>
>> In practice, yes, that will probably result in a second call to free()
>> on the same pointer, which invokes undefined behavior.
>>
>> In principle, it's the second test ("if (array)") that invokes
>> undefined behavior. Once a pointer value has been passed to free(),
>> that pointer value becomes indeterminate, and any reference to it is
>> undefined.
>>
>> (Which means that, given some smarts from the compiler, the value of
>> "array" *could* be set to NULL after the call to free() -- though
>> I don't think the compiler would be doing you any favors if it
>> performed this kind of magic.)
>
> Not within reason, it couldn't, since the bit pattern, as seen through
> an unsigned char pointer, has to remain the same.
[...]

Yes, good point, I should have remembered that.

Ok, if the compiler can prove that you don't happen to examine the
pointer's representation as an array of characters, then it can do
anything it likes with it. For that matter, if it can prove that you
do examine it as a pointer, it can do anything it likes with it, since
the program's behavior is undefined.

For example, this program:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int main(void)
{
int *p = malloc(sizeof *p);
int *null_pointer = NULL;
if (p == NULL) {
exit(EXIT_FAILURE);
}
free(p);
if (p == NULL) {
puts("Hmm, that's odd");
}
if (memcmp(&p, &null_pointer, sizeof p) == 0) {
puts("Hmm, that's really odd");
}
return 0;
}

can legally print one or both of "Hmm, that's odd" and "Hmm, that's
really odd". Of course, it can also legally print "Whiskey Tango
Foxtrot".

But I think I've gone rather far afield from the actual point, which
is that once you attempt to access the value of a free()d pointer,
all bets are off. In practice, non-DS9K compilers aren't likely to
go out of their way to make such a program behave in an unexpected
manner, but optimization can still do some very surprising things.

io_x

unread,
Mar 16, 2010, 1:45:26 PM3/16/10
to

"Richard Bos" <ral...@xs4all.nl> ha scritto nel messaggio
news:4b9ebd96...@news.xs4all.nl...

i have not understand enought but

> int *pointer;
> unsigned memory_null[sizeof pointer], memory_ptr[sizeof pointer];

here in my machine, above memory_null and memory_ptr are both 4 unsigned
for a total memory of 4*4+4*4 chars
but here sizeof pointer should be 4 chars you use more space than enought
the enought should be 4 + 4 chars

than who says that int* has the same aligned of unsigned ?
why not easily write just "int *pointer, *memory_null, *memory_ptr;" ?

all the remain seems ok

David Thompson

unread,
Mar 31, 2010, 2:55:44 AM3/31/10
to
On 14 Mar 2010 18:09:33 GMT, Seebs <usenet...@seebs.net> wrote:

> On 2010-03-14, Sheldon <shej...@gmail.com> wrote:
> > I have a fairly large program that has a problem that I cannot seem to
> > solve.
> > The program runs smoothly but when it tries to free allocated memory
> > before exiting I get a segmentation fault.

<snip>


> A thousand lines isn't all that much. A couple thoughts:
>
> 1. Go through removing functional code to strip it down. If you find
> that, after removing something, you suddenly don't have a problem, then
> that code is what did it.

Warning: If it's a memory clobber bug, removing some (good) code that
doesn't contain the bug may merely hide the symptom(s) because the
clobber moves to a different place that doesn't (currently) matter, or
disguise them because the clobber only causes slightly wrong output
rather than a crash. For example if it now damages your PRNG and
(only) causes the output to be less random than it should, you have
practically no chance of noticing that by looking at output. This can
mislead you into looking at the good code for a bug that isn't there.
You may need to test a whole lattice of cuts to find, with any
confidence, which ones actually remove the problem.

> 2. Once you get it stripped down, it should be a lot shorter.
>

That's true. If you are able to find a small(er) version that DOES
still have the symptom, you have a smaller haystack to deal with.

> My guess would be that you're overrunning something -- allocating the
> wrong size of space. The last time I had this, I'd allocated enough space
> for a pointer, rather than enough space for the thing pointed to...
>

As others have said, the easiest way to look for this is one of the
several tools that actually checks for overruns. I don't know which
(or if any) also catch wild or stale pointers, though.

0 new messages