putget.m4

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

	if(varp->ndims == 0) /* scalar variable */
	{
		return( putNCv_$1(ncp, varp, start, 1, value) );
	}

	if(IS_RECVAR(varp))
	{
		status = NCvnrecs(ncp, *start + *edges);
		if(status != NC_NOERR)
			return status;

		if(varp->ndims == 1
			&& ncp->recsize <= varp->len)
		{
			/* one dimensional && the only record variable  */
			return( putNCv_$1(ncp, varp, start, *edges, value) );
		}
	}

	/*
	 * find max contiguous
	 *   and accumulate max count for a single io operation
	 */
	ii = NCiocount(ncp, varp, edges, &iocount);

	if(ii == -1)
	{
		return( putNCv_$1(ncp, varp, start, iocount, value) );
	}

	assert(ii >= 0);


	{ /* inline */
	ALLOC_ONSTACK(coord, size_t, varp->ndims);
	ALLOC_ONSTACK(upper, size_t, varp->ndims);
	const size_t index = ii;

	/* copy in starting indices */
	(void) memcpy(coord, start, varp->ndims * sizeof(size_t));

	/* set up in maximum indices */
	set_upper(upper, start, edges, &upper[varp->ndims]);

	/* ripple counter */
	while(*coord < *upper)
	{
		const int lstatus = putNCv_$1(ncp, varp, coord, iocount,
				 value);
		if(lstatus != NC_NOERR)
		{
			if(lstatus != NC_ERANGE)
			{
				status = lstatus;
				/* fatal for the loop */
				break;
			}
			/* else NC_ERANGE, not fatal for the loop */
			if(status == NC_NOERR)
				status = lstatus;
		}
		value += iocount;
		odo1(start, upper, coord, &upper[index], &coord[index]);
	}

	FREE_ONSTACK(upper);
	FREE_ONSTACK(coord);
	} /* end inline */

	return status;
}
')dnl

NCPUTVARA(text, char)

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


dnl
dnl NCGETVARA(Abbrv, Type)
dnl
define(`NCGETVARA',dnl
`dnl
int
nc_get_vara_$1(int ncid, int varid,
	 const size_t *start, const size_t *edges, $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_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;

	if(varp->ndims == 0) /* scalar variable */
	{
		return( getNCv_$1(ncp, varp, start, 1, value) );
	}

	if(IS_RECVAR(varp))
	{
		if(*start + *edges > NC_get_numrecs(ncp))
			return NC_EEDGE;
		if(varp->ndims == 1 && ncp->recsize <= varp->len)
		{
			/* one dimensional && the only record variable  */
			return( getNCv_$1(ncp, varp, start, *edges, value) );
		}
	}

	/*
	 * find max contiguous
	 *   and accumulate max count for a single io operation
	 */
	ii = NCiocount(ncp, varp, edges, &iocount);

	if(ii == -1)
	{
		return( getNCv_$1(ncp, varp, start, iocount, value) );
	}

	assert(ii >= 0);


	{ /* inline */
	ALLOC_ONSTACK(coord, size_t, varp->ndims);
	ALLOC_ONSTACK(upper, size_t, varp->ndims);
	const size_t index = ii;

	/* copy in starting indices */
	(void) memcpy(coord, start, varp->ndims * sizeof(size_t));

	/* set up in maximum indices */
	set_upper(upper, start, edges, &upper[varp->ndims]);

	/* ripple counter */
	while(*coord < *upper)
	{
		const int lstatus = getNCv_$1(ncp, varp, coord, iocount,
				value);
		if(lstatus != NC_NOERR)
		{
			if(lstatus != NC_ERANGE)
			{
				status = lstatus;
				/* fatal for the loop */
				break;
			}
			/* else NC_ERANGE, not fatal for the loop */
			if(status == NC_NOERR)
				status = lstatus;
		}
		value += iocount;
		odo1(start, upper, coord, &upper[index], &coord[index]);
	}

	FREE_ONSTACK(upper);
	FREE_ONSTACK(coord);
	} /* end inline */

	return status;
}
')dnl

NCGETVARA(text, char)

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


#if defined(__cplusplus)
/* C++ consts default to internal linkage and must be initialized */
const size_t coord_zero[NC_MAX_VAR_DIMS] = {0};
#else
static const size_t coord_zero[NC_MAX_VAR_DIMS];
#endif

dnl
dnl NCPUTVAR(Abbrev, Type)
dnl
define(`NCPUTVAR',dnl
`dnl
int
nc_put_var_$1(int ncid, int varid, const $2 *value)
{
	int status = NC_NOERR;
	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)
	if(varp->ndims == 0) /* scalar variable */
	{
		const size_t zed = 0;
		return( putNCv_$1(ncp, varp, &zed, 1, value) );
	}

	if(!IS_RECVAR(varp))
	{
		return(putNCv_$1(ncp, varp, coord_zero, *varp->dsizes, value));
	}
	/* else */

	if(varp->ndims == 1
			&& ncp->recsize <= varp->len)
	{
		/* one dimensional && the only record variable  */
		return(putNCv_$1(ncp, varp, coord_zero, NC_get_numrecs(ncp),
			value));
	}
	/* else */

	{
	ALLOC_ONSTACK(coord, size_t, varp->ndims);
	size_t elemsPerRec = 1;
	const size_t nrecs = NC_get_numrecs(ncp);
	(void) memset(coord, 0, varp->ndims * sizeof(size_t));
	/* TODO: fix dsizes to avoid this nonsense */
	if(varp->ndims > 1)
		elemsPerRec = varp->dsizes[1];
	while(*coord < nrecs)
	{
		const int lstatus = putNCv_$1(ncp, varp, coord, elemsPerRec,
				 value);
		if(lstatus != NC_NOERR)
		{
			if(lstatus != NC_ERANGE)
			{
				status = lstatus;
				/* fatal for the loop */
				break;
			}
			/* else NC_ERANGE, not fatal for the loop */
			if(status == NC_NOERR)
				status = lstatus;
		}
		value += elemsPerRec;
		(*coord)++;
	}
	FREE_ONSTACK(coord);
	} /* elemsPerRec */

	return status;
}
')dnl

NCPUTVAR(text, char)

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


dnl
dnl NCGETVAR(Abbrv, Type)
dnl
define(`NCGETVAR',dnl
`dnl
int
nc_get_var_$1(int ncid, int varid, $2 *value)
{
	int status = NC_NOERR;
	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 */

	if(varp->ndims == 0) /* scalar variable */
	{
		const size_t zed = 0;
		return( getNCv_$1(ncp, varp, &zed, 1, value) );
	}

NCTEXTCOND($1)

	if(!IS_RECVAR(varp))
	{
		return(getNCv_$1(ncp, varp, coord_zero, *varp->dsizes, value));
	}
	/* else */

	if(varp->ndims == 1
			&& ncp->recsize <= varp->len)
	{
		/* one dimensional && the only record variable  */
		return(getNCv_$1(ncp, varp, coord_zero, NC_get_numrecs(ncp),
			value));
	}
	/* else */

	{
	ALLOC_ONSTACK(coord, size_t, varp->ndims);
	size_t elemsPerRec = 1;
	const size_t nrecs = NC_get_numrecs(ncp);
	(void) memset(coord, 0, varp->ndims * sizeof(size_t));
	/* TODO: fix dsizes to avoid this nonsense */
	if(varp->ndims > 1)
		elemsPerRec = varp->dsizes[1];
	while(*coord < nrecs)
	{
		const int lstatus = getNCv_$1(ncp, varp, coord, elemsPerRec,
				value);
		if(lstatus != NC_NOERR)
		{
			if(lstatus != NC_ERANGE)
			{
				status = lstatus;
				/* fatal for the loop */
				break;
			}
			/* else NC_ERANGE, not fatal for the loop */
			if(status == NC_NOERR)
				status = lstatus;
		}
		value += elemsPerRec;
		(*coord)++;
	}
	FREE_ONSTACK(coord);
	} /* elemsPerRec */

	return status;
}
')dnl

NCGETVAR(text, char)

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


/* Begin putgetg.c */

dnl
dnl  NC_VARM_Upper_Body(void)
dnl
define(`NC_VARM_Upper_Body',dnl
`dnl
	int status = ENOERR;
	NC *ncp;
	NC_var *varp;
	int maxidim;	/* maximum dimensional index */

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

	if (NC_indef (ncp))
	{
		return NC_EINDEFINE;
	}
')dnl
dnl
dnl  NC_VARM_Mid_Body(put_or_get, Abbrv)
dnl
define(`NC_VARM_Mid_Body',dnl
`dnl
	varp = NC_lookupvar (ncp, varid);
	if (varp == NULL)
		return NC_ENOTVAR;

NCTEXTCOND($2)
	maxidim = (int) varp->ndims - 1;

	if (maxidim < 0)
	{
		/*
		 * The variable is a scalar; consequently,
		 * there s only one thing to get and only one place to put it.
		 * (Why was I called?)
		 */
		return $1 (ncp, varp, start, 1, value);
	}
	
	/*
	 * else
	 * The variable is an array.
	 */
	{
		int idim;
		size_t *mystart = NULL;
		size_t *myedges;
		size_t *iocount;	/* count vector */
		size_t *stop;	/* stop indexes */
		size_t *length;	/* edge lengths in bytes */
		ptrdiff_t *mystride;
		ptrdiff_t *mymap;

		/*
		 * Verify stride argument.
		 */
		for (idim = 0; idim <= maxidim; ++idim)
		{
			if (stride != NULL
				&& (stride[idim] == 0
		/* cast needed for braindead systems with signed size_t */
				|| (unsigned long) stride[idim] >= X_INT_MAX))
			{
				return NC_ESTRIDE;
			}
		}

		/* assert(sizeof(ptrdiff_t) >= sizeof(size_t)); */
		mystart = (size_t *)calloc(varp->ndims * 7, sizeof(ptrdiff_t));
		if(mystart == NULL)
			return NC_ENOMEM;
		myedges = mystart + varp->ndims;
		iocount = myedges + varp->ndims;
		stop = iocount + varp->ndims;
		length = stop + varp->ndims;
		mystride = (ptrdiff_t *)(length + varp->ndims);
		mymap = mystride + varp->ndims;

		/*
		 * Initialize I/O parameters.
		 */
		for (idim = maxidim; idim >= 0; --idim)
		{
			mystart[idim] = start != NULL
				? start[idim]
				: 0;

			if (edges[idim] == 0)
			{
				status = NC_NOERR;	/* read/write no data */
				goto done;
			}

			myedges[idim] = edges != NULL
				? edges[idim]
				: idim == 0 && IS_RECVAR (varp)
				? NC_get_numrecs(ncp) - mystart[idim]
				: varp->shape[idim] - mystart[idim];
			mystride[idim] = stride != NULL
				? stride[idim]
				: 1;
			mymap[idim] = map != NULL
				? map[idim]
				: idim == maxidim
				? 1
				: mymap[idim + 1] * (ptrdiff_t) myedges[idim + 1];

			iocount[idim] = 1;
			length[idim] = mymap[idim] * myedges[idim];
			stop[idim] = mystart[idim] + myedges[idim] * mystride[idim];
		}
')dnl
dnl
dnl  NC_VARM_Lower_Body(put_or_get)
dnl
define(`NC_VARM_Lower_Body',dnl
`dnl
		/*
		 * As an optimization, adjust I/O parameters when the fastest 
		 * dimension has unity stride both externally and internally.
		 * In this case, the user could have called a simpler routine
		 * (i.e. ncvar$1()
		 */
		if (mystride[maxidim] == 1
			&& mymap[maxidim] == 1)
		{
			iocount[maxidim] = myedges[maxidim];
			mystride[maxidim] = (ptrdiff_t) myedges[maxidim];
			mymap[maxidim] = (ptrdiff_t) length[maxidim];
		}

		/*
		 * Perform I/O.  Exit when done.
		 */
		for (;;)
		{
			/* TODO: */
			int lstatus = $1 (ncid, varid, mystart, iocount,
						value);
			if (lstatus != NC_NOERR 
				&& (status == NC_NOERR || lstatus != NC_ERANGE))
				status = lstatus;

			/*
			 * The following code permutes through the variable s
			 * external start-index space and it s internal address
			 * space.  At the UPC, this algorithm is commonly
			 * called "odometer code".
			 */
			idim = maxidim;
		carry:
			value += mymap[idim];
			mystart[idim] += mystride[idim];
			if (mystart[idim] == stop[idim])
			{
				mystart[idim] = start[idim];
				value -= length[idim];
				if (--idim < 0)
					break; /* normal return */
				goto carry;
			}
		} /* I/O loop */
	done:
		free(mystart);
	} /* variable is array */
	return status;
')dnl

dnl
dnl NCGETVARS(Abbrv, Type)
dnl
define(`NCGETVARS',dnl
`dnl
int
nc_get_vars_$1 (
	int ncid,
	int varid,
	const size_t * start,
	const size_t * edges,
	const ptrdiff_t * stride,
	$2 *value)
{
	return nc_get_varm_$1 (ncid, varid, start, edges,
			 stride, 0, value);
}
')dnl

NCGETVARS(text, char)

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


dnl
dnl NCPUTVARS(Abbrv, Type)
dnl
define(`NCPUTVARS',dnl
`dnl
int
nc_put_vars_$1 (
	int ncid,
	int varid,
	const size_t * start,
	const size_t * edges,
	const ptrdiff_t * stride,
	const $2 *value)
{
	return nc_put_varm_$1 (ncid, varid, start, edges,
			 stride, 0, value);
}
')dnl

NCPUTVARS(text, char)

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


/*
 * Generalized hyperslab input.
 */
dnl
dnl NCGETVARM(Abbrv, Type)
dnl
define(`NCGETVARM',dnl
`dnl
int
nc_get_varm_$1(int ncid, int varid,
	const size_t *start, const size_t *edges,
	const ptrdiff_t *stride,
	const ptrdiff_t *map,
	$2 *value)
{
NC_VARM_Upper_Body()
NC_VARM_Mid_Body(getNCv_$1, $1)
		/*
		 * Check start, edges
		 */
		for (idim = maxidim; idim >= 0; --idim)
		{
			size_t dimlen =