nc.c

一个用来实现偏微分方程中网格的计算库

/*
 *	Copyright 1996, University Corporation for Atmospheric Research
 *      See netcdf/COPYRIGHT file for copying and redistribution conditions.
 */
/* $Id: nc.c 2501 2007-11-20 02:33:29Z benkirk $ */

#include "nc.h"
#include "rnd.h"
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "ncx.h"
#if defined(LOCKNUMREC) /* && _CRAYMPP */
#  include <mpp/shmem.h>
#  include <intrinsics.h>
#endif

/* list of open netcdf's */
static NC *NClist = NULL;

/* This is the default create format for nc_create and nc__create. */
int default_create_format = NC_FORMAT_CLASSIC;

/* These have to do with version numbers. */
#define MAGIC_NUM_LEN 4
#define VER_CLASSIC 1
#define VER_64BIT_OFFSET 2
#define VER_HDF5 3

static void
add_to_NCList(NC *ncp)
{
	assert(ncp != NULL);

	ncp->prev = NULL;
	if(NClist != NULL)
		NClist->prev = ncp;
	ncp->next = NClist;
	NClist = ncp;
}

static void
del_from_NCList(NC *ncp)
{
	assert(ncp != NULL);

	if(NClist == ncp)
	{
		assert(ncp->prev == NULL);
		NClist = ncp->next;
	}
	else
	{
		assert(ncp->prev != NULL);
		ncp->prev->next = ncp->next;
	}

	if(ncp->next != NULL)
		ncp->next->prev = ncp->prev;

	ncp->next = NULL;
	ncp->prev = NULL;
}


int
NC_check_id(int ncid, NC **ncpp)
{
	NC *ncp;

	if(ncid >= 0)
	{
		for(ncp = NClist; ncp != NULL; ncp = ncp->next)
		{
			if(ncp->nciop->fd == ncid)
			{
				*ncpp = ncp;
				return NC_NOERR; /* normal return */
			}
		}
	}

	/* else, not found */
	return NC_EBADID;
}


static void
free_NC(NC *ncp)
{
	if(ncp == NULL)
		return;
	free_NC_dimarrayV(&ncp->dims);
	free_NC_attrarrayV(&ncp->attrs);
	free_NC_vararrayV(&ncp->vars);
#if _CRAYMPP && defined(LOCKNUMREC)
	shfree(ncp);
#else
	free(ncp);
#endif /* _CRAYMPP && LOCKNUMREC */
}


static NC *
new_NC(const size_t *chunkp)
{
	NC *ncp;

#if _CRAYMPP && defined(LOCKNUMREC)
	ncp = (NC *) shmalloc(sizeof(NC));
#else
	ncp = (NC *) malloc(sizeof(NC));
#endif /* _CRAYMPP && LOCKNUMREC */
	if(ncp == NULL)
		return NULL;
	(void) memset(ncp, 0, sizeof(NC));

	ncp->xsz = MIN_NC_XSZ;
	assert(ncp->xsz == ncx_len_NC(ncp,0));
	
	ncp->chunk = chunkp != NULL ? *chunkp : NC_SIZEHINT_DEFAULT;

	return ncp;
}


static NC *
dup_NC(const NC *ref)
{
	NC *ncp;

#if _CRAYMPP && defined(LOCKNUMREC)
	ncp = (NC *) shmalloc(sizeof(NC));
#else
	ncp = (NC *) malloc(sizeof(NC));
#endif /* _CRAYMPP && LOCKNUMREC */
	if(ncp == NULL)
		return NULL;
	(void) memset(ncp, 0, sizeof(NC));

	if(dup_NC_dimarrayV(&ncp->dims, &ref->dims) != NC_NOERR)
		goto err;
	if(dup_NC_attrarrayV(&ncp->attrs, &ref->attrs) != NC_NOERR)
		goto err;
	if(dup_NC_vararrayV(&ncp->vars, &ref->vars) != NC_NOERR)
		goto err;

	ncp->xsz = ref->xsz;
	ncp->begin_var = ref->begin_var;
	ncp->begin_rec = ref->begin_rec;
	ncp->recsize = ref->recsize;
	NC_set_numrecs(ncp, NC_get_numrecs(ref));
	return ncp;
err:
	free_NC(ncp);
	return NULL;
}


/*
 *  Verify that this is a user nc_type
 * Formerly
NCcktype()
 * Sense of the return is changed.
 */
int
nc_cktype(nc_type type)
{
	switch((int)type){
	case NC_BYTE:
	case NC_CHAR:
	case NC_SHORT:
	case NC_INT:
	case NC_FLOAT:
	case NC_DOUBLE:
		return(NC_NOERR);
	}
	return(NC_EBADTYPE);
}


/*
 * How many objects of 'type'
 * will fit into xbufsize?
 */
size_t
ncx_howmany(nc_type type, size_t xbufsize)
{
	switch(type){
	case NC_BYTE:
	case NC_CHAR:
		return xbufsize;
	case NC_SHORT:
		return xbufsize/X_SIZEOF_SHORT;
	case NC_INT:
		return xbufsize/X_SIZEOF_INT;
	case NC_FLOAT:
		return xbufsize/X_SIZEOF_FLOAT;
	case NC_DOUBLE:
		return xbufsize/X_SIZEOF_DOUBLE;
	}
	assert("ncx_howmany: Bad type" == 0);
	return(0);
}

#define	D_RNDUP(x, align) _RNDUP(x, (off_t)(align))

/*
 * Compute each variable's 'begin' offset,
 * update 'begin_rec' as well.
 */
static int
NC_begins(NC *ncp,
	size_t h_minfree, size_t v_align,
	size_t v_minfree, size_t r_align)
{
	size_t ii;
	int sizeof_off_t;
	off_t index = 0;
	NC_var **vpp;
	NC_var *last = NULL;

	if(v_align == NC_ALIGN_CHUNK)
		v_align = ncp->chunk;
	if(r_align == NC_ALIGN_CHUNK)
		r_align = ncp->chunk;

	if (fIsSet(ncp->flags, NC_64BIT_OFFSET)) {
	  sizeof_off_t = 8;
	} else {
	  sizeof_off_t = 4;
	}
	
	ncp->xsz = ncx_len_NC(ncp,sizeof_off_t);

	if(ncp->vars.nelems == 0) 
		return NC_NOERR;

	/* only (re)calculate begin_var if there is not sufficient space in header
	   or start of non-record variables is not aligned as requested by valign */
	if (ncp->begin_var < ncp->xsz + h_minfree ||
	    ncp->begin_var != D_RNDUP(ncp->begin_var, v_align) ) 
	{
	  index = (off_t) ncp->xsz;
	  ncp->begin_var = D_RNDUP(index, v_align);
	  if(ncp->begin_var < index + h_minfree)
	  {
	    ncp->begin_var = D_RNDUP(index + (off_t)h_minfree, v_align);
	  }
	}
	index = ncp->begin_var;

	/* loop thru vars, first pass is for the 'non-record' vars */
	vpp = ncp->vars.value;
	for(ii = 0; ii < ncp->vars.nelems ; ii++, vpp++)
	{
		if( IS_RECVAR(*vpp) )
		{
			/* skip record variables on this pass */
			continue;
		}
#if 0
fprintf(stderr, "    VAR %d %s: %ld\n", ii, (*vpp)->name->cp, (long)index);
#endif
                if( sizeof_off_t == 4 && (index > X_OFF_MAX || index < 0) ) 
		{
		    return NC_EVARSIZE;
                }
		(*vpp)->begin = index;
		index += (*vpp)->len;
	}

	/* only (re)calculate begin_rec if there is not sufficient
	   space at end of non-record variables or if start of record
	   variables is not aligned as requested by r_align */
	if (ncp->begin_rec < index + v_minfree ||
	    ncp->begin_rec != D_RNDUP(ncp->begin_rec, r_align) )
	{
	  ncp->begin_rec = D_RNDUP(index, r_align);
	  if(ncp->begin_rec < index + v_minfree)
	  {
	    ncp->begin_rec = D_RNDUP(index + (off_t)v_minfree, r_align);
	  }
	}
	index = ncp->begin_rec;

	ncp->recsize = 0;

	/* loop thru vars, second pass is for the 'record' vars */
	vpp = (NC_var **)ncp->vars.value;
	for(ii = 0; ii < ncp->vars.nelems; ii++, vpp++)
	{
		if( !IS_RECVAR(*vpp) )
		{
			/* skip non-record variables on this pass */
			continue;
		}

#if 0
fprintf(stderr, "    REC %d %s: %ld\n", ii, (*vpp)->name->cp, (long)index);
#endif
                if( sizeof_off_t == 4 && (index > X_OFF_MAX || index < 0) ) 
		{
		    return NC_EVARSIZE;
                }
		(*vpp)->begin = index;
		index += (*vpp)->len;
		/* check if record size must fit in 32-bits */
#if SIZEOF_OFF_T == SIZEOF_SIZE_T && SIZEOF_SIZE_T == 4
		if( ncp->recsize > X_UINT_MAX - (*vpp)->len )
		{
		    return NC_EVARSIZE;
		}
#endif
		ncp->recsize += (*vpp)->len;
		last = (*vpp);
	}

	/*
	 * for special case of exactly one record variable, pack value
	 */
	if(last != NULL && ncp->recsize == last->len)
		ncp->recsize = *last->dsizes * last->xsz;

	if(NC_IsNew(ncp))
		NC_set_numrecs(ncp, 0);
	return NC_NOERR;
}


/*
 * Read just the numrecs member.
 * (A relatively expensive way to do things.)
 */
int
read_numrecs(NC *ncp)
{
	int status = NC_NOERR;
	const void *xp = NULL;
	size_t nrecs = NC_get_numrecs(ncp);

	assert(!NC_indef(ncp));

#define NC_NUMRECS_OFFSET 4
#define NC_NUMRECS_EXTENT 4
	status = ncp->nciop->get(ncp->nciop,
		 NC_NUMRECS_OFFSET, NC_NUMRECS_EXTENT, 0, (void **)&xp);
					/* cast away const */
	if(status != NC_NOERR)
		return status;

	status = ncx_get_size_t(&xp, &nrecs);

	(void) ncp->nciop->rel(ncp->nciop, NC_NUMRECS_OFFSET, 0);

	if(status == NC_NOERR)
	{
		NC_set_numrecs(ncp, nrecs);
		fClr(ncp->flags, NC_NDIRTY);
	}

	return status;
}


/*
 * Write out just the numrecs member.
 * (A relatively expensive way to do things.)
 */
int
write_numrecs(NC *ncp)
{
	int status = NC_NOERR;
	void *xp = NULL;

	assert(!NC_readonly(ncp));
	assert(!NC_indef(ncp));

	status = ncp->nciop->get(ncp->nciop,
		 NC_NUMRECS_OFFSET, NC_NUMRECS_EXTENT, RGN_WRITE, &xp);
	if(status != NC_NOERR)
		return status;

	{
		const size_t nrecs = NC_get_numrecs(ncp);
		status = ncx_put_size_t(&xp, &nrecs);
	}

	(void) ncp->nciop->rel(ncp->nciop, NC_NUMRECS_OFFSET, RGN_MODIFIED);

	if(status == NC_NOERR)
		fClr(ncp->flags, NC_NDIRTY);

	return status;
}


/*
 * Read in the header
 * It is expensive.
 */
static int
read_NC(NC *ncp)
{
	int status = NC_NOERR;

	free_NC_dimarrayV(&ncp->dims);
	free_NC_attrarrayV(&ncp->attrs);
	free_NC_vararrayV(&ncp->vars);

	status = nc_get_NC(ncp);

	if(status == NC_NOERR)
		fClr(ncp->flags, NC_NDIRTY | NC_HDIRTY);

	return status;
}


/*
 * Write out the header
 */
static int
write_NC(NC *ncp)
{
	int status = NC_NOERR;

	assert(!NC_readonly(ncp));

	status = ncx_put_NC(ncp, NULL, 0, 0);

	if(status == NC_NOERR)
		fClr(ncp->flags, NC_NDIRTY | NC_HDIRTY);

	return status;
}


/*
 * Write the header or the numrecs if necessary.
 */
int
NC_sync(NC *ncp)
{
	assert(!NC_readonly(ncp));

	if(NC_hdirty(ncp))
	{
		return write_NC(ncp);
	}
	/* else */

	if(NC_ndirty(ncp))
	{
		return write_numrecs(ncp);
	}
	/* else */

	return NC_NOERR;
}


/*
 * Initialize the 'non-record' variables.
 */
static int
fillerup(NC *ncp)
{
	int status = NC_NOERR;
	size_t ii;
	NC_var **varpp;

	assert(!NC_readonly(ncp));
	assert(NC_dofill(ncp));

	/* loop thru vars */
	varpp = ncp->vars.value;
	for(ii = 0; ii < ncp->vars.nelems; ii++, varpp++)
	{
		if(IS_RECVAR(*varpp))
		{
			/* skip record variables */
			continue;
		}

		status = fill_NC_var(ncp, *varpp, (*varpp)->len, 0);
		if(status != NC_NOERR)
			break;
	}
	return status;
}

/* Begin endef */

/*
 */
static int
fill_added_recs(NC *gnu, NC *old)
{
	NC_var ** const gnu_varpp = (NC_var **)gnu->vars.value;
	NC_var *const *const gnu_end = &gnu_varpp[gnu->vars.nelems];

	const int old_nrecs = (int) NC_get_numrecs(old);
	int recno = 0;
	NC_var **vpp = gnu_varpp;
	NC_var *const *const end = &vpp[gnu->vars.nelems];
	int numrecvars = 0;
	NC_var *recvarp = NULL;

	/* Determine if there is only one record variable.  If so, we
	   must treat as a special case because there's no record padding */
	for(; vpp < end; vpp++) {
	    if(IS_RECVAR(*vpp)) {
		recvarp = *vpp;
		numrecvars++;
	    }
	}

	for(; recno < old_nrecs; recno++)
	    {
		int varid = (int)old->vars.nelems;
		for(; varid < (int)gnu->vars.nelems; varid++)
		    {
			const NC_var *const gnu_varp = *(gnu_varpp + varid);
			if(!IS_RECVAR(gnu_varp))
			    {
				/* skip non-record variables */
				continue;
			    }
			/* else */
			{
			    size_t varsize = numrecvars == 1 ? gnu->recsize :  gnu_varp->len;
			    const int status = fill_NC_var(gnu, gnu_varp, varsize, recno);
			    if(status != NC_NOERR)
				return status;
			}
		    }
	    }
	return NC_NOERR;
}

/*
 */