putget.c

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

/* Do not edit this file. It is produced from the corresponding .m4 source */
/*
 *	Copyright 1996, University Corporation for Atmospheric Research
 *      See netcdf/COPYRIGHT file for copying and redistribution conditions.
 */
/* $Id: putget.c 2501 2007-11-20 02:33:29Z benkirk $ */

#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



/*
 * Next 6 type specific functions
 * Fill a some memory with the default special value.
 * Formerly
NC_arrayfill()
 */
static int
NC_fill_schar(
	void **xpp,
	size_t nelems)	/* how many */
{
	schar fillp[NFILL * sizeof(double)/X_SIZEOF_CHAR];

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

	{
		schar *vp = fillp;	/* lower bound of area to be filled */
		const schar *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_BYTE;
		}
	}
	return ncx_putn_schar_schar(xpp, nelems, fillp);
}

static int
NC_fill_char(
	void **xpp,
	size_t nelems)	/* how many */
{
	char fillp[NFILL * sizeof(double)/X_SIZEOF_CHAR];

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

	{
		char *vp = fillp;	/* lower bound of area to be filled */
		const char *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_CHAR;
		}
	}
	return ncx_putn_char_char(xpp, nelems, fillp);
}

static int
NC_fill_short(
	void **xpp,
	size_t nelems)	/* how many */
{
	short fillp[NFILL * sizeof(double)/X_SIZEOF_SHORT];

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

	{
		short *vp = fillp;	/* lower bound of area to be filled */
		const short *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_SHORT;
		}
	}
	return ncx_putn_short_short(xpp, nelems, fillp);
}


#if (SIZEOF_INT >= X_SIZEOF_INT)
static int
NC_fill_int(
	void **xpp,
	size_t nelems)	/* how many */
{
	int fillp[NFILL * sizeof(double)/X_SIZEOF_INT];

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

	{
		int *vp = fillp;	/* lower bound of area to be filled */
		const int *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_INT;
		}
	}
	return ncx_putn_int_int(xpp, nelems, fillp);
}

#elif SIZEOF_LONG == X_SIZEOF_INT
static int
NC_fill_int(
	void **xpp,
	size_t nelems)	/* how many */
{
	long fillp[NFILL * sizeof(double)/X_SIZEOF_INT];

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

	{
		long *vp = fillp;	/* lower bound of area to be filled */
		const long *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_INT;
		}
	}
	return ncx_putn_int_long(xpp, nelems, fillp);
}

#else
#error "NC_fill_int implementation"
#endif

static int
NC_fill_float(
	void **xpp,
	size_t nelems)	/* how many */
{
	float fillp[NFILL * sizeof(double)/X_SIZEOF_FLOAT];

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

	{
		float *vp = fillp;	/* lower bound of area to be filled */
		const float *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_FLOAT;
		}
	}
	return ncx_putn_float_float(xpp, nelems, fillp);
}

static int
NC_fill_double(
	void **xpp,
	size_t nelems)	/* how many */
{
	double fillp[NFILL * sizeof(double)/X_SIZEOF_DOUBLE];

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

	{
		double *vp = fillp;	/* lower bound of area to be filled */
		const double *const end = vp + nelems;
		while(vp < end)
		{
			*vp++ = NC_FILL_DOUBLE;
		}
	}
	return ncx_putn_double_double(xpp, nelems, fillp);
}





/* 
 * 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))