X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fxdrf%2Flibxdrf.m4;fp=source%2Fwham%2Fsrc-NEWSC%2Fxdrf%2Flibxdrf.m4;h=aecb5b54829074e1234c9a2bf5c4e02380f91568;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC/xdrf/libxdrf.m4 b/source/wham/src-NEWSC/xdrf/libxdrf.m4 new file mode 100644 index 0000000..aecb5b5 --- /dev/null +++ b/source/wham/src-NEWSC/xdrf/libxdrf.m4 @@ -0,0 +1,1233 @@ +/*____________________________________________________________________________ + | + | libxdrf - portable fortran interface to xdr. some xdr routines + | are C routines for compressed coordinates + | + | version 1.1 + | + | This collection of routines is intended to write and read + | data in a portable way to a file, so data written on one type + | of machine can be read back on a different type. + | + | all fortran routines use an integer 'xdrid', which is an id to the + | current xdr file, and is set by xdrfopen. + | most routines have in integer 'ret' which is the return value. + | The value of 'ret' is zero on failure, and most of the time one + | on succes. + | + | There are three routines useful for C users: + | xdropen(), xdrclose(), xdr3dfcoord(). + | The first two replace xdrstdio_create and xdr_destroy, and *must* be + | used when you plan to use xdr3dfcoord(). (they are also a bit + | easier to interface). For writing data other than compressed coordinates + | you should use the standard C xdr routines (see xdr man page) + | + | xdrfopen(xdrid, filename, mode, ret) + | character *(*) filename + | character *(*) mode + | + | this will open the file with the given filename (string) + | and the given mode, it returns an id in xdrid, which is + | to be used in all other calls to xdrf routines. + | mode is 'w' to create, or update an file, for all other + | values of mode the file is opened for reading + | + | you need to call xdrfclose to flush the output and close + | the file. + | Note that you should not use xdrstdio_create, which comes with the + | standard xdr library + | + | xdrfclose(xdrid, ret) + | flush the data to the file, and closes the file; + | You should not use xdr_destroy (which comes standard with + | the xdr libraries. + | + | xdrfbool(xdrid, bp, ret) + | integer pb + | + | This filter produces values of either 1 or 0 + | + | xdrfchar(xdrid, cp, ret) + | character cp + | + | filter that translate between characters and their xdr representation + | Note that the characters in not compressed and occupies 4 bytes. + | + | xdrfdouble(xdrid, dp, ret) + | double dp + | + | read/write a double. + | + | xdrffloat(xdrid, fp, ret) + | float fp + | + | read/write a float. + | + | xdrfint(xdrid, ip, ret) + | integer ip + | + | read/write integer. + | + | xdrflong(xdrid, lp, ret) + | integer lp + | + | this routine has a possible portablility problem due to 64 bits longs. + | + | xdrfshort(xdrid, sp, ret) + | integer *2 sp + | + | xdrfstring(xdrid, sp, maxsize, ret) + | character *(*) + | integer maxsize + | + | read/write a string, with maximum length given by maxsize + | + | xdrfwrapstring(xdris, sp, ret) + | character *(*) + | + | read/write a string (it is the same as xdrfstring accept that it finds + | the stringlength itself. + | + | xdrfvector(xdrid, cp, size, xdrfproc, ret) + | character *(*) + | integer size + | external xdrfproc + | + | read/write an array pointed to by cp, with number of elements + | defined by 'size'. the routine 'xdrfproc' is the name + | of one of the above routines to read/write data (like xdrfdouble) + | In contrast with the c-version you don't need to specify the + | byte size of an element. + | xdrfstring is not allowed here (it is in the c version) + | + | xdrf3dfcoord(xdrid, fp, size, precision, ret) + | real (*) fp + | real precision + | integer size + | + | this is *NOT* a standard xdr routine. I named it this way, because + | it invites people to use the other xdr routines. + | It is introduced to store specifically 3d coordinates of molecules + | (as found in molecular dynamics) and it writes it in a compressed way. + | It starts by multiplying all numbers by precision and + | rounding the result to integer. effectively converting + | all floating point numbers to fixed point. + | it uses an algorithm for compression that is optimized for + | molecular data, but could be used for other 3d coordinates + | as well. There is subtantial overhead involved, so call this + | routine only if you have a large number of coordinates to read/write + | + | ________________________________________________________________________ + | + | Below are the routines to be used by C programmers. Use the 'normal' + | xdr routines to write integers, floats, etc (see man xdr) + | + | int xdropen(XDR *xdrs, const char *filename, const char *type) + | This will open the file with the given filename and the + | given mode. You should pass it an allocated XDR struct + | in xdrs, to be used in all other calls to xdr routines. + | Mode is 'w' to create, or update an file, and for all + | other values of mode the file is opened for reading. + | You need to call xdrclose to flush the output and close + | the file. + | + | Note that you should not use xdrstdio_create, which + | comes with the standard xdr library. + | + | int xdrclose(XDR *xdrs) + | Flush the data to the file, and close the file; + | You should not use xdr_destroy (which comes standard + | with the xdr libraries). + | + | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) + | This is \fInot\fR a standard xdr routine. I named it this + | way, because it invites people to use the other xdr + | routines. + | + | (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl +*/ + + +#include +#include +#include +#include +#include +#include +#include +#include "xdrf.h" + +int ftocstr(char *, int, char *, int); +int ctofstr(char *, int, char *); + +#define MAXID 20 +static FILE *xdrfiles[MAXID]; +static XDR *xdridptr[MAXID]; +static char xdrmodes[MAXID]; +static unsigned int cnt; + +typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *); + +void +FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret') +int *xdrid, *ret; +int *pb; +{ + *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb); + cnt += sizeof(int); +} + +void +FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret') +int *xdrid, *ret; +char *cp; +{ + *ret = xdr_char(xdridptr[*xdrid], cp); + cnt += sizeof(char); +} + +void +FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret') +int *xdrid, *ret; +double *dp; +{ + *ret = xdr_double(xdridptr[*xdrid], dp); + cnt += sizeof(double); +} + +void +FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret') +int *xdrid, *ret; +float *fp; +{ + *ret = xdr_float(xdridptr[*xdrid], fp); + cnt += sizeof(float); +} + +void +FUNCTION(xdrfint) ARGS(`xdrid, ip, ret') +int *xdrid, *ret; +int *ip; +{ + *ret = xdr_int(xdridptr[*xdrid], ip); + cnt += sizeof(int); +} + +void +FUNCTION(xdrflong) ARGS(`xdrid, lp, ret') +int *xdrid, *ret; +long *lp; +{ + *ret = xdr_long(xdridptr[*xdrid], lp); + cnt += sizeof(long); +} + +void +FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret') +int *xdrid, *ret; +short *sp; +{ + *ret = xdr_short(xdridptr[*xdrid], sp); + cnt += sizeof(sp); +} + +void +FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret') +int *xdrid, *ret; +char *ucp; +{ + *ret = xdr_u_char(xdridptr[*xdrid], ucp); + cnt += sizeof(char); +} + +void +FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret') +int *xdrid, *ret; +unsigned long *ulp; +{ + *ret = xdr_u_long(xdridptr[*xdrid], ulp); + cnt += sizeof(unsigned long); +} + +void +FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret') +int *xdrid, *ret; +unsigned short *usp; +{ + *ret = xdr_u_short(xdridptr[*xdrid], usp); + cnt += sizeof(unsigned short); +} + +void +FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret') +int *xdrid, *ret; +float *fp; +int *size; +float *precision; +{ + *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision); +} + +void +FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret') +int *xdrid, *ret; +STRING_ARG_DECL(sp); +int *maxsize; +{ + char *tsp; + + tsp = (char*) malloc(((STRING_LEN(sp)) + 1) * sizeof(char)); + if (tsp == NULL) { + *ret = -1; + return; + } + if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) { + *ret = -1; + free(tsp); + return; + } + *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize); + ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp); + cnt += *maxsize; + free(tsp); +} + +void +FUNCTION(xdrfwrapstring) ARGS(`xdrid, STRING_ARG(sp), ret') +int *xdrid, *ret; +STRING_ARG_DECL(sp); +{ + char *tsp; + int maxsize; + maxsize = (STRING_LEN(sp)) + 1; + tsp = (char*) malloc(maxsize * sizeof(char)); + if (tsp == NULL) { + *ret = -1; + return; + } + if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) { + *ret = -1; + free(tsp); + return; + } + *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize); + ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp); + cnt += maxsize; + free(tsp); +} + +void +FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret') +int *xdrid, *ret; +caddr_t *cp; +int *ccnt; +{ + *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt); + cnt += *ccnt; +} + +void +FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret') +int *xdrid, *ret; +int *pos; +{ + *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos); +} + +void +FUNCTION(xdrf) ARGS(`xdrid, pos') +int *xdrid, *pos; +{ + *pos = xdr_getpos(xdridptr[*xdrid]); +} + +void +FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret') +int *xdrid, *ret; +char *cp; +int *size; +FUNCTION(xdrfproc) elproc; +{ + int lcnt; + cnt = 0; + for (lcnt = 0; lcnt < *size; lcnt++) { + elproc(xdrid, (cp+cnt) , ret); + } +} + + +void +FUNCTION(xdrfclose) ARGS(`xdrid, ret') +int *xdrid; +int *ret; +{ + *ret = xdrclose(xdridptr[*xdrid]); + cnt = 0; +} + +void +FUNCTION(xdrfopen) ARGS(`xdrid, STRING_ARG(fp), STRING_ARG(mode), ret') +int *xdrid; +STRING_ARG_DECL(fp); +STRING_ARG_DECL(mode); +int *ret; +{ + char fname[512]; + char fmode[3]; + + if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) { + *ret = 0; + } + if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode), + STRING_LEN(mode))) { + *ret = 0; + } + + *xdrid = xdropen(NULL, fname, fmode); + if (*xdrid == 0) + *ret = 0; + else + *ret = 1; +} + +/*___________________________________________________________________________ + | + | what follows are the C routines for opening, closing xdr streams + | and the routine to read/write compressed coordinates together + | with some routines to assist in this task (those are marked + | static and cannot be called from user programs) +*/ +#define MAXABS INT_MAX-2 + +#ifndef MIN +#define MIN(x,y) ((x) < (y) ? (x):(y)) +#endif +#ifndef MAX +#define MAX(x,y) ((x) > (y) ? (x):(y)) +#endif +#ifndef SQR +#define SQR(x) ((x)*(x)) +#endif +static int magicints[] = { + 0, 0, 0, 0, 0, 0, 0, 0, 0, + 8, 10, 12, 16, 20, 25, 32, 40, 50, 64, + 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, + 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501, + 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536, + 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561, + 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042, + 8388607, 10568983, 13316085, 16777216 }; + +#define FIRSTIDX 9 +/* note that magicints[FIRSTIDX-1] == 0 */ +#define LASTIDX (sizeof(magicints) / sizeof(*magicints)) + + +/*__________________________________________________________________________ + | + | xdropen - open xdr file + | + | This versions differs from xdrstdio_create, because I need to know + | the state of the file (read or write) so I can use xdr3dfcoord + | in eigther read or write mode, and the file descriptor + | so I can close the file (something xdr_destroy doesn't do). + | +*/ + +int xdropen(XDR *xdrs, const char *filename, const char *type) { + static int init_done = 0; + enum xdr_op lmode; + const char *type1; + int xdrid; + + if (init_done == 0) { + for (xdrid = 1; xdrid < MAXID; xdrid++) { + xdridptr[xdrid] = NULL; + } + init_done = 1; + } + xdrid = 1; + while (xdrid < MAXID && xdridptr[xdrid] != NULL) { + xdrid++; + } + if (xdrid == MAXID) { + return 0; + } + if (*type == 'w' || *type == 'W') { + type = "w+"; + type1 = "a+"; + lmode = XDR_ENCODE; + } else { + type = "r"; + type1 = "r"; + lmode = XDR_DECODE; + } + xdrfiles[xdrid] = fopen(filename, type1); + if (xdrfiles[xdrid] == NULL) { + xdrs = NULL; + return 0; + } + xdrmodes[xdrid] = *type; + /* next test isn't usefull in the case of C language + * but is used for the Fortran interface + * (C users are expected to pass the address of an already allocated + * XDR staructure) + */ + if (xdrs == NULL) { + xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR)); + xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode); + } else { + xdridptr[xdrid] = xdrs; + xdrstdio_create(xdrs, xdrfiles[xdrid], lmode); + } + return xdrid; +} + +/*_________________________________________________________________________ + | + | xdrclose - close a xdr file + | + | This will flush the xdr buffers, and destroy the xdr stream. + | It also closes the associated file descriptor (this is *not* + | done by xdr_destroy). + | +*/ + +int xdrclose(XDR *xdrs) { + int xdrid; + + if (xdrs == NULL) { + fprintf(stderr, "xdrclose: passed a NULL pointer\n"); + exit(1); + } + for (xdrid = 1; xdrid < MAXID; xdrid++) { + if (xdridptr[xdrid] == xdrs) { + + xdr_destroy(xdrs); + fclose(xdrfiles[xdrid]); + xdridptr[xdrid] = NULL; + return 1; + } + } + fprintf(stderr, "xdrclose: no such open xdr file\n"); + exit(1); + +} + +/*____________________________________________________________________________ + | + | sendbits - encode num into buf using the specified number of bits + | + | This routines appends the value of num to the bits already present in + | the array buf. You need to give it the number of bits to use and you + | better make sure that this number of bits is enough to hold the value + | Also num must be positive. + | +*/ + +static void sendbits(int buf[], int num_of_bits, int num) { + + unsigned int cnt, lastbyte; + int lastbits; + unsigned char * cbuf; + + cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); + cnt = (unsigned int) buf[0]; + lastbits = buf[1]; + lastbyte =(unsigned int) buf[2]; + while (num_of_bits >= 8) { + lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/); + cbuf[cnt++] = lastbyte >> lastbits; + num_of_bits -= 8; + } + if (num_of_bits > 0) { + lastbyte = (lastbyte << num_of_bits) | num; + lastbits += num_of_bits; + if (lastbits >= 8) { + lastbits -= 8; + cbuf[cnt++] = lastbyte >> lastbits; + } + } + buf[0] = cnt; + buf[1] = lastbits; + buf[2] = lastbyte; + if (lastbits>0) { + cbuf[cnt] = lastbyte << (8 - lastbits); + } +} + +/*_________________________________________________________________________ + | + | sizeofint - calculate bitsize of an integer + | + | return the number of bits needed to store an integer with given max size + | +*/ + +static int sizeofint(const int size) { + unsigned int num = 1; + int num_of_bits = 0; + + while (size >= num && num_of_bits < 32) { + num_of_bits++; + num <<= 1; + } + return num_of_bits; +} + +/*___________________________________________________________________________ + | + | sizeofints - calculate 'bitsize' of compressed ints + | + | given the number of small unsigned integers and the maximum value + | return the number of bits needed to read or write them with the + | routines receiveints and sendints. You need this parameter when + | calling these routines. Note that for many calls I can use + | the variable 'smallidx' which is exactly the number of bits, and + | So I don't need to call 'sizeofints for those calls. +*/ + +static int sizeofints( const int num_of_ints, unsigned int sizes[]) { + int i, num; + unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp; + num_of_bytes = 1; + bytes[0] = 1; + num_of_bits = 0; + for (i=0; i < num_of_ints; i++) { + tmp = 0; + for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) { + tmp = bytes[bytecnt] * sizes[i] + tmp; + bytes[bytecnt] = tmp & 0xff; + tmp >>= 8; + } + while (tmp != 0) { + bytes[bytecnt++] = tmp & 0xff; + tmp >>= 8; + } + num_of_bytes = bytecnt; + } + num = 1; + num_of_bytes--; + while (bytes[num_of_bytes] >= num) { + num_of_bits++; + num *= 2; + } + return num_of_bits + num_of_bytes * 8; + +} + +/*____________________________________________________________________________ + | + | sendints - send a small set of small integers in compressed format + | + | this routine is used internally by xdr3dfcoord, to send a set of + | small integers to the buffer. + | Multiplication with fixed (specified maximum ) sizes is used to get + | to one big, multibyte integer. Allthough the routine could be + | modified to handle sizes bigger than 16777216, or more than just + | a few integers, this is not done, because the gain in compression + | isn't worth the effort. Note that overflowing the multiplication + | or the byte buffer (32 bytes) is unchecked and causes bad results. + | + */ + +static void sendints(int buf[], const int num_of_ints, const int num_of_bits, + unsigned int sizes[], unsigned int nums[]) { + + int i; + unsigned int bytes[32], num_of_bytes, bytecnt, tmp; + + tmp = nums[0]; + num_of_bytes = 0; + do { + bytes[num_of_bytes++] = tmp & 0xff; + tmp >>= 8; + } while (tmp != 0); + + for (i = 1; i < num_of_ints; i++) { + if (nums[i] >= sizes[i]) { + fprintf(stderr,"major breakdown in sendints num %d doesn't " + "match size %d\n", nums[i], sizes[i]); + exit(1); + } + /* use one step multiply */ + tmp = nums[i]; + for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) { + tmp = bytes[bytecnt] * sizes[i] + tmp; + bytes[bytecnt] = tmp & 0xff; + tmp >>= 8; + } + while (tmp != 0) { + bytes[bytecnt++] = tmp & 0xff; + tmp >>= 8; + } + num_of_bytes = bytecnt; + } + if (num_of_bits >= num_of_bytes * 8) { + for (i = 0; i < num_of_bytes; i++) { + sendbits(buf, 8, bytes[i]); + } + sendbits(buf, num_of_bits - num_of_bytes * 8, 0); + } else { + for (i = 0; i < num_of_bytes-1; i++) { + sendbits(buf, 8, bytes[i]); + } + sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]); + } +} + + +/*___________________________________________________________________________ + | + | receivebits - decode number from buf using specified number of bits + | + | extract the number of bits from the array buf and construct an integer + | from it. Return that value. + | +*/ + +static int receivebits(int buf[], int num_of_bits) { + + int cnt, num; + unsigned int lastbits, lastbyte; + unsigned char * cbuf; + int mask = (1 << num_of_bits) -1; + + cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf); + cnt = buf[0]; + lastbits = (unsigned int) buf[1]; + lastbyte = (unsigned int) buf[2]; + + num = 0; + while (num_of_bits >= 8) { + lastbyte = ( lastbyte << 8 ) | cbuf[cnt++]; + num |= (lastbyte >> lastbits) << (num_of_bits - 8); + num_of_bits -=8; + } + if (num_of_bits > 0) { + if (lastbits < num_of_bits) { + lastbits += 8; + lastbyte = (lastbyte << 8) | cbuf[cnt++]; + } + lastbits -= num_of_bits; + num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1); + } + num &= mask; + buf[0] = cnt; + buf[1] = lastbits; + buf[2] = lastbyte; + return num; +} + +/*____________________________________________________________________________ + | + | receiveints - decode 'small' integers from the buf array + | + | this routine is the inverse from sendints() and decodes the small integers + | written to buf by calculating the remainder and doing divisions with + | the given sizes[]. You need to specify the total number of bits to be + | used from buf in num_of_bits. + | +*/ + +static void receiveints(int buf[], const int num_of_ints, int num_of_bits, + unsigned int sizes[], int nums[]) { + int bytes[32]; + int i, j, num_of_bytes, p, num; + + bytes[1] = bytes[2] = bytes[3] = 0; + num_of_bytes = 0; + while (num_of_bits > 8) { + bytes[num_of_bytes++] = receivebits(buf, 8); + num_of_bits -= 8; + } + if (num_of_bits > 0) { + bytes[num_of_bytes++] = receivebits(buf, num_of_bits); + } + for (i = num_of_ints-1; i > 0; i--) { + num = 0; + for (j = num_of_bytes-1; j >=0; j--) { + num = (num << 8) | bytes[j]; + p = num / sizes[i]; + bytes[j] = p; + num = num - p * sizes[i]; + } + nums[i] = num; + } + nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24); +} + +/*____________________________________________________________________________ + | + | xdr3dfcoord - read or write compressed 3d coordinates to xdr file. + | + | this routine reads or writes (depending on how you opened the file with + | xdropen() ) a large number of 3d coordinates (stored in *fp). + | The number of coordinates triplets to write is given by *size. On + | read this number may be zero, in which case it reads as many as were written + | or it may specify the number if triplets to read (which should match the + | number written). + | Compression is achieved by first converting all floating numbers to integer + | using multiplication by *precision and rounding to the nearest integer. + | Then the minimum and maximum value are calculated to determine the range. + | The limited range of integers so found, is used to compress the coordinates. + | In addition the differences between succesive coordinates is calculated. + | If the difference happens to be 'small' then only the difference is saved, + | compressing the data even more. The notion of 'small' is changed dynamically + | and is enlarged or reduced whenever needed or possible. + | Extra compression is achieved in the case of GROMOS and coordinates of + | water molecules. GROMOS first writes out the Oxygen position, followed by + | the two hydrogens. In order to make the differences smaller (and thereby + | compression the data better) the order is changed into first one hydrogen + | then the oxygen, followed by the other hydrogen. This is rather special, but + | it shouldn't harm in the general case. + | + */ + +int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) { + + + static int *ip = NULL; + static int oldsize; + static int *buf; + + int minint[3], maxint[3], mindiff, *lip, diff; + int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx; + int minidx, maxidx; + unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip; + int flag, k; + int small, smaller, larger, i, is_small, is_smaller, run, prevrun; + float *lfp, lf; + int tmp, *thiscoord, prevcoord[3]; + unsigned int tmpcoord[30]; + + int bufsize, xdrid, lsize; + unsigned int bitsize; + float inv_precision; + int errval = 1; + + /* find out if xdrs is opened for reading or for writing */ + xdrid = 0; + while (xdridptr[xdrid] != xdrs) { + xdrid++; + if (xdrid >= MAXID) { + fprintf(stderr, "xdr error. no open xdr stream\n"); + exit (1); + } + } + if (xdrmodes[xdrid] == 'w') { + + /* xdrs is open for writing */ + + if (xdr_int(xdrs, size) == 0) + return 0; + size3 = *size * 3; + /* when the number of coordinates is small, don't try to compress; just + * write them as floats using xdr_vector + */ + if (*size <= 9 ) { + return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp), + (xdrproc_t)xdr_float)); + } + + xdr_float(xdrs, precision); + if (ip == NULL) { + ip = (int *)malloc(size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = size3 * 1.2; + buf = (int *)malloc(bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } else if (*size > oldsize) { + ip = (int *)realloc(ip, size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = size3 * 1.2; + buf = (int *)realloc(buf, bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } + /* buf[0-2] are special and do not contain actual data */ + buf[0] = buf[1] = buf[2] = 0; + minint[0] = minint[1] = minint[2] = INT_MAX; + maxint[0] = maxint[1] = maxint[2] = INT_MIN; + prevrun = -1; + lfp = fp; + lip = ip; + mindiff = INT_MAX; + oldlint1 = oldlint2 = oldlint3 = 0; + while(lfp < fp + size3 ) { + /* find nearest integer */ + if (*lfp >= 0.0) + lf = *lfp * *precision + 0.5; + else + lf = *lfp * *precision - 0.5; + if (fabs(lf) > MAXABS) { + /* scaling would cause overflow */ + errval = 0; + } + lint1 = lf; + if (lint1 < minint[0]) minint[0] = lint1; + if (lint1 > maxint[0]) maxint[0] = lint1; + *lip++ = lint1; + lfp++; + if (*lfp >= 0.0) + lf = *lfp * *precision + 0.5; + else + lf = *lfp * *precision - 0.5; + if (fabs(lf) > MAXABS) { + /* scaling would cause overflow */ + errval = 0; + } + lint2 = lf; + if (lint2 < minint[1]) minint[1] = lint2; + if (lint2 > maxint[1]) maxint[1] = lint2; + *lip++ = lint2; + lfp++; + if (*lfp >= 0.0) + lf = *lfp * *precision + 0.5; + else + lf = *lfp * *precision - 0.5; + if (fabs(lf) > MAXABS) { + /* scaling would cause overflow */ + errval = 0; + } + lint3 = lf; + if (lint3 < minint[2]) minint[2] = lint3; + if (lint3 > maxint[2]) maxint[2] = lint3; + *lip++ = lint3; + lfp++; + diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3); + if (diff < mindiff && lfp > fp + 3) + mindiff = diff; + oldlint1 = lint1; + oldlint2 = lint2; + oldlint3 = lint3; + } + xdr_int(xdrs, &(minint[0])); + xdr_int(xdrs, &(minint[1])); + xdr_int(xdrs, &(minint[2])); + + xdr_int(xdrs, &(maxint[0])); + xdr_int(xdrs, &(maxint[1])); + xdr_int(xdrs, &(maxint[2])); + + if ((float)maxint[0] - (float)minint[0] >= MAXABS || + (float)maxint[1] - (float)minint[1] >= MAXABS || + (float)maxint[2] - (float)minint[2] >= MAXABS) { + /* turning value in unsigned by subtracting minint + * would cause overflow + */ + errval = 0; + } + sizeint[0] = maxint[0] - minint[0]+1; + sizeint[1] = maxint[1] - minint[1]+1; + sizeint[2] = maxint[2] - minint[2]+1; + + /* check if one of the sizes is to big to be multiplied */ + if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) { + bitsizeint[0] = sizeofint(sizeint[0]); + bitsizeint[1] = sizeofint(sizeint[1]); + bitsizeint[2] = sizeofint(sizeint[2]); + bitsize = 0; /* flag the use of large sizes */ + } else { + bitsize = sizeofints(3, sizeint); + } + lip = ip; + luip = (unsigned int *) ip; + smallidx = FIRSTIDX; + while (smallidx < LASTIDX && magicints[smallidx] < mindiff) { + smallidx++; + } + xdr_int(xdrs, &smallidx); + maxidx = MIN(LASTIDX, smallidx + 8) ; + minidx = maxidx - 8; /* often this equal smallidx */ + smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2; + small = magicints[smallidx] / 2; + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + larger = magicints[maxidx] / 2; + i = 0; + while (i < *size) { + is_small = 0; + thiscoord = (int *)(luip) + i * 3; + if (smallidx < maxidx && i >= 1 && + abs(thiscoord[0] - prevcoord[0]) < larger && + abs(thiscoord[1] - prevcoord[1]) < larger && + abs(thiscoord[2] - prevcoord[2]) < larger) { + is_smaller = 1; + } else if (smallidx > minidx) { + is_smaller = -1; + } else { + is_smaller = 0; + } + if (i + 1 < *size) { + if (abs(thiscoord[0] - thiscoord[3]) < small && + abs(thiscoord[1] - thiscoord[4]) < small && + abs(thiscoord[2] - thiscoord[5]) < small) { + /* interchange first with second atom for better + * compression of water molecules + */ + tmp = thiscoord[0]; thiscoord[0] = thiscoord[3]; + thiscoord[3] = tmp; + tmp = thiscoord[1]; thiscoord[1] = thiscoord[4]; + thiscoord[4] = tmp; + tmp = thiscoord[2]; thiscoord[2] = thiscoord[5]; + thiscoord[5] = tmp; + is_small = 1; + } + + } + tmpcoord[0] = thiscoord[0] - minint[0]; + tmpcoord[1] = thiscoord[1] - minint[1]; + tmpcoord[2] = thiscoord[2] - minint[2]; + if (bitsize == 0) { + sendbits(buf, bitsizeint[0], tmpcoord[0]); + sendbits(buf, bitsizeint[1], tmpcoord[1]); + sendbits(buf, bitsizeint[2], tmpcoord[2]); + } else { + sendints(buf, 3, bitsize, sizeint, tmpcoord); + } + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + thiscoord = thiscoord + 3; + i++; + + run = 0; + if (is_small == 0 && is_smaller == -1) + is_smaller = 0; + while (is_small && run < 8*3) { + if (is_smaller == -1 && ( + SQR(thiscoord[0] - prevcoord[0]) + + SQR(thiscoord[1] - prevcoord[1]) + + SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) { + is_smaller = 0; + } + + tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small; + tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small; + tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small; + + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + + i++; + thiscoord = thiscoord + 3; + is_small = 0; + if (i < *size && + abs(thiscoord[0] - prevcoord[0]) < small && + abs(thiscoord[1] - prevcoord[1]) < small && + abs(thiscoord[2] - prevcoord[2]) < small) { + is_small = 1; + } + } + if (run != prevrun || is_smaller != 0) { + prevrun = run; + sendbits(buf, 1, 1); /* flag the change in run-length */ + sendbits(buf, 5, run+is_smaller+1); + } else { + sendbits(buf, 1, 0); /* flag the fact that runlength did not change */ + } + for (k=0; k < run; k+=3) { + sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]); + } + if (is_smaller != 0) { + smallidx += is_smaller; + if (is_smaller < 0) { + small = smaller; + smaller = magicints[smallidx-1] / 2; + } else { + smaller = small; + small = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx]; + } + } + if (buf[1] != 0) buf[0]++;; + xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */ + return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0])); + } else { + + /* xdrs is open for reading */ + + if (xdr_int(xdrs, &lsize) == 0) + return 0; + if (*size != 0 && lsize != *size) { + fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; " + "%d arg vs %d in file", *size, lsize); + } + *size = lsize; + size3 = *size * 3; + if (*size <= 9) { + return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp), + (xdrproc_t)xdr_float)); + } + xdr_float(xdrs, precision); + if (ip == NULL) { + ip = (int *)malloc(size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = size3 * 1.2; + buf = (int *)malloc(bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } else if (*size > oldsize) { + ip = (int *)realloc(ip, size3 * sizeof(*ip)); + if (ip == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + bufsize = size3 * 1.2; + buf = (int *)realloc(buf, bufsize * sizeof(*buf)); + if (buf == NULL) { + fprintf(stderr,"malloc failed\n"); + exit(1); + } + oldsize = *size; + } + buf[0] = buf[1] = buf[2] = 0; + + xdr_int(xdrs, &(minint[0])); + xdr_int(xdrs, &(minint[1])); + xdr_int(xdrs, &(minint[2])); + + xdr_int(xdrs, &(maxint[0])); + xdr_int(xdrs, &(maxint[1])); + xdr_int(xdrs, &(maxint[2])); + + sizeint[0] = maxint[0] - minint[0]+1; + sizeint[1] = maxint[1] - minint[1]+1; + sizeint[2] = maxint[2] - minint[2]+1; + + /* check if one of the sizes is to big to be multiplied */ + if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) { + bitsizeint[0] = sizeofint(sizeint[0]); + bitsizeint[1] = sizeofint(sizeint[1]); + bitsizeint[2] = sizeofint(sizeint[2]); + bitsize = 0; /* flag the use of large sizes */ + } else { + bitsize = sizeofints(3, sizeint); + } + + xdr_int(xdrs, &smallidx); + maxidx = MIN(LASTIDX, smallidx + 8) ; + minidx = maxidx - 8; /* often this equal smallidx */ + smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2; + small = magicints[smallidx] / 2; + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + larger = magicints[maxidx]; + + /* buf[0] holds the length in bytes */ + + if (xdr_int(xdrs, &(buf[0])) == 0) + return 0; + if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0) + return 0; + buf[0] = buf[1] = buf[2] = 0; + + lfp = fp; + inv_precision = 1.0 / * precision; + run = 0; + i = 0; + lip = ip; + while ( i < lsize ) { + thiscoord = (int *)(lip) + i * 3; + + if (bitsize == 0) { + thiscoord[0] = receivebits(buf, bitsizeint[0]); + thiscoord[1] = receivebits(buf, bitsizeint[1]); + thiscoord[2] = receivebits(buf, bitsizeint[2]); + } else { + receiveints(buf, 3, bitsize, sizeint, thiscoord); + } + + i++; + thiscoord[0] += minint[0]; + thiscoord[1] += minint[1]; + thiscoord[2] += minint[2]; + + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + + + flag = receivebits(buf, 1); + is_smaller = 0; + if (flag == 1) { + run = receivebits(buf, 5); + is_smaller = run % 3; + run -= is_smaller; + is_smaller--; + } + if (run > 0) { + thiscoord += 3; + for (k = 0; k < run; k+=3) { + receiveints(buf, 3, smallidx, sizesmall, thiscoord); + i++; + thiscoord[0] += prevcoord[0] - small; + thiscoord[1] += prevcoord[1] - small; + thiscoord[2] += prevcoord[2] - small; + if (k == 0) { + /* interchange first with second atom for better + * compression of water molecules + */ + tmp = thiscoord[0]; thiscoord[0] = prevcoord[0]; + prevcoord[0] = tmp; + tmp = thiscoord[1]; thiscoord[1] = prevcoord[1]; + prevcoord[1] = tmp; + tmp = thiscoord[2]; thiscoord[2] = prevcoord[2]; + prevcoord[2] = tmp; + *lfp++ = prevcoord[0] * inv_precision; + *lfp++ = prevcoord[1] * inv_precision; + *lfp++ = prevcoord[2] * inv_precision; + } else { + prevcoord[0] = thiscoord[0]; + prevcoord[1] = thiscoord[1]; + prevcoord[2] = thiscoord[2]; + } + *lfp++ = thiscoord[0] * inv_precision; + *lfp++ = thiscoord[1] * inv_precision; + *lfp++ = thiscoord[2] * inv_precision; + } + } else { + *lfp++ = thiscoord[0] * inv_precision; + *lfp++ = thiscoord[1] * inv_precision; + *lfp++ = thiscoord[2] * inv_precision; + } + smallidx += is_smaller; + if (is_smaller < 0) { + small = smaller; + if (smallidx > FIRSTIDX) { + smaller = magicints[smallidx - 1] /2; + } else { + smaller = 0; + } + } else if (is_smaller > 0) { + smaller = small; + small = magicints[smallidx] / 2; + } + sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ; + } + } + return 1; +} + + +