putget.m4

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

		lcoord *= varp->xsz;
		
		if(IS_RECVAR(varp))
			lcoord += (off_t)(*coord) * ncp->recsize;
		
		lcoord += varp->begin;
		return lcoord;
	}
}


dnl
dnl Output 'nelems' items of contiguous data of type "Type"
dnl for variable 'varp' at 'start'.
dnl "Xtype" had better match 'varp->type'.
dnl---
dnl
dnl PUTNCVX(Xtype, Type)
dnl
define(`PUTNCVX',dnl
`dnl
static int
putNCvx_$1_$2(NC *ncp, const NC_var *varp,
		 const size_t *start, size_t nelems, const $2 *value)
{
	off_t offset = NC_varoffset(ncp, varp, start);
	size_t remaining = varp->xsz * nelems;
	int status = NC_NOERR;
	void *xp;

	if(nelems == 0)
		return NC_NOERR;

	assert(value != NULL);

	for(;;)
	{
		size_t extent = MIN(remaining, ncp->chunk);
		size_t nput = ncx_howmany(varp->type, extent);

		int lstatus = ncp->nciop->get(ncp->nciop, offset, extent,
				 RGN_WRITE, &xp);	
		if(lstatus != NC_NOERR)
			return lstatus;
		
		lstatus = ncx_putn_$1_$2(&xp, nput, value);
		if(lstatus != NC_NOERR && status == NC_NOERR)
		{
			/* not fatal to the loop */
			status = lstatus;
		}

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

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

	}

	return status;
}
')dnl

PUTNCVX(char, char)

PUTNCVX(schar, schar)
PUTNCVX(schar, uchar)
PUTNCVX(schar, short)
PUTNCVX(schar, int)
PUTNCVX(schar, long)
PUTNCVX(schar, float)
PUTNCVX(schar, double)

PUTNCVX(short, schar)
PUTNCVX(short, uchar)
PUTNCVX(short, short)
PUTNCVX(short, int)
PUTNCVX(short, long)
PUTNCVX(short, float)
PUTNCVX(short, double)

PUTNCVX(int, schar)
PUTNCVX(int, uchar)
PUTNCVX(int, short)
PUTNCVX(int, int)
PUTNCVX(int, long)
PUTNCVX(int, float)
PUTNCVX(int, double)

PUTNCVX(float, schar)
PUTNCVX(float, uchar)
PUTNCVX(float, short)
PUTNCVX(float, int)
PUTNCVX(float, long)
PUTNCVX(float, float)
PUTNCVX(float, double)

PUTNCVX(double, schar)
PUTNCVX(double, uchar)
PUTNCVX(double, short)
PUTNCVX(double, int)
PUTNCVX(double, long)
PUTNCVX(double, float)
PUTNCVX(double, double)


dnl
dnl PUTNCV(Type)
dnl
define(`PUTNCV',dnl
`dnl
static int
putNCv_$1(NC *ncp, const NC_var *varp,
		 const size_t *start, size_t nelems, const $1 *value)
{
	switch(varp->type){
	case NC_CHAR:
		return NC_ECHAR;
	case NC_BYTE:
		return putNCvx_schar_$1(ncp, varp, start, nelems,
			value);
	case NC_SHORT:
		return putNCvx_short_$1(ncp, varp, start, nelems,
			value);
	case NC_INT:
		return putNCvx_int_$1(ncp, varp, start, nelems,
			value);
	case NC_FLOAT:
		return putNCvx_float_$1(ncp, varp, start, nelems,
			value);
	case NC_DOUBLE: 
		return putNCvx_double_$1(ncp, varp, start, nelems,
			value);
	}
	return NC_EBADTYPE;
}
')dnl

static int
putNCv_text(NC *ncp, const NC_var *varp,
		 const size_t *start, size_t nelems, const char *value)
{
	if(varp->type != NC_CHAR)
		return NC_ECHAR;
	return putNCvx_char_char(ncp, varp, start, nelems, value);
}

PUTNCV(schar)
PUTNCV(uchar)
PUTNCV(short)
PUTNCV(int)
PUTNCV(long)
PUTNCV(float)
PUTNCV(double)


dnl
dnl GETNCVX(XType, Type)
dnl
define(`GETNCVX',dnl
`dnl
static int
getNCvx_$1_$2(const NC *ncp, const NC_var *varp,
		 const size_t *start, size_t nelems, $2 *value)
{
	off_t offset = NC_varoffset(ncp, varp, start);
	size_t remaining = varp->xsz * nelems;
	int status = NC_NOERR;
	const void *xp;

	if(nelems == 0)
		return NC_NOERR;

	assert(value != NULL);

	for(;;)
	{
		size_t extent = MIN(remaining, ncp->chunk);
		size_t nget = ncx_howmany(varp->type, extent);

		int lstatus = ncp->nciop->get(ncp->nciop, offset, extent,
				 0, (void **)&xp);	/* cast away const */
		if(lstatus != NC_NOERR)
			return lstatus;
		
		lstatus = ncx_getn_$1_$2(&xp, nget, value);
		if(lstatus != NC_NOERR && status == NC_NOERR)
			status = lstatus;

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

		remaining -= extent;
		if(remaining == 0)
			break; /* normal loop exit */
		offset += extent;
		value += nget;
	}

	return status;
}
')dnl

GETNCVX(char, char)

GETNCVX(schar, schar)
GETNCVX(schar, uchar)
GETNCVX(schar, short)
GETNCVX(schar, int)
GETNCVX(schar, long)
GETNCVX(schar, float)
GETNCVX(schar, double)

GETNCVX(short, schar)
GETNCVX(short, uchar)
GETNCVX(short, short)
GETNCVX(short, int)
GETNCVX(short, long)
GETNCVX(short, float)
GETNCVX(short, double)

GETNCVX(int, schar)
GETNCVX(int, uchar)
GETNCVX(int, short)
GETNCVX(int, int)
GETNCVX(int, long)
GETNCVX(int, float)
GETNCVX(int, double)

GETNCVX(float, schar)
GETNCVX(float, uchar)
GETNCVX(float, short)
GETNCVX(float, int)
GETNCVX(float, long)
GETNCVX(float, float)
GETNCVX(float, double)

GETNCVX(double, schar)
GETNCVX(double, uchar)
GETNCVX(double, short)
GETNCVX(double, int)
GETNCVX(double, long)
GETNCVX(double, float)
GETNCVX(double, double)


dnl
dnl GETNCV(Type)
dnl
define(`GETNCV',dnl
`dnl
static int
getNCv_$1(const NC *ncp, const NC_var *varp,
		 const size_t *start, size_t nelems, $1 *value)
{
	switch(varp->type){
	case NC_CHAR:
		return NC_ECHAR;
	case NC_BYTE:
		return getNCvx_schar_$1(ncp, varp, start, nelems,
			value);
	case NC_SHORT:
		return getNCvx_short_$1(ncp, varp, start, nelems,
			value);
	case NC_INT:
		return getNCvx_int_$1(ncp, varp, start, nelems,
			value);
	case NC_FLOAT:
		return getNCvx_float_$1(ncp, varp, start, nelems,
			value);
	case NC_DOUBLE: 
		return getNCvx_double_$1(ncp, varp, start, nelems,
			value);
	}
	return NC_EBADTYPE;
}
')dnl

GETNCV(schar)
GETNCV(uchar)
GETNCV(short)
GETNCV(int)
GETNCV(long)
GETNCV(float)
GETNCV(double)


static int
getNCv_text(const NC *ncp, const NC_var *varp,
		 const size_t *start, size_t nelems, char *value)
{
	if(varp->type != NC_CHAR)
		return NC_ECHAR;
	return getNCvx_char_char(ncp, varp, start, nelems, value);
}


/*
 * Copy 'nbytes' contiguous external values
 * from ('inncp', invp', inncoord')
 * to   ('outncp', 'outvp', 'outcoord')
 * 'inncp' shouldn't be the same as 'outncp'.
 * Used only by ncvarcopy()
 */
static int
NCxvarcpy(NC *inncp, NC_var *invp, size_t *incoord,
	NC *outncp, NC_var *outvp, size_t *outcoord, size_t nbytes)
{
	int status;
	off_t inoffset = NC_varoffset(inncp, invp, incoord);
	off_t outoffset = NC_varoffset(outncp, outvp, outcoord);
	void *inxp;
	void *outxp;
	const size_t chunk = MIN(inncp->chunk, outncp->chunk);

	do {
		const size_t extent = MIN(nbytes, chunk);

		status = inncp->nciop->get(inncp->nciop, inoffset, extent,
				 0, &inxp);	
		if(status != NC_NOERR)
			return status;

		status = outncp->nciop->get(outncp->nciop, outoffset, extent,
				 RGN_WRITE, &outxp);	
		if(status != NC_NOERR)
		{
			(void) inncp->nciop->rel(inncp->nciop, inoffset, 0);	
			break;
		}

		(void) memcpy(outxp, inxp, extent);

		status = outncp->nciop->rel(outncp->nciop, outoffset,
			 RGN_MODIFIED);
		(void) inncp->nciop->rel(inncp->nciop, inoffset, 0);	

		nbytes -= extent;
		if(nbytes == 0)
			break; /* normal loop exit */
		inoffset += extent;
		outoffset += extent;
		
	} while (status == NC_NOERR);

	return status;
}


/*
 *  For ncvar{put,get},
 *  find the largest contiguous block from within 'edges'.
 *  returns the index to the left of this (which may be -1).
 *  Compute the number of contiguous elements and return
 *  that in *iocountp.
 *  The presence of "record" variables makes this routine
 *  overly subtle.
 */
static int
NCiocount(const NC *const ncp, const NC_var *const varp,
	const size_t *const edges,
	size_t *const iocountp)
{
	const size_t *edp0 = edges;
	const size_t *edp = edges + varp->ndims;
	const size_t *shp = varp->shape + varp->ndims;

	if(IS_RECVAR(varp))
	{
		if(varp->ndims == 1 && ncp->recsize <= varp->len)
		{
			/* one dimensional && the only 'record' variable */
			*iocountp = *edges;
			return(0);
		}
		/* else */
		edp0++;
	}

	assert(edges != NULL);

	/* find max contiguous */
	while(edp > edp0)
	{
		shp--; edp--;
		if(*edp < *shp )
		{
			const size_t *zedp = edp;
			while(zedp >= edp0)
			{
				if(*zedp == 0)
				{
					*iocountp = 0;
					goto done;
				}
				/* Tip of the hat to segmented architectures */
				if(zedp == edp0)
					break;
				zedp--;
			}
			break;
		}
		assert(*edp == *shp);
	}

	/*
	 * edp, shp reference rightmost index s.t. *(edp +1) == *(shp +1)
	 *
	 * Or there is only one dimension.
	 * If there is only one dimension and it is 'non record' dimension,
	 * 	edp is &edges[0] and we will return -1.
	 * If there is only one dimension and and it is a "record dimension",
	 *	edp is &edges[1] (out of bounds) and we will return 0;
	 */
	assert(shp >= varp->shape + varp->ndims -1 
		|| *(edp +1) == *(shp +1));

	/* now accumulate max count for a single io operation */
	for(*iocountp = 1, edp0 = edp;
		 	edp0 < edges + varp->ndims;
			edp0++)
	{
		*iocountp *= *edp0;
	}

done:
	return((int)(edp - edges) - 1);
}


/*
 * Set the elements of the array 'upp' to
 * the sum of the corresponding elements of
 * 'stp' and 'edp'. 'end' should be &stp[nelems].
 */
static void
set_upper(size_t *upp, /* modified on return */
	const size_t *stp,
	const size_t *edp,
	const size_t *const end)
{
	while(upp < end) {
		*upp++ = *stp++ + *edp++;
	}
}


/*
 * The infamous and oft-discussed odometer code.
 *
 * 'start[]' is the starting coordinate.
 * 'upper[]' is the upper bound s.t. start[ii] < upper[ii].
 * 'coord[]' is the register, the current coordinate value.
 * For some ii,
 * upp == &upper[ii]
 * cdp == &coord[ii]
 * 
 * Running this routine increments *cdp.
 *
 * If after the increment, *cdp is equal to *upp
 * (and cdp is not the leftmost dimension),
 * *cdp is "zeroed" to the starting value and
 * we need to "carry", eg, increment one place to
 * the left.
 * 
 * TODO: Some architectures hate recursion?
 * 	Reimplement non-recursively.
 */
static void
odo1(const size_t *const start, const size_t *const upper,
	size_t *const coord, /* modified on return */
	const size_t *upp,
	size_t *cdp)
{
	assert(coord <= cdp && cdp <= coord + NC_MAX_VAR_DIMS);
	assert(upper <= upp && upp <= upper + NC_MAX_VAR_DIMS);
	assert(upp - upper == cdp - coord);
	
	assert(*cdp <= *upp);

	(*cdp)++;
	if(cdp != coord && *cdp >= *upp)
	{
		*cdp = start[cdp - coord];
		odo1(start, upper, coord, upp -1, cdp -1);
	}
}
#ifdef _CRAYC
#pragma _CRI noinline odo1
#endif


dnl
dnl NCTEXTCOND(Abbrv)
dnl This is used inside the NC{PUT,GET} macros below
dnl
define(`NCTEXTCOND',dnl
`dnl
ifelse($1, text,dnl
`dnl
	if(varp->type != NC_CHAR)
		return NC_ECHAR;
',dnl
`dnl
	if(varp->type == NC_CHAR)
		return NC_ECHAR;
')dnl
')dnl


/* Public */

dnl
dnl NCPUTVAR1(Abbrev, Type)
dnl
define(`NCPUTVAR1',dnl
`dnl
int
nc_put_var1_$1(int ncid, int varid, const size_t *coord,
	const $2 *value)
{
	int status;
	NC *ncp;
	const NC_var *varp;

	status = NC_check_id(ncid, &ncp); 
	if(status != NC_NOERR)
		return status;

	if(NC_readonly(ncp))
		return NC_EPERM;

	if(NC_indef(ncp))
		return NC_EINDEFINE;

	varp = NC_lookupvar(ncp, varid);
	if(varp == NULL)
		return NC_ENOTVAR; /* TODO: lost NC_EGLOBAL */

NCTEXTCOND($1)
	status = NCcoordck(ncp, varp, coord);
	if(status != NC_NOERR)
		return status;

	if(IS_RECVAR(varp))
	{
		status = NCvnrecs(ncp, *coord +1);
		if(status != NC_NOERR)
			return status;
	}

	return putNCv_$1(ncp, varp, coord, 1, value);
}
')dnl

NCPUTVAR1(text, char)

NCPUTVAR1(uchar, uchar)
NCPUTVAR1(schar, schar)
NCPUTVAR1(short, short)
NCPUTVAR1(int, int)
NCPUTVAR1(long, long)
NCPUTVAR1(float, float)
NCPUTVAR1(double, double)

dnl
dnl NCGETVAR1(Abbrv, Type)
dnl
define(`NCGETVAR1',dnl
`dnl
int
nc_get_var1_$1(int ncid, int varid, const size_t *coord, $2 *value)
{
	int status;
	NC *ncp;
	const NC_var *varp;

	status = NC_check_id(ncid, &ncp); 
	if(status != NC_NOERR)
		return status;

	if(NC_indef(ncp))
		return NC_EINDEFINE;

	varp = NC_lookupvar(ncp, varid);
	if(varp == NULL)
		return NC_ENOTVAR; /* TODO: lost NC_EGLOBAL */

NCTEXTCOND($1)
	status = NCcoordck(ncp, varp, coord);
	if(status != NC_NOERR)
		return status;

	return getNCv_$1(ncp, varp, coord, 1, value);
}
')dnl

NCGETVAR1(text, char)

NCGETVAR1(uchar, uchar)
NCGETVAR1(schar, schar)
NCGETVAR1(short, short)
NCGETVAR1(int, int)
NCGETVAR1(long, long)
NCGETVAR1(float, float)
NCGETVAR1(double, double)

dnl
dnl NCPUTVARA(Abbrv, Type)
dnl
define(`NCPUTVARA',dnl
`dnl
int
nc_put_vara_$1(int ncid, int varid,
	 const size_t *start, const size_t *edges, const $2 *value)
{
	int status = NC_NOERR;
	NC *ncp;
	const NC_var *varp;
	int ii;
	size_t iocount;

	status = NC_check_id(ncid, &ncp); 
	if(status != NC_NOERR)
		return status;

	if(NC_readonly(ncp))
		return NC_EPERM;

	if(NC_indef(ncp))
		return NC_EINDEFINE;

	varp = NC_lookupvar(ncp, varid);
	if(varp == NULL)
		return NC_ENOTVAR; /* TODO: lost NC_EGLOBAL */

NCTEXTCOND($1)
	status = NCcoordck(ncp, varp, start);
	if(status != NC_NOERR)
		return status;
	status = NCedgeck(ncp, varp, start, edges);
	if(status != NC_NOERR)
		return status;