ncx.c

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

/* What IEEE double precision floating point looks like on a Vax */
struct	ieee_double {
	unsigned int	exp_hi   : 7;
	unsigned int	sign     : 1;
	unsigned int 	mant_6   : 4;
	unsigned int	exp_lo   : 4;
	unsigned int	mant_5   : 8;
	unsigned int	mant_4   : 8;

	unsigned int	mant_lo  : 32;
};

/* Vax double precision floating point */
struct  vax_double {
	unsigned int	mantissa1 : 7;
	unsigned int	exp       : 8;
	unsigned int	sign      : 1;
	unsigned int	mantissa2 : 16;
	unsigned int	mantissa3 : 16;
	unsigned int	mantissa4 : 16;
};

#define VAX_DBL_BIAS	0x81
#define IEEE_DBL_BIAS	0x3ff
#define MASK(nbits)	((1 << nbits) - 1)

static const struct dbl_limits {
	struct	vax_double d;
	struct	ieee_double ieee;
} dbl_limits[2] = {
	{{ 0x7f, 0xff, 0x0, 0xffff, 0xffff, 0xffff },	/* Max Vax */
	{ 0x7f, 0x0, 0x0, 0xf, 0x0, 0x0, 0x0}}, /* Max IEEE */
	{{ 0x0, 0x0, 0x0, 0x0, 0x0, 0x0},		/* Min Vax */
	{ 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0}}, /* Min IEEE */
};


static void
get_ix_double(const void *xp, double *ip)
{
	struct vax_double *const vdp =
			 (struct vax_double *)ip;
	const struct ieee_double *const idp =
			 (const struct ieee_double *) xp;
	{
		const struct dbl_limits *lim;
		int ii;
		for (ii = 0, lim = dbl_limits;
			ii < sizeof(dbl_limits)/sizeof(struct dbl_limits);
			ii++, lim++)
		{
			if ((idp->mant_lo == lim->ieee.mant_lo)
				&& (idp->mant_4 == lim->ieee.mant_4)
				&& (idp->mant_5 == lim->ieee.mant_5)
				&& (idp->mant_6 == lim->ieee.mant_6)
				&& (idp->exp_lo == lim->ieee.exp_lo)
				&& (idp->exp_hi == lim->ieee.exp_hi)
				)
			{
				*vdp = lim->d;
				goto doneit;
			}
		}
	}
	{
		unsigned exp = idp->exp_hi << 4 | idp->exp_lo;
		vdp->exp = exp - IEEE_DBL_BIAS + VAX_DBL_BIAS;
	}
	{
		unsigned mant_hi = ((idp->mant_6 << 16)
				 | (idp->mant_5 << 8)
				 | idp->mant_4);
		unsigned mant_lo = SWAP4(idp->mant_lo);
		vdp->mantissa1 = (mant_hi >> 13);
		vdp->mantissa2 = ((mant_hi & MASK(13)) << 3)
				| (mant_lo >> 29);
		vdp->mantissa3 = (mant_lo >> 13);
		vdp->mantissa4 = (mant_lo << 3);
	}
	doneit:
		vdp->sign = idp->sign;

}


static void
put_ix_double(void *xp, const double *ip)
{
	const struct vax_double *const vdp = 
			(const struct vax_double *)ip;
	struct ieee_double *const idp =
			 (struct ieee_double *) xp;

	if ((vdp->mantissa4 > (dbl_limits[0].d.mantissa4 - 3)) &&
		(vdp->mantissa3 == dbl_limits[0].d.mantissa3) &&
		(vdp->mantissa2 == dbl_limits[0].d.mantissa2) &&
		(vdp->mantissa1 == dbl_limits[0].d.mantissa1) &&
		(vdp->exp == dbl_limits[0].d.exp))
	{
		*idp = dbl_limits[0].ieee;
		goto shipit;
	}
	if ((vdp->mantissa4 == dbl_limits[1].d.mantissa4) &&
		(vdp->mantissa3 == dbl_limits[1].d.mantissa3) &&
		(vdp->mantissa2 == dbl_limits[1].d.mantissa2) &&
		(vdp->mantissa1 == dbl_limits[1].d.mantissa1) &&
		(vdp->exp == dbl_limits[1].d.exp))
	{
		*idp = dbl_limits[1].ieee;
		goto shipit;
	}

	{
		unsigned exp = vdp->exp - VAX_DBL_BIAS + IEEE_DBL_BIAS;

		unsigned mant_lo = ((vdp->mantissa2 & MASK(3)) << 29) |
			(vdp->mantissa3 << 13) |
			((vdp->mantissa4 >> 3) & MASK(13));

		unsigned mant_hi = (vdp->mantissa1 << 13)
				 | (vdp->mantissa2 >> 3);

		if((vdp->mantissa4 & 7) > 4)
		{
			/* round up */
			mant_lo++;
			if(mant_lo == 0)
			{
				mant_hi++;
				if(mant_hi > 0xffffff)
				{
					mant_hi = 0;
					exp++;
				}
			}
		}

		idp->mant_lo = SWAP4(mant_lo);
		idp->mant_6 = mant_hi >> 16;
		idp->mant_5 = (mant_hi & 0xff00) >> 8;
		idp->mant_4 = mant_hi;
		idp->exp_hi = exp >> 4;
		idp->exp_lo = exp;
	}
		
	shipit:
		idp->sign = vdp->sign;

}

	/* vax */
#elif defined(_CRAY) && !defined(__crayx1)

static void
get_ix_double(const void *xp, double *ip)
{
	const ieee_double *idp = (const ieee_double *) xp;
	cray_single *csp = (cray_single *) ip;

	if(idp->exp == 0)
	{
		/* ieee subnormal */
		*ip = (double)idp->mant;
		if(idp->mant != 0)
		{
			csp->exp -= (ieee_double_bias + 51);
		}
	}
	else
	{
		csp->exp  = idp->exp + cs_id_bias + 1;
		csp->mant = idp->mant >> (52 - 48 + 1);
		csp->mant |= (1 << (48 - 1));
	}
	csp->sign = idp->sign;
}

static void
put_ix_double(void *xp, const double *ip)
{
	ieee_double *idp = (ieee_double *) xp;
	const cray_single *csp = (const cray_single *) ip;

	int ieee_exp = csp->exp - cs_id_bias -1;

	idp->sign = csp->sign;

	if(ieee_exp >= 0x7ff)
	{
		/* NC_ERANGE => ieee Inf */
		idp->exp = 0x7ff;
		idp->mant = 0x0;
	}
	else if(ieee_exp > 0)
	{
		/* normal ieee representation */
		idp->exp  = ieee_exp;
		/* assumes cray rep is in normal form */
		assert(csp->mant & 0x800000000000);
		idp->mant = (((csp->mant << 1) &
				0xffffffffffff) << (52 - 48));
	}
	else if(ieee_exp >= (-(52 -48)))
	{
		/* ieee subnormal, left  */
		const int lshift = (52 - 48) + ieee_exp;
		idp->mant = csp->mant << lshift;
		idp->exp  = 0;
	}
	else if(ieee_exp >= -52)
	{
		/* ieee subnormal, right  */
		const int rshift = (- (52 - 48) - ieee_exp);

		idp->mant = csp->mant >> rshift;

#if 0
		if(csp->mant & (1 << (rshift -1)))
		{
			/* round up */
			idp->mant++;
		}
#endif

		idp->exp  = 0;
	}
	else
	{
		/* smaller than ieee can represent */
		idp->exp = 0;
		idp->mant = 0;
	}
}
#elif _SX && _FLOAT2
static void
get_ix_double(const void *xp, double *ip)
{
	const int ncnv = ie3_fl2(xp, ip, 8, 8, 1);
}

static void
put_ix_double(void *xp, const double *ip)
{
	const int ncnv = fl2_ie3(ip, xp, 8, 8, 1);
}
#else
#error "ix_double implementation"
#endif

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

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

int
ncx_get_double_short(const void *xp, short *ip)
{
	double xx;
	get_ix_double(xp, &xx);
	*ip = (short) xx;
	if(xx > SHORT_MAX || xx < SHORT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_double_int(const void *xp, int *ip)
{
	double xx;
	get_ix_double(xp, &xx);
	*ip = (int) xx;
	if(xx > INT_MAX || xx < INT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_double_long(const void *xp, long *ip)
{
	double xx;
	get_ix_double(xp, &xx);
	*ip = (long) xx;
	if(xx > LONG_MAX || xx < LONG_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_double_float(const void *xp, float *ip)
{
	double xx;
	get_ix_double(xp, &xx);
	if(xx > FLT_MAX || xx < (-FLT_MAX))
	{
		*ip = FLT_MAX;
		return NC_ERANGE;
	}
	if(xx < (-FLT_MAX))
	{
		*ip = (-FLT_MAX);
		return NC_ERANGE;
	}
	*ip = (float) xx;
	return ENOERR;
}

int
ncx_get_double_double(const void *xp, double *ip)
{
	/* TODO */
	get_ix_double(xp, ip);
	return ENOERR;
}


int
ncx_put_double_schar(void *xp, const schar *ip)
{
	double xx = (double) *ip;
	put_ix_double(xp, &xx);
	return ENOERR;
}

int
ncx_put_double_uchar(void *xp, const uchar *ip)
{
	double xx = (double) *ip;
	put_ix_double(xp, &xx);
	return ENOERR;
}

int
ncx_put_double_short(void *xp, const short *ip)
{
	double xx = (double) *ip;
	put_ix_double(xp, &xx);
#if 0	/* TODO: figure this out */
	if((double)(*ip) > X_DOUBLE_MAX || (double)(*ip) < X_DOUBLE_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_double_int(void *xp, const int *ip)
{
	double xx = (double) *ip;
	put_ix_double(xp, &xx);
#if 0	/* TODO: figure this out */
	if((double)(*ip) > X_DOUBLE_MAX || (double)(*ip) < X_DOUBLE_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_double_long(void *xp, const long *ip)
{
	double xx = (double) *ip;
	put_ix_double(xp, &xx);
#if 1	/* TODO: figure this out */
	if((double)(*ip) > X_DOUBLE_MAX || (double)(*ip) < X_DOUBLE_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_double_float(void *xp, const float *ip)
{
	double xx = (double) *ip;
	put_ix_double(xp, &xx);
#if 1	/* TODO: figure this out */
	if((double)(*ip) > X_DOUBLE_MAX || (double)(*ip) < X_DOUBLE_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_double_double(void *xp, const double *ip)
{
	put_ix_double(xp, ip);
#ifdef NO_IEEE_FLOAT
	if(*ip > X_DOUBLE_MAX || *ip < X_DOUBLE_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}


/* x_size_t */

#if SIZEOF_SIZE_T < X_SIZEOF_SIZE_T
#error "x_size_t implementation"
/* netcdf requires size_t which can hold a values from 0 to 2^32 -1 */
#endif

int
ncx_put_size_t(void **xpp, const size_t *ulp)
{
	/* similar to put_ix_int() */
	uchar *cp = (uchar *) *xpp;
	assert(*ulp <= X_SIZE_MAX);

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

	*xpp = (void *)((char *)(*xpp) + X_SIZEOF_SIZE_T);
	return ENOERR;
}

int
ncx_get_size_t(const void **xpp,  size_t *ulp)
{
	/* similar to get_ix_int */
	const uchar *cp = (const uchar *) *xpp;

	*ulp = (unsigned)(*cp++ << 24);
	*ulp |= (*cp++ << 16);
	*ulp |= (*cp++ << 8);
	*ulp |= *cp; 

	*xpp = (const void *)((const char *)(*xpp) + X_SIZEOF_SIZE_T);
	return ENOERR;
}

/* x_off_t */

int
ncx_put_off_t(void **xpp, const off_t *lp, size_t sizeof_off_t)
{
	/* similar to put_ix_int() */
	uchar *cp = (uchar *) *xpp;
		/* No negative offsets stored in netcdf */
	if (*lp < 0) {
	  /* Assume this is an overflow of a 32-bit int... */
	  return ERANGE;
	}
	  
	assert(sizeof_off_t == 4 || sizeof_off_t == 8);

	if (sizeof_off_t == 4) {
		*cp++ = (uchar) ((*lp)               >> 24);
		*cp++ = (uchar)(((*lp) & 0x00ff0000) >> 16);
		*cp++ = (uchar)(((*lp) & 0x0000ff00) >>  8);
		*cp   = (uchar)( (*lp) & 0x000000ff);
	} else {
#if SIZEOF_OFF_T == 4
/* Write a 64-bit offset on a system with only a 32-bit offset */
		*cp++ = (uchar)0;
		*cp++ = (uchar)0;
		*cp++ = (uchar)0;
		*cp++ = (uchar)0;

		*cp++ = (uchar)(((*lp) & 0xff000000) >> 24);
		*cp++ = (uchar)(((*lp) & 0x00ff0000) >> 16);
		*cp++ = (uchar)(((*lp) & 0x0000ff00) >>  8);
		*cp   = (uchar)( (*lp) & 0x000000ff);
#else
		*cp++ = (uchar) ((*lp)                          >> 56);
		*cp++ = (uchar)(((*lp) & 0x00ff000000000000ULL) >> 48);
		*cp++ = (uchar)(((*lp) & 0x0000ff0000000000ULL) >> 40);
		*cp++ = (uchar)(((*lp) & 0x000000ff00000000ULL) >> 32);
		*cp++ = (uchar)(((*lp) & 0x00000000ff000000ULL) >> 24);
		*cp++ = (uchar)(((*lp) & 0x0000000000ff0000ULL) >> 16);
		*cp++ = (uchar)(((*lp) & 0x000000000000ff00ULL) >>  8);
		*cp   = (uchar)( (*lp) & 0x00000000000000ffULL);
#endif
	}
	*xpp = (void *)((char *)(*xpp) + sizeof_off_t);
	return ENOERR;
}

int
ncx_get_off_t(const void **xpp, off_t *lp, size_t sizeof_off_t)
{
	/* similar to get_ix_int() */
	const uchar *cp = (const uchar *) *xpp;
	assert(sizeof_off_t == 4 || sizeof_off_t == 8);

 	if (sizeof_off_t == 4) {
		*lp = *cp++ << 24;
		*lp |= (*cp++ << 16);
		*lp |= (*cp++ <<  8);
		*lp |= *cp; 
	} else {
#if SIZEOF_OFF_T == 4
/* Read a 64-bit offset on a system with only a 32-bit offset */
/* If the offset overflows, set an error code and return */
		*lp =  ((off_t)(*cp++) << 24);
		*lp |= ((off_t)(*cp++) << 16);
		*lp |= ((off_t)(*cp++) <<  8);
		*lp |= ((off_t)(*cp++));
/*
 * lp now contains the upper 32-bits of the 64-bit offset.  if lp is
 * not zero, then the dataset is larger than can be represented
 * on this system.  Set an error code and return.
 */
		if (*lp != 0) {
		  return ERANGE;
		}

		*lp  = ((off_t)(*cp++) << 24);
		*lp |= ((off_t)(*cp++) << 16);
		*lp |= ((off_t)(*cp++) <<  8);
		*lp |=  (off_t)*cp;

		if (*lp < 0) {
		  /*
		   * If this fails, then the offset is >2^31, but less
		   * than 2^32 which is not allowed, but is not caught
		   * by the previous check
		   */
		  return ERANGE;
		}
#else
		*lp =  ((off_t)(*cp++) << 56);
		*lp |= ((off_t)(*cp++) << 48);
		*lp |= ((off_t)(*cp++) << 40);
		*lp |= ((off_t)(*cp++) << 32);
		*lp |= ((off_t)(*cp++) << 24);
		*lp |= ((off_t)(*cp++) << 16);
		*lp |= ((off_t)(*cp++) <<  8);
		*lp |=  (off_t)*cp;
#endif
	}
	*xpp = (const void *)((const char *)(*xpp) + sizeof_off_t);
	return ENOERR;
}


/*
 * Aggregate numeric conversion functions.
 */



/* schar */

int
ncx_getn_schar_schar(const void **xpp, size_t nelems, schar *tp)
{
		(void) memcpy(tp, *xpp, nelems);
	*xpp = (void *)((char *)(*xpp) + nelems);
	return ENOERR;

}
int
ncx_getn_schar_uchar(const void **xpp, size_t nelems, uchar *tp)
{
		(void) memcpy(tp, *xpp, nelems);
	*xpp = (void *)((char *)(*xpp) + nelems);
	return ENOERR;

}
int
ncx_getn_schar_short(const void **xpp, size_t nelems, short *tp)
{
	schar *xp = (schar *)(*xpp);

	while(nelems-- != 0)
	{
		*tp++ = *xp++;
	}

	*xpp = (const void *)xp;
	return ENOERR;
}

int
ncx_getn_schar_int(const void **xpp, size_t nelems, int *tp)
{
	schar *xp = (schar *)(*xpp);

	while(nelems-- != 0)
	{
		*tp++ = *xp++;
	}

	*xpp = (const void *)xp;
	return ENOERR;
}

int
ncx_getn_schar_long(const void **xpp, size_t nelems, long *tp)
{
	schar *xp = (schar *)(*xpp);

	while(nelems-- != 0)
	{
		*tp++ = *xp++;
	}

	*xpp = (const void *)xp;
	return ENOERR;
}

int
ncx_getn_schar_float(const void **xpp, size_t nelems, float *tp)
{
	schar *xp = (schar *)(*xpp);

	while(nelems-- != 0)
	{
		*tp++ = *xp++;
	}

	*xpp = (const void *)xp;
	return ENOERR;
}

int
ncx_getn_schar_double(const void **xpp, size_t nelems, double *tp)
{
	schar *xp = (schar *)(*xpp);

	while(nelems-- != 0)
	{
		*tp++ = *xp++;
	}

	*xpp = (const void *)xp;
	return ENOERR;
}


int
ncx_pad_getn_schar_schar(const void **xpp, size_t nelems, schar *tp)
{
		size_t rndup = nelems % X_ALIGN;

	if(rndup)
		rndup = X_ALIGN - rndup;

	(void) memcpy(tp, *xpp, nelems);
	*xpp = (void *)((char *)(*xpp) + nelems + rndup);

	return ENOERR;

}
int
ncx_pad_getn_schar_uchar(const void **xpp, size_t nelems, uchar *tp)
{
		size_t rndup = nelems % X_ALIGN;

	if(rndup)
		rndup = X_ALIGN - rndup;

	(void) memcpy(tp, *xpp, nelems);
	*xpp = (void *)((char *)(*xpp) + nelems + rndup);

	return ENOERR;

}
int
ncx_pad_getn_schar_short(const void **xpp, size_t nelems, short *tp)
{
	size_t rndup = nelems % X_ALIGN;
	schar *xp = (schar *) *xpp;

	if(rndup)
		rndup = X_ALIGN - rndup;

	while(nelems-- != 0)
	{
		*tp++ = *xp++;
	}

	*xpp = (void *)(xp + rndup);
	return ENOERR;
}

int