ncx.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.
 * 	
 * 	This file contains some routines derived from code
 *	which is copyrighted by Sun Microsystems, Inc.
 *	The "#ifdef vax" versions of
 *		 ncx_put_float_float()
 *		 ncx_get_float_float()
 *		 ncx_put_double_double()
 *		 ncx_get_double_double()
 *		 ncx_putn_float_float()
 *		 ncx_getn_float_float()
 *		 ncx_putn_double_double()
 *		 ncx_getn_double_double()
 * 	are derived from xdr_float() and xdr_double() routines
 *	in the freely available, copyrighted Sun RPCSRC 3.9
 *	distribution, xdr_float.c.
 * 	Our "value added" is that these are always memory to memory,
 *	they handle IEEE subnormals properly, and their "n" versions
 *	operate speedily on arrays.
 */
/* $Id: ncx.c 2501 2007-11-20 02:33:29Z benkirk $ */

/*
 * An external data representation interface.
 */

#include "ncx.h"
#include <string.h>
#include <limits.h>
/* alias poorly named limits.h macros */
#define  SHORT_MAX  SHRT_MAX
#define  SHORT_MIN  SHRT_MIN
#define USHORT_MAX USHRT_MAX
#include <float.h>
#ifndef FLT_MAX /* This POSIX macro missing on some systems */
# ifndef NO_IEEE_FLOAT
# define FLT_MAX 3.40282347e+38f
# else
# error "You will need to define FLT_MAX"
# endif
#endif
#include <assert.h>

/*
 * If the machine's float domain is "smaller" than the external one
 * use the machine domain
 */
#if defined(FLT_MAX_EXP) && FLT_MAX_EXP < 128 /* 128 is X_FLT_MAX_EXP */
#undef X_FLOAT_MAX
# define X_FLOAT_MAX FLT_MAX
#undef X_FLOAT_MIN
# define X_FLOAT_MIN (-X_FLOAT_MAX)
#endif

#if _SX /* NEC SUPER UX */
#if _INT64
#undef  INT_MAX /* workaround cpp bug */
#define INT_MAX  X_INT_MAX
#undef  INT_MIN /* workaround cpp bug */
#define INT_MIN  X_INT_MIN
#undef  LONG_MAX /* workaround cpp bug */
#define LONG_MAX  X_INT_MAX
#undef  LONG_MIN /* workaround cpp bug */
#define LONG_MIN  X_INT_MIN
#elif _LONG64
#undef  LONG_MAX /* workaround cpp bug */
#define LONG_MAX  4294967295L
#undef  LONG_MIN /* workaround cpp bug */
#define LONG_MIN -4294967295L
#endif
#endif /* _SX */

static const char nada[X_ALIGN] = {0, 0, 0, 0};

#ifndef WORDS_BIGENDIAN
/* LITTLE_ENDIAN: DEC and intel */
/*
 * Routines to convert to BIGENDIAN.
 * Optimize the swapn?b() and swap?b() routines aggressivly.
 */

#define SWAP2(a) ( (((a) & 0xff) << 8) | \
		(((a) >> 8) & 0xff) )

#define SWAP4(a) ( ((a) << 24) | \
		(((a) <<  8) & 0x00ff0000) | \
		(((a) >>  8) & 0x0000ff00) | \
		(((a) >> 24) & 0x000000ff) )


static void
swapn2b(void *dst, const void *src, size_t nn)
{
	char *op = dst;
	const char *ip = src;

/* unroll the following to reduce loop overhead
 *
 *	while(nn-- != 0)
 *	{
 *		*op++ = *(++ip);
 *		*op++ = *(ip++ -1);
 *	}                                       
 */
	while(nn > 3)
	{
		*op++ = *(++ip);
		*op++ = *(ip++ -1);
		*op++ = *(++ip);
		*op++ = *(ip++ -1);
		*op++ = *(++ip);
		*op++ = *(ip++ -1);
		*op++ = *(++ip);
		*op++ = *(ip++ -1);
		nn -= 4;
	}
	while(nn-- != 0)
	{
		*op++ = *(++ip);
		*op++ = *(ip++ -1);
	}
}

# ifndef vax
static void
swap4b(void *dst, const void *src)
{
	char *op = dst;
	const char *ip = src;
	op[0] = ip[3];
	op[1] = ip[2];
	op[2] = ip[1];
	op[3] = ip[0];
}
# endif /* !vax */

static void
swapn4b(void *dst, const void *src, size_t nn)
{
	char *op = dst;
	const char *ip = src;

/* unroll the following to reduce loop overhead
 *	while(nn-- != 0)
 *	{
 *		op[0] = ip[3];
 *		op[1] = ip[2];
 *		op[2] = ip[1];
 *		op[3] = ip[0];
 *		op += 4;
 *		ip += 4;
 *	}
 */
	while(nn > 3)
	{
		op[0] = ip[3];
		op[1] = ip[2];
		op[2] = ip[1];
		op[3] = ip[0];
		op[4] = ip[7];
		op[5] = ip[6];
		op[6] = ip[5];
		op[7] = ip[4];
		op[8] = ip[11];
		op[9] = ip[10];
		op[10] = ip[9];
		op[11] = ip[8];
		op[12] = ip[15];
		op[13] = ip[14];
		op[14] = ip[13];
		op[15] = ip[12];
		op += 16;
		ip += 16;
		nn -= 4;
	}
	while(nn-- != 0)
	{
		op[0] = ip[3];
		op[1] = ip[2];
		op[2] = ip[1];
		op[3] = ip[0];
		op += 4;
		ip += 4;
	}
}

# ifndef vax
static void
swap8b(void *dst, const void *src)
{
	char *op = dst;
	const char *ip = src;
	op[0] = ip[7];
	op[1] = ip[6];
	op[2] = ip[5];
	op[3] = ip[4];
	op[4] = ip[3];
	op[5] = ip[2];
	op[6] = ip[1];
	op[7] = ip[0];
}
# endif /* !vax */

# ifndef vax
static void
swapn8b(void *dst, const void *src, size_t nn)
{
	char *op = dst;
	const char *ip = src;

/* unroll the following to reduce loop overhead
 *	while(nn-- != 0)
 *	{
 *		op[0] = ip[7];
 *		op[1] = ip[6];
 *		op[2] = ip[5];
 *		op[3] = ip[4];
 *		op[4] = ip[3];
 *		op[5] = ip[2];
 *		op[6] = ip[1];
 *		op[7] = ip[0];
 *		op += 8;
 *		ip += 8;
 *	}
 */
	while(nn > 1)
	{
		op[0] = ip[7];
		op[1] = ip[6];
		op[2] = ip[5];
		op[3] = ip[4];
		op[4] = ip[3];
		op[5] = ip[2];
		op[6] = ip[1];
		op[7] = ip[0];
		op[8] = ip[15];
		op[9] = ip[14];
		op[10] = ip[13];
		op[11] = ip[12];
		op[12] = ip[11];
		op[13] = ip[10];
		op[14] = ip[9];
		op[15] = ip[8];
		op += 16;
		ip += 16;
		nn -= 2;
	}
	while(nn-- != 0)
	{
		op[0] = ip[7];
		op[1] = ip[6];
		op[2] = ip[5];
		op[3] = ip[4];
		op[4] = ip[3];
		op[5] = ip[2];
		op[6] = ip[1];
		op[7] = ip[0];
		op += 8;
		ip += 8;
	}
}
# endif /* !vax */

#endif /* LITTLE_ENDIAN */


/*
 * Primitive numeric conversion functions.
 */

/* x_schar */

 /* We don't implement and x_schar primitives. */


/* x_short */

#if SHORT_MAX == X_SHORT_MAX
typedef short ix_short;
#define SIZEOF_IX_SHORT SIZEOF_SHORT
#define IX_SHORT_MAX SHORT_MAX
#elif INT_MAX >= X_SHORT_MAX
typedef int ix_short;
#define SIZEOF_IX_SHORT SIZEOF_INT
#define IX_SHORT_MAX INT_MAX
#elif LONG_MAX >= X_SHORT_MAX
typedef long ix_short;
#define SIZEOF_IX_SHORT SIZEOF_LONG
#define IX_SHORT_MAX LONG_MAX
#else
#error "ix_short implementation"
#endif

static void
get_ix_short(const void *xp, ix_short *ip)
{
	const uchar *cp = (const uchar *) xp;
	*ip = *cp++ << 8;
#if SIZEOF_IX_SHORT > X_SIZEOF_SHORT
	if(*ip & 0x8000)
	{
		/* extern is negative */
		*ip |= (~(0xffff)); /* N.B. Assumes "twos complement" */
	}
#endif
	*ip |= *cp; 
}

static void
put_ix_short(void *xp, const ix_short *ip)
{
	uchar *cp = (uchar *) xp;
	*cp++ = (*ip) >> 8;
	*cp = (*ip) & 0xff;
}


int
ncx_get_short_schar(const void *xp, schar *ip)
{
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
	if(xx > SCHAR_MAX || xx < SCHAR_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_short_uchar(const void *xp, uchar *ip)
{
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
	if(xx > UCHAR_MAX || xx < 0)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_short_short(const void *xp, short *ip)
{
#if SIZEOF_IX_SHORT == SIZEOF_SHORT && IX_SHORT_MAX == SHORT_MAX
	get_ix_short(xp, (ix_short *)ip);
	return ENOERR;
#else
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
#   if IX_SHORT_MAX > SHORT_MAX
	if(xx > SHORT_MAX || xx < SHORT_MIN)
		return NC_ERANGE;
#   endif
	return ENOERR;
#endif
}

int
ncx_get_short_int(const void *xp, int *ip)
{
#if SIZEOF_IX_SHORT == SIZEOF_INT && IX_SHORT_MAX == INT_MAX
	get_ix_short(xp, (ix_short *)ip);
	return ENOERR;
#else
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
#   if IX_SHORT_MAX > INT_MAX
	if(xx > INT_MAX || xx < INT_MIN)
		return NC_ERANGE;
#   endif
	return ENOERR;
#endif
}

int
ncx_get_short_long(const void *xp, long *ip)
{
#if SIZEOF_IX_SHORT == SIZEOF_LONG && IX_SHORT_MAX == LONG_MAX
	get_ix_short(xp, (ix_short *)ip);
	return ENOERR;
#else
	/* assert(LONG_MAX >= X_SHORT_MAX); */
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
	return ENOERR;
#endif
}

int
ncx_get_short_float(const void *xp, float *ip)
{
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
#if 0	/* TODO: determine when necessary */
	if(xx > FLT_MAX || xx < (-FLT_MAX))
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_get_short_double(const void *xp, double *ip)
{
	/* assert(DBL_MAX >= X_SHORT_MAX); */
	ix_short xx;
	get_ix_short(xp, &xx);
	*ip = xx;
	return ENOERR;
}

int
ncx_put_short_schar(void *xp, const schar *ip)
{
	uchar *cp = (uchar *) xp;
	if(*ip & 0x80)
		*cp++ = 0xff;
	else
		*cp++ = 0;
	*cp = (uchar)*ip;
	return ENOERR;
}

int
ncx_put_short_uchar(void *xp, const uchar *ip)
{
	uchar *cp = (uchar *) xp;
	*cp++ = 0;
	*cp = *ip;
	return ENOERR;
}

int
ncx_put_short_short(void *xp, const short *ip)
{
#if SIZEOF_IX_SHORT == SIZEOF_SHORT && X_SHORT_MAX == SHORT_MAX
	put_ix_short(xp, (const ix_short *)ip);
	return ENOERR;
#else
	ix_short xx = (ix_short)*ip;
	put_ix_short(xp, &xx);
# if X_SHORT_MAX < SHORT_MAX
	if(*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
		return NC_ERANGE;
# endif
	return ENOERR;
#endif
}

int
ncx_put_short_int(void *xp, const int *ip)
{
#if SIZEOF_IX_SHORT == SIZEOF_INT && X_SHORT_MAX == INT_MAX
	put_ix_short(xp, (const ix_short *)ip);
	return ENOERR;
#else
	ix_short xx = (ix_short)*ip;
	put_ix_short(xp, &xx);
# if X_SHORT_MAX < INT_MAX
	if(*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
		return NC_ERANGE;
# endif
	return ENOERR;
#endif
}

int
ncx_put_short_long(void *xp, const long *ip)
{
#if SIZEOF_IX_SHORT == SIZEOF_LONG && X_SHORT_MAX == LONG_MAX
	put_ix_short(xp, (const ix_short *)ip);
	return ENOERR;
#else
	ix_short xx = (ix_short)*ip;
	put_ix_short(xp, &xx);
# if X_SHORT_MAX < LONG_MAX
	if(*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
		return NC_ERANGE;
# endif
	return ENOERR;
#endif
}

int
ncx_put_short_float(void *xp, const float *ip)
{
	ix_short xx = *ip;
	put_ix_short(xp, &xx);
	if(*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_put_short_double(void *xp, const double *ip)
{
	ix_short xx = *ip;
	put_ix_short(xp, &xx);
	if(*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

/* x_int */

#if SHORT_MAX == X_INT_MAX
typedef short ix_int;
#define SIZEOF_IX_INT SIZEOF_SHORT
#define IX_INT_MAX SHORT_MAX
#elif INT_MAX  >= X_INT_MAX
typedef int ix_int;
#define SIZEOF_IX_INT SIZEOF_INT
#define IX_INT_MAX INT_MAX
#elif LONG_MAX  >= X_INT_MAX
typedef long ix_int;
#define SIZEOF_IX_INT SIZEOF_LONG
#define IX_INT_MAX LONG_MAX
#else
#error "ix_int implementation"
#endif


static void
get_ix_int(const void *xp, ix_int *ip)
{
	const uchar *cp = (const uchar *) xp;

	*ip = *cp++ << 24;
#if SIZEOF_IX_INT > X_SIZEOF_INT
	if(*ip & 0x80000000)
	{
		/* extern is negative */
		*ip |= (~(0xffffffff)); /* N.B. Assumes "twos complement" */
	}
#endif
	*ip |= (*cp++ << 16);
	*ip |= (*cp++ << 8);
	*ip |= *cp; 
}

static void
put_ix_int(void *xp, const ix_int *ip)
{
	uchar *cp = (uchar *) xp;

	*cp++ = (*ip) >> 24;
	*cp++ = ((*ip) & 0x00ff0000) >> 16;
	*cp++ = ((*ip) & 0x0000ff00) >>  8;
	*cp   = ((*ip) & 0x000000ff);
}


int
ncx_get_int_schar(const void *xp, schar *ip)
{
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
	if(xx > SCHAR_MAX || xx < SCHAR_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_int_uchar(const void *xp, uchar *ip)
{
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
	if(xx > UCHAR_MAX || xx < 0)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_int_short(const void *xp, short *ip)
{
#if SIZEOF_IX_INT == SIZEOF_SHORT && IX_INT_MAX == SHORT_MAX
	get_ix_int(xp, (ix_int *)ip);
	return ENOERR;
#else
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
#  if IX_INT_MAX > SHORT_MAX
	if(xx > SHORT_MAX || xx < SHORT_MIN)
		return NC_ERANGE;
#  endif
	return ENOERR;
#endif
}

int
ncx_get_int_int(const void *xp, int *ip)
{
#if SIZEOF_IX_INT == SIZEOF_INT && IX_INT_MAX == INT_MAX
	get_ix_int(xp, (ix_int *)ip);
	return ENOERR;
#else
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
#  if IX_INT_MAX > INT_MAX
	if(xx > INT_MAX || xx < INT_MIN)
		return NC_ERANGE;
#  endif
	return ENOERR;
#endif
}

int
ncx_get_int_long(const void *xp, long *ip)
{
#if SIZEOF_IX_INT == SIZEOF_LONG && IX_INT_MAX == LONG_MAX
	get_ix_int(xp, (ix_int *)ip);
	return ENOERR;
#else
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
#  if IX_INT_MAX > LONG_MAX	/* unlikely */
	if(xx > LONG_MAX || xx < LONG_MIN)
		return NC_ERANGE;
#  endif
	return ENOERR;
#endif
}

int
ncx_get_int_float(const void *xp, float *ip)
{
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
#if 0	/* TODO: determine when necessary */
	if(xx > FLT_MAX || xx < (-FLT_MAX))
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_get_int_double(const void *xp, double *ip)
{
	/* assert((DBL_MAX >= X_INT_MAX); */
	ix_int xx;
	get_ix_int(xp, &xx);
	*ip = xx;
	return ENOERR;
}

int
ncx_put_int_schar(void *xp, const schar *ip)
{
	uchar *cp = (uchar *) xp;
	if(*ip & 0x80)
	{
		*cp++ = 0xff;
		*cp++ = 0xff;
		*cp++ = 0xff;
	}
	else
	{
		*cp++ = 0x00;
		*cp++ = 0x00;
		*cp++ = 0x00;
	}
	*cp = (uchar)*ip;
	return ENOERR;
}

int
ncx_put_int_uchar(void *xp, const uchar *ip)
{
	uchar *cp = (uchar *) xp;
	*cp++ = 0x00;
	*cp++ = 0x00;
	*cp++ = 0x00;
	*cp   = *ip;
	return ENOERR;
}

int
ncx_put_int_short(void *xp, const short *ip)
{
#if SIZEOF_IX_INT == SIZEOF_SHORT && IX_INT_MAX == SHORT_MAX
	put_ix_int(xp, (ix_int *)ip);
	return ENOERR;
#else