putget.m4

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

dnl This is m4 source.
dnl Process using m4 to produce 'C' language file.
dnl
undefine(`begin')dnl
undefine(`index')dnl
undefine(`len')dnl
dnl
dnl If you see this line, you can ignore the next one.
/* Do not edit this file. It is produced from the corresponding .m4 source */
dnl
/*
 *	Copyright 1996, University Corporation for Atmospheric Research
 *      See netcdf/COPYRIGHT file for copying and redistribution conditions.
 */
/* $Id: putget.m4,v 2.62 2006/06/20 15:05:14 ed Exp $ */

#include "nc.h"
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include "ncx.h"
#include "fbits.h"
#include "onstack.h"
#ifdef LOCKNUMREC
#  include <mpp/shmem.h>	/* for SGI/Cray SHMEM routines */
#  ifdef LN_TEST
#    include <stdio.h>
#  endif
#endif

#undef MIN  /* system may define MIN somewhere and complain */
#define MIN(mm,nn) (((mm) < (nn)) ? (mm) : (nn))

/* #define ODEBUG 1 */

#if ODEBUG
#include <stdio.h>
/*
 * Print the values of an array of size_t
 */
void
arrayp(const char *label, size_t count, const size_t *array)
{
	(void) fprintf(stderr, "%s", label);
	(void) fputc('\t',stderr);	
	for(; count > 0; count--, array++)
		(void) fprintf(stderr," %lu", (unsigned long)*array);
	(void) fputc('\n',stderr);	
}
#endif /* ODEBUG */

/*
 *  This is how much space is required by the user, as in
 *
 *   vals = malloc(nel * nctypelen(var.type));
 *   ncvarget(cdfid, varid, cor, edg, vals);
 */
int
nctypelen(nc_type type) 
{
	switch(type){
	case NC_BYTE :
	case NC_CHAR :
		return((int)sizeof(char));
	case NC_SHORT :
		return(int)(sizeof(short));
	case NC_INT :
		return((int)sizeof(int));
	case NC_FLOAT :
		return((int)sizeof(float));
	case NC_DOUBLE : 
		return((int)sizeof(double));
	}

	return -1;
}


/* Begin fill */
/*
 * This is tunable parameter.
 * It essentially controls the tradeoff between the number of times
 * memcpy() gets called to copy the external data to fill 
 * a large buffer vs the number of times its called to
 * prepare the external data.
 */
#if	_SX
/* NEC SX specific optimization */
#define	NFILL	2048
#else
#define	NFILL	16
#endif


dnl
dnl NCFILL(Type, Xtype, XSize, Fill)
dnl
define(`NCFILL',dnl
`dnl
static int
NC_fill_$2(
	void **xpp,
	size_t nelems)	/* how many */
{
	$1 fillp[NFILL * sizeof(double)/$3];

	assert(nelems <= sizeof(fillp)/sizeof(fillp[0]));

	{
		$1 *vp = fillp;	/* lower bound of area to be filled */
		const $1 *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = $4;
		}
	}
	return ncx_putn_$2_$1(xpp, nelems, fillp);
}
')dnl

/*
 * Next 6 type specific functions
 * Fill a some memory with the default special value.
 * Formerly
NC_arrayfill()
 */
NCFILL(schar, schar, X_SIZEOF_CHAR, NC_FILL_BYTE)
NCFILL(char, char, X_SIZEOF_CHAR, NC_FILL_CHAR)
NCFILL(short, short, X_SIZEOF_SHORT, NC_FILL_SHORT)

#if (SIZEOF_INT >= X_SIZEOF_INT)
NCFILL(int, int, X_SIZEOF_INT, NC_FILL_INT)
#elif SIZEOF_LONG == X_SIZEOF_INT
NCFILL(long, int, X_SIZEOF_INT, NC_FILL_INT)
#else
#error "NC_fill_int implementation"
#endif

NCFILL(float, float, X_SIZEOF_FLOAT, NC_FILL_FLOAT)
NCFILL(double, double, X_SIZEOF_DOUBLE, NC_FILL_DOUBLE)




/* 
 * Fill the external space for variable 'varp' values at 'recno' with
 * the appropriate value. If 'varp' is not a record variable, fill the
 * whole thing.  For the special case when 'varp' is the only record
 * variable and it is of type byte, char, or short, varsize should be
 * ncp->recsize, otherwise it should be varp->len.
 * Formerly
xdr_NC_fill()
 */
int
fill_NC_var(NC *ncp, const NC_var *varp, size_t varsize, size_t recno)
{
	char xfillp[NFILL * X_SIZEOF_DOUBLE];
	const size_t step = varp->xsz;
	const size_t nelems = sizeof(xfillp)/step;
	const size_t xsz = varp->xsz * nelems;
	NC_attr **attrpp = NULL;
	off_t offset;
	size_t remaining = varsize;

	void *xp;
	int status = NC_NOERR;

	/*
	 * Set up fill value
	 */
	attrpp = NC_findattr(&varp->attrs, _FillValue);
	if( attrpp != NULL )
	{
		/* User defined fill value */
		if( (*attrpp)->type != varp->type || (*attrpp)->nelems != 1 )
		{
			return NC_EBADTYPE;
		}
		else
		{
			/* Use the user defined value */
			char *cp = xfillp;
			const char *const end = &xfillp[sizeof(xfillp)];

			assert(step <= (*attrpp)->xsz);

			for( /*NADA*/; cp < end; cp += step)
			{
				(void) memcpy(cp, (*attrpp)->xvalue, step);
			}
		}
	}
	else
	{
		/* use the default */
		
		assert(xsz % X_ALIGN == 0);
		assert(xsz <= sizeof(xfillp));
	
		xp = xfillp;
	
		switch(varp->type){
		case NC_BYTE :
			status = NC_fill_schar(&xp, nelems);
			break;
		case NC_CHAR :
			status = NC_fill_char(&xp, nelems);
			break;
		case NC_SHORT :
			status = NC_fill_short(&xp, nelems);
			break;
		case NC_INT :
			status = NC_fill_int(&xp, nelems);
			break;
		case NC_FLOAT :
			status = NC_fill_float(&xp, nelems);
			break;
		case NC_DOUBLE : 
			status = NC_fill_double(&xp, nelems);
			break;
		default :
			assert("fill_NC_var invalid type" == 0);
			status = NC_EBADTYPE;
			break;
		}
		if(status != NC_NOERR)
			return status;
	
		assert(xp == xfillp + xsz);
	}

	/*
	 * copyout:
	 * xfillp now contains 'nelems' elements of the fill value
	 * in external representation.
	 */

	/*
	 * Copy it out.
	 */

	offset = varp->begin;
	if(IS_RECVAR(varp))
	{
		offset += (off_t)ncp->recsize * recno;
	}

	assert(remaining > 0);
	for(;;)
	{
		const size_t chunksz = MIN(remaining, ncp->chunk);
		size_t ii;

		status = ncp->nciop->get(ncp->nciop, offset, chunksz,
				 RGN_WRITE, &xp);	
		if(status != NC_NOERR)
		{
			return status;
		}

		/*
		 * fill the chunksz buffer in units  of xsz
		 */
		for(ii = 0; ii < chunksz/xsz; ii++)
		{
			(void) memcpy(xp, xfillp, xsz);
			xp = (char *)xp + xsz;
		}
		/*
		 * Deal with any remainder
		 */
		{
			const size_t rem = chunksz % xsz;
			if(rem != 0)
			{
				(void) memcpy(xp, xfillp, rem);
				/* xp = (char *)xp + xsz; */
			}

		}

		status = ncp->nciop->rel(ncp->nciop, offset, RGN_MODIFIED);

		if(status != NC_NOERR)
		{
			break;
		}

		remaining -= chunksz;
		if(remaining == 0)
			break;	/* normal loop exit */
		offset += chunksz;

	}

	return status;
}
/* End fill */


/*
 * Add a record containing the fill values.
 */
static int
NCfillrecord(NC *ncp, const NC_var *const *varpp, size_t recno)
{
	size_t ii = 0;
	for(; ii < ncp->vars.nelems; ii++, varpp++)
	{
		if( !IS_RECVAR(*varpp) )
		{
			continue;	/* skip non-record variables */
		}
		{
		const int status = fill_NC_var(ncp, *varpp, (*varpp)->len, recno);
		if(status != NC_NOERR)
			return status;
		}
	}
	return NC_NOERR;
}


/*
 * Add a record containing the fill values in the special case when
 * there is exactly one record variable, where we don't require each
 * record to be four-byte aligned (no record padding).
 */
static int
NCfillspecialrecord(NC *ncp, const NC_var *varp, size_t recno)
{
    int status;
    assert(IS_RECVAR(varp));
    status = fill_NC_var(ncp, varp, ncp->recsize, recno);
    if(status != NC_NOERR)
	return status;
    return NC_NOERR;
}


/*
 * It is advantageous to
 * #define TOUCH_LAST
 * when using memory mapped io.
 */
#if TOUCH_LAST
/*
 * Grow the file to a size which can contain recno
 */
static int
NCtouchlast(NC *ncp, const NC_var *const *varpp, size_t recno)
{
	int status = NC_NOERR;
	const NC_var *varp = NULL;
	
	{
	size_t ii = 0;
	for(; ii < ncp->vars.nelems; ii++, varpp++)
	{
		if( !IS_RECVAR(*varpp) )
		{
			continue;	/* skip non-record variables */
		}
		varp = *varpp;
	}
	}
	assert(varp != NULL);
	assert( IS_RECVAR(varp) );
	{
		const off_t offset = varp->begin
				+ (off_t)(recno-1) * (off_t)ncp->recsize
				+ (off_t)(varp->len - varp->xsz);
		void *xp;


		status = ncp->nciop->get(ncp->nciop, offset, varp->xsz,
				 RGN_WRITE, &xp);	
		if(status != NC_NOERR)
			return status;
		(void)memset(xp, 0, varp->xsz);
		status = ncp->nciop->rel(ncp->nciop, offset, RGN_MODIFIED);
	}
	return status;
}
#endif /* TOUCH_LAST */


/*
 * Ensure that the netcdf file has 'numrecs' records,
 * add records and fill as neccessary.
 */
static int
NCvnrecs(NC *ncp, size_t numrecs)
{
	int status = NC_NOERR;
#ifdef LOCKNUMREC
	ushmem_t myticket = 0, nowserving = 0;
	ushmem_t numpe = (ushmem_t) _num_pes();

	/* get ticket and wait */
	myticket = shmem_short_finc((shmem_t *) ncp->lock + LOCKNUMREC_LOCK,
		ncp->lock[LOCKNUMREC_BASEPE]);
#ifdef LN_TEST
		fprintf(stderr,"%d of %d : ticket = %hu\n",
			_my_pe(), _num_pes(), myticket);
#endif
	do {
		shmem_short_get((shmem_t *) &nowserving,
			(shmem_t *) ncp->lock + LOCKNUMREC_SERVING, 1,
			ncp->lock[LOCKNUMREC_BASEPE]);
#ifdef LN_TEST
		fprintf(stderr,"%d of %d : serving = %hu\n",
			_my_pe(), _num_pes(), nowserving);
#endif
		/* work-around for non-unique tickets */
		if (nowserving > myticket && nowserving < myticket + numpe ) {
			/* get a new ticket ... you've been bypassed */ 
			/* and handle the unlikely wrap-around effect */
			myticket = shmem_short_finc(
				(shmem_t *) ncp->lock + LOCKNUMREC_LOCK,
				ncp->lock[LOCKNUMREC_BASEPE]);
#ifdef LN_TEST
				fprintf(stderr,"%d of %d : new ticket = %hu\n",
					_my_pe(), _num_pes(), myticket);
#endif
		}
	} while(nowserving != myticket);
	/* now our turn to check & update value */
#endif

	if(numrecs > NC_get_numrecs(ncp))
	{


#if TOUCH_LAST
		status = NCtouchlast(ncp,
			(const NC_var *const*)ncp->vars.value,
			numrecs);
		if(status != NC_NOERR)
			goto common_return;
#endif /* TOUCH_LAST */

		set_NC_ndirty(ncp);

		if(!NC_dofill(ncp))
		{
			/* Simply set the new numrecs value */
			NC_set_numrecs(ncp, numrecs);
		}
		else
		{
		    /* Treat two cases differently: 
		        - exactly one record variable (no padding)
                        - multiple record variables (each record padded 
                          to 4-byte alignment)
		    */
		    NC_var **vpp = (NC_var **)ncp->vars.value;
		    NC_var *const *const end = &vpp[ncp->vars.nelems];
		    NC_var *recvarp = NULL;	/* last record var */
		    int numrecvars = 0;
		    size_t cur_nrecs;
		    
		    /* determine how many record variables */
		    for( /*NADA*/; vpp < end; vpp++) {
			if(IS_RECVAR(*vpp)) {
			    recvarp = *vpp;
			    numrecvars++;
			}
		    }
		    
		    if (numrecvars != 1) { /* usual case */
			/* Fill each record out to numrecs */
			while((cur_nrecs = NC_get_numrecs(ncp)) < numrecs)
			    {
				status = NCfillrecord(ncp,
					(const NC_var *const*)ncp->vars.value,
					cur_nrecs);
				if(status != NC_NOERR)
				{
					break;
				}
				NC_increase_numrecs(ncp, cur_nrecs +1);
			}
			if(status != NC_NOERR)
				goto common_return;
		    } else {	/* special case */
			/* Fill each record out to numrecs */
			while((cur_nrecs = NC_get_numrecs(ncp)) < numrecs)
			    {
				status = NCfillspecialrecord(ncp,
					recvarp,
					cur_nrecs);
				if(status != NC_NOERR)
				{
					break;
				}
				NC_increase_numrecs(ncp, cur_nrecs +1);
			}
			if(status != NC_NOERR)
				goto common_return;
			
		    }
		}

		if(NC_doNsync(ncp))
		{
			status = write_numrecs(ncp);
		}

	}
common_return:
#ifdef LOCKNUMREC
	/* finished with our lock - increment serving number */
	(void) shmem_short_finc((shmem_t *) ncp->lock + LOCKNUMREC_SERVING,
		ncp->lock[LOCKNUMREC_BASEPE]);
#endif
	return status;
}


/* 
 * Check whether 'coord' values are valid for the variable.
 */
static int
NCcoordck(NC *ncp, const NC_var *varp, const size_t *coord)
{
	const size_t *ip;
	size_t *up;

	if(varp->ndims == 0)
		return NC_NOERR;	/* 'scalar' variable */

	if(IS_RECVAR(varp))
	{
		if(*coord > X_INT_MAX)
			return NC_EINVALCOORDS; /* sanity check */
		if(NC_readonly(ncp) && *coord >= NC_get_numrecs(ncp))
		{
			if(!NC_doNsync(ncp))
				return NC_EINVALCOORDS;
			/* else */
			{
				/* Update from disk and check again */
				const int status = read_numrecs(ncp);
				if(status != NC_NOERR)
					return status;
				if(*coord >= NC_get_numrecs(ncp))
					return NC_EINVALCOORDS;
			}
		}
		ip = coord + 1;
		up = varp->shape + 1;
	}
	else
	{
		ip = coord;
		up = varp->shape;
	}
	
#ifdef CDEBUG
fprintf(stderr,"	NCcoordck: coord %ld, count %d, ip %ld\n",
		coord, varp->ndims, ip );
#endif /* CDEBUG */

	for(; ip < coord + varp->ndims; ip++, up++)
	{

#ifdef CDEBUG
fprintf(stderr,"	NCcoordck: ip %p, *ip %ld, up %p, *up %lu\n",
			ip, *ip, up, *up );
#endif /* CDEBUG */

		/* cast needed for braindead systems with signed size_t */
		if((unsigned long) *ip >= (unsigned long) *up )
			return NC_EINVALCOORDS;
	}

	return NC_NOERR;
}


/* 
 * Check whether 'edges' are valid for the variable and 'start'
 */
/*ARGSUSED*/
static int
NCedgeck(const NC *ncp, const NC_var *varp,
	 const size_t *start, const size_t *edges)
{
	const size_t *const end = start + varp->ndims;
	const size_t *shp = varp->shape;

	if(varp->ndims == 0)
		return NC_NOERR;	/* 'scalar' variable */

	if(IS_RECVAR(varp))
	{
		start++;
		edges++;
		shp++;
	}

	for(; start < end; start++, edges++, shp++)
	{
		/* cast needed for braindead systems with signed size_t */
		if((unsigned long) *edges > *shp ||
			(unsigned long) *start + (unsigned long) *edges > *shp)
		{
			return(NC_EEDGE);
		}
	}
	return NC_NOERR;
}


/* 
 * Translate the (variable, coord) pair into a seek index
 */
static off_t
NC_varoffset(const NC *ncp, const NC_var *varp, const size_t *coord)
{
	if(varp->ndims == 0) /* 'scalar' variable */
		return varp->begin;

	if(varp->ndims == 1)
	{
		if(IS_RECVAR(varp))
			return varp->begin +
				 (off_t)(*coord) * (off_t)ncp->recsize;
		/* else */
		return varp->begin + (off_t)(*coord) * (off_t)varp->xsz;
	}
	/* else */
	{
		off_t lcoord = (off_t)coord[varp->ndims -1];

		size_t *up = varp->dsizes +1;
		const size_t *ip = coord;
		const size_t *const end = varp->dsizes + varp->ndims;
		
		if(IS_RECVAR(varp))
			up++, ip++;

		for(; up < end; up++, ip++)
			lcoord += *up * *ip;