1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
5 * Copyright (c) Erik Lindahl, David van der Spoel 2003,2004.
6 * Coordinate compression (c) by Frans van Hoesel.
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public License
10 * as published by the Free Software Foundation; either version 3
11 * of the License, or (at your option) any later version.
14 /* Get HAVE_RPC_XDR_H, F77_FUNC from config.h if available */
25 #define _FILE_OFFSET_BITS 64
27 /* get fixed-width types if we are using ANSI C99 */
30 #elif (defined HAVE_INTTYPES_H)
31 # include <inttypes.h>
42 /* Default FORTRAN name mangling is: lower case name, append underscore */
44 #define F77_FUNC(name,NAME) name ## _
47 char *exdr_message
[exdrNR
] = {
55 "Compressed 3D coordinate",
64 * Declare our own XDR routines statically if no libraries are present.
65 * Actual implementation is at the end of this file.
67 * We don't want the low-level XDR implementation as part of the Gromacs
68 * documentation, so skip it for doxygen too...
70 #if (!defined HAVE_RPC_XDR_H && !defined DOXYGEN)
80 /* We need integer types that are guaranteed to be 4 bytes wide.
81 * If ANSI C99 headers were included they are already defined
82 * as int32_t and uint32_t. Check, and if not define them ourselves.
83 * Since it is just our workaround for missing ANSI C99 types, avoid adding
84 * it to the doxygen documentation.
86 #if !(defined INT32_MAX || defined DOXYGEN)
87 # if (INT_MAX == 2147483647)
89 # define uint32_t unsigned int
90 # define INT32_MAX 2147483647
91 # elif (LONG_MAX == 2147483647)
93 # define uint32_t unsigned long
94 # define INT32_MAX 2147483647L
96 # error ERROR: No 32 bit wide integer type found!
97 # error Use system XDR libraries instead, or update xdrfile.c
101 typedef struct XDR XDR
;
108 int (*x_getlong
) (XDR
*__xdrs
, int32_t *__lp
);
109 int (*x_putlong
) (XDR
*__xdrs
, int32_t *__lp
);
110 int (*x_getbytes
) (XDR
*__xdrs
, char *__addr
, unsigned int __len
);
111 int (*x_putbytes
) (XDR
*__xdrs
, char *__addr
, unsigned int __len
);
112 /* two next routines are not 64-bit IO safe - don't use! */
113 unsigned int (*x_getpostn
) (XDR
*__xdrs
);
114 int (*x_setpostn
) (XDR
*__xdrs
, unsigned int __pos
);
115 void (*x_destroy
) (XDR
*__xdrs
);
121 static int xdr_char (XDR
*xdrs
, char *ip
);
122 static int xdr_u_char (XDR
*xdrs
, unsigned char *ip
);
123 static int xdr_short (XDR
*xdrs
, short *ip
);
124 static int xdr_u_short (XDR
*xdrs
, unsigned short *ip
);
125 static int xdr_int (XDR
*xdrs
, int *ip
);
126 static int xdr_u_int (XDR
*xdrs
, unsigned int *ip
);
127 static int xdr_float (XDR
*xdrs
, float *ip
);
128 static int xdr_double (XDR
*xdrs
, double *ip
);
129 static int xdr_string (XDR
*xdrs
, char **ip
, unsigned int maxsize
);
130 static int xdr_opaque (XDR
*xdrs
, char *cp
, unsigned int cnt
);
131 static void xdrstdio_create (XDR
*xdrs
, FILE *fp
, enum xdr_op xop
);
133 #define xdr_getpos(xdrs) \
134 (*(xdrs)->x_ops->x_getpostn)(xdrs)
135 #define xdr_setpos(xdrs, pos) \
136 (*(xdrs)->x_ops->x_setpostn)(xdrs, pos)
137 #define xdr_destroy(xdrs) \
139 if ((xdrs)->x_ops->x_destroy) \
140 (*(xdrs)->x_ops->x_destroy)(xdrs); \
142 #endif /* end of our own XDR declarations */
148 /** Contents of the abstract XDRFILE data structure.
152 * This structure is used to provide an XDR file interface that is
153 * virtual identical to the standard UNIX fopen/fread/fwrite/fclose.
157 FILE * fp
; /**< pointer to standard C library file handle */
158 XDR
* xdr
; /**< pointer to corresponding XDR handle */
159 char mode
; /**< r=read, w=write, a=append */
160 int * buf1
; /**< Buffer for internal use */
161 int buf1size
; /**< Current allocated length of buf1 */
162 int * buf2
; /**< Buffer for internal use */
163 int buf2size
; /**< Current allocated length of buf2 */
169 /*************************************************************
170 * Implementation of higher-level routines to read/write *
171 * portable data based on the XDR standard. These should be *
172 * called from C - see further down for Fortran77 wrappers. *
173 *************************************************************/
176 xdrfile_open(const char *path
, const char *mode
)
182 /* make sure XDR files are opened in binary mode... */
183 if(*mode
=='w' || *mode
=='W')
185 sprintf(newmode
,"wb+");
187 } else if(*mode
== 'a' || *mode
== 'A')
189 sprintf(newmode
,"ab+");
190 xdrmode
= XDR_ENCODE
;
191 } else if(*mode
== 'r' || *mode
== 'R')
193 sprintf(newmode
,"rb");
194 xdrmode
= XDR_DECODE
;
195 } else /* cannot determine mode */
198 if((xfp
=(XDRFILE
*)malloc(sizeof(XDRFILE
)))==NULL
)
200 if((xfp
->fp
=fopen(path
,newmode
))==NULL
)
205 if((xfp
->xdr
=(XDR
*)malloc(sizeof(XDR
)))==NULL
)
212 xdrstdio_create((XDR
*)(xfp
->xdr
),xfp
->fp
,xdrmode
);
213 xfp
->buf1
= xfp
->buf2
= NULL
;
214 xfp
->buf1size
= xfp
->buf2size
= 0;
219 xdrfile_close(XDRFILE
*xfp
)
224 /* flush and destroy XDR stream */
226 xdr_destroy((XDR
*)(xfp
->xdr
));
236 return ret
; /* return 0 if ok */
242 xdrfile_read_int(int *ptr
, int ndata
, XDRFILE
* xfp
)
246 /* read write is encoded in the XDR struct */
247 while(i
<ndata
&& xdr_int((XDR
*)(xfp
->xdr
),ptr
+i
))
254 xdrfile_write_int(int *ptr
, int ndata
, XDRFILE
* xfp
)
258 /* read write is encoded in the XDR struct */
259 while(i
<ndata
&& xdr_int((XDR
*)(xfp
->xdr
),ptr
+i
))
266 xdrfile_read_uint(unsigned int *ptr
, int ndata
, XDRFILE
* xfp
)
270 /* read write is encoded in the XDR struct */
271 while(i
<ndata
&& xdr_u_int((XDR
*)(xfp
->xdr
),ptr
+i
))
278 xdrfile_write_uint(unsigned int *ptr
, int ndata
, XDRFILE
* xfp
)
282 /* read write is encoded in the XDR struct */
283 while(i
<ndata
&& xdr_u_int((XDR
*)(xfp
->xdr
),ptr
+i
))
289 xdrfile_read_char(char *ptr
, int ndata
, XDRFILE
* xfp
)
293 /* read write is encoded in the XDR struct */
294 while(i
<ndata
&& xdr_char((XDR
*)(xfp
->xdr
),ptr
+i
))
301 xdrfile_write_char(char *ptr
, int ndata
, XDRFILE
* xfp
)
305 /* read write is encoded in the XDR struct */
306 while(i
<ndata
&& xdr_char((XDR
*)(xfp
->xdr
),ptr
+i
))
313 xdrfile_read_uchar(unsigned char *ptr
, int ndata
, XDRFILE
* xfp
)
317 /* read write is encoded in the XDR struct */
318 while(i
<ndata
&& xdr_u_char((XDR
*)(xfp
->xdr
),ptr
+i
))
325 xdrfile_write_uchar(unsigned char *ptr
, int ndata
, XDRFILE
* xfp
)
329 /* read write is encoded in the XDR struct */
330 while(i
<ndata
&& xdr_u_char((XDR
*)(xfp
->xdr
),ptr
+i
))
336 xdrfile_read_short(short *ptr
, int ndata
, XDRFILE
* xfp
)
340 /* read write is encoded in the XDR struct */
341 while(i
<ndata
&& xdr_short((XDR
*)(xfp
->xdr
),ptr
+i
))
348 xdrfile_write_short(short *ptr
, int ndata
, XDRFILE
* xfp
)
352 /* read write is encoded in the XDR struct */
353 while(i
<ndata
&& xdr_short((XDR
*)(xfp
->xdr
),ptr
+i
))
360 xdrfile_read_ushort(unsigned short *ptr
, int ndata
, XDRFILE
* xfp
)
364 /* read write is encoded in the XDR struct */
365 while(i
<ndata
&& xdr_u_short((XDR
*)(xfp
->xdr
),ptr
+i
))
372 xdrfile_write_ushort(unsigned short *ptr
, int ndata
, XDRFILE
* xfp
)
376 /* read write is encoded in the XDR struct */
377 while(i
<ndata
&& xdr_u_short((XDR
*)(xfp
->xdr
),ptr
+i
))
383 xdrfile_read_float(float *ptr
, int ndata
, XDRFILE
* xfp
)
386 /* read write is encoded in the XDR struct */
387 while(i
<ndata
&& xdr_float((XDR
*)(xfp
->xdr
),ptr
+i
))
393 xdrfile_write_float(float *ptr
, int ndata
, XDRFILE
* xfp
)
396 /* read write is encoded in the XDR struct */
397 while(i
<ndata
&& xdr_float((XDR
*)(xfp
->xdr
),ptr
+i
))
403 xdrfile_read_double(double *ptr
, int ndata
, XDRFILE
* xfp
)
406 /* read write is encoded in the XDR struct */
407 while(i
<ndata
&& xdr_double((XDR
*)(xfp
->xdr
),ptr
+i
))
413 xdrfile_write_double(double *ptr
, int ndata
, XDRFILE
* xfp
)
416 /* read write is encoded in the XDR struct */
417 while(i
<ndata
&& xdr_double((XDR
*)(xfp
->xdr
),ptr
+i
))
423 xdrfile_read_string(char *ptr
, int maxlen
, XDRFILE
* xfp
)
426 if(xdr_string((XDR
*)(xfp
->xdr
),&ptr
,maxlen
)) {
428 while(i
<maxlen
&& ptr
[i
]!=0)
439 xdrfile_write_string(char *ptr
, XDRFILE
* xfp
)
441 int len
=strlen(ptr
)+1;
443 if(xdr_string((XDR
*)(xfp
->xdr
),&ptr
,len
))
451 xdrfile_read_opaque(char *ptr
, int cnt
, XDRFILE
* xfp
)
453 if(xdr_opaque((XDR
*)(xfp
->xdr
),ptr
,cnt
))
461 xdrfile_write_opaque(char *ptr
, int cnt
, XDRFILE
* xfp
)
463 if(xdr_opaque((XDR
*)(xfp
->xdr
),ptr
,cnt
))
470 /* Internal support routines for reading/writing compressed coordinates
471 * sizeofint - calculate smallest number of bits necessary
472 * to represent a certain integer.
475 sizeofint(int size
) {
476 unsigned int num
= 1;
479 while (size
>= num
&& num_of_bits
< 32)
489 * sizeofints - calculate 'bitsize' of compressed ints
491 * given a number of small unsigned integers and the maximum value
492 * return the number of bits needed to read or write them with the
493 * routines encodeints/decodeints. You need this parameter when
494 * calling those routines.
495 * (However, in some cases we can just use the variable 'smallidx'
496 * which is the exact number of bits, and them we dont need to call
500 sizeofints(int num_of_ints
, unsigned int sizes
[])
503 unsigned int num_of_bytes
, num_of_bits
, bytes
[32], bytecnt
, tmp
;
507 for (i
=0; i
< num_of_ints
; i
++)
510 for (bytecnt
= 0; bytecnt
< num_of_bytes
; bytecnt
++)
512 tmp
= bytes
[bytecnt
] * sizes
[i
] + tmp
;
513 bytes
[bytecnt
] = tmp
& 0xff;
518 bytes
[bytecnt
++] = tmp
& 0xff;
521 num_of_bytes
= bytecnt
;
525 while (bytes
[num_of_bytes
] >= num
)
530 return num_of_bits
+ num_of_bytes
* 8;
536 * encodebits - encode num into buf using the specified number of bits
538 * This routines appends the value of num to the bits already present in
539 * the array buf. You need to give it the number of bits to use and you had
540 * better make sure that this number of bits is enough to hold the value.
541 * Num must also be positive.
544 encodebits(int buf
[], int num_of_bits
, int num
)
547 unsigned int cnt
, lastbyte
;
549 unsigned char * cbuf
;
551 cbuf
= ((unsigned char *)buf
) + 3 * sizeof(*buf
);
552 cnt
= (unsigned int) buf
[0];
554 lastbyte
=(unsigned int) buf
[2];
555 while (num_of_bits
>= 8)
557 lastbyte
= (lastbyte
<< 8) | ((num
>> (num_of_bits
-8)) /* & 0xff*/);
558 cbuf
[cnt
++] = lastbyte
>> lastbits
;
563 lastbyte
= (lastbyte
<< num_of_bits
) | num
;
564 lastbits
+= num_of_bits
;
568 cbuf
[cnt
++] = lastbyte
>> lastbits
;
576 cbuf
[cnt
] = lastbyte
<< (8 - lastbits
);
581 * encodeints - encode a small set of small integers in compressed format
583 * this routine is used internally by xdr3dfcoord, to encode a set of
584 * small integers to the buffer for writing to a file.
585 * Multiplication with fixed (specified maximum) sizes is used to get
586 * to one big, multibyte integer. Allthough the routine could be
587 * modified to handle sizes bigger than 16777216, or more than just
588 * a few integers, this is not done because the gain in compression
589 * isn't worth the effort. Note that overflowing the multiplication
590 * or the byte buffer (32 bytes) is unchecked and whould cause bad results.
591 * THese things are checked in the calling routines, so make sure not
592 * to remove those checks...
596 encodeints(int buf
[], int num_of_ints
, int num_of_bits
,
597 unsigned int sizes
[], unsigned int nums
[])
601 unsigned int bytes
[32], num_of_bytes
, bytecnt
, tmp
;
607 bytes
[num_of_bytes
++] = tmp
& 0xff;
611 for (i
= 1; i
< num_of_ints
; i
++)
613 if (nums
[i
] >= sizes
[i
])
615 fprintf(stderr
,"major breakdown in encodeints - num %u doesn't "
616 "match size %u\n", nums
[i
], sizes
[i
]);
619 /* use one step multiply */
621 for (bytecnt
= 0; bytecnt
< num_of_bytes
; bytecnt
++)
623 tmp
= bytes
[bytecnt
] * sizes
[i
] + tmp
;
624 bytes
[bytecnt
] = tmp
& 0xff;
629 bytes
[bytecnt
++] = tmp
& 0xff;
632 num_of_bytes
= bytecnt
;
634 if (num_of_bits
>= num_of_bytes
* 8)
636 for (i
= 0; i
< num_of_bytes
; i
++)
638 encodebits(buf
, 8, bytes
[i
]);
640 encodebits(buf
, num_of_bits
- num_of_bytes
* 8, 0);
644 for (i
= 0; i
< num_of_bytes
-1; i
++)
646 encodebits(buf
, 8, bytes
[i
]);
648 encodebits(buf
, num_of_bits
- (num_of_bytes
-1) * 8, bytes
[i
]);
654 * decodebits - decode number from buf using specified number of bits
656 * extract the number of bits from the array buf and construct an integer
657 * from it. Return that value.
662 decodebits(int buf
[], int num_of_bits
)
666 unsigned int lastbits
, lastbyte
;
667 unsigned char * cbuf
;
668 int mask
= (1 << num_of_bits
) -1;
670 cbuf
= ((unsigned char *)buf
) + 3 * sizeof(*buf
);
672 lastbits
= (unsigned int) buf
[1];
673 lastbyte
= (unsigned int) buf
[2];
676 while (num_of_bits
>= 8)
678 lastbyte
= ( lastbyte
<< 8 ) | cbuf
[cnt
++];
679 num
|= (lastbyte
>> lastbits
) << (num_of_bits
- 8);
684 if (lastbits
< num_of_bits
)
687 lastbyte
= (lastbyte
<< 8) | cbuf
[cnt
++];
689 lastbits
-= num_of_bits
;
690 num
|= (lastbyte
>> lastbits
) & ((1 << num_of_bits
) -1);
700 * decodeints - decode 'small' integers from the buf array
702 * this routine is the inverse from encodeints() and decodes the small integers
703 * written to buf by calculating the remainder and doing divisions with
704 * the given sizes[]. You need to specify the total number of bits to be
705 * used from buf in num_of_bits.
710 decodeints(int buf
[], int num_of_ints
, int num_of_bits
,
711 unsigned int sizes
[], int nums
[])
715 int i
, j
, num_of_bytes
, p
, num
;
717 bytes
[1] = bytes
[2] = bytes
[3] = 0;
719 while (num_of_bits
> 8)
721 bytes
[num_of_bytes
++] = decodebits(buf
, 8);
726 bytes
[num_of_bytes
++] = decodebits(buf
, num_of_bits
);
728 for (i
= num_of_ints
-1; i
> 0; i
--)
731 for (j
= num_of_bytes
-1; j
>=0; j
--)
733 num
= (num
<< 8) | bytes
[j
];
736 num
= num
- p
* sizes
[i
];
740 nums
[0] = bytes
[0] | (bytes
[1] << 8) | (bytes
[2] << 16) | (bytes
[3] << 24);
744 static const int magicints
[] =
746 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
747 80, 101, 128, 161, 203, 256, 322, 406, 512, 645, 812, 1024, 1290,
748 1625, 2048, 2580, 3250, 4096, 5060, 6501, 8192, 10321, 13003,
749 16384, 20642, 26007, 32768, 41285, 52015, 65536,82570, 104031,
750 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
751 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021,
752 4194304, 5284491, 6658042, 8388607, 10568983, 13316085, 16777216
756 /* note that magicints[FIRSTIDX-1] == 0 */
757 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
759 /* Compressed coordinate routines - modified from the original
760 * implementation by Frans v. Hoesel to make them threadsafe.
763 xdrfile_decompress_coord_float(float *ptr
,
768 int minint
[3], maxint
[3], *lip
;
769 int smallidx
, minidx
, maxidx
;
770 unsigned sizeint
[3], sizesmall
[3], bitsizeint
[3], size3
;
771 int k
, *buf1
, *buf2
, lsize
, flag
;
772 int smallnum
, smaller
, larger
, i
, is_smaller
, run
;
773 float *lfp
, inv_precision
;
774 int tmp
, *thiscoord
, prevcoord
[3];
775 unsigned int bitsize
;
781 if(xfp
==NULL
|| ptr
==NULL
)
783 tmp
=xdrfile_read_int(&lsize
,1,xfp
);
785 return -1; /* return if we could not read size */
788 fprintf(stderr
, "Requested to decompress %d coords, file contains %d\n",
794 if(size3
>xfp
->buf1size
)
796 if((xfp
->buf1
=(int *)malloc(sizeof(int)*size3
))==NULL
)
798 fprintf(stderr
,"Cannot allocate memory for decompressing coordinates.\n");
802 xfp
->buf2size
=size3
*1.2;
803 if((xfp
->buf2
=(int *)malloc(sizeof(int)*xfp
->buf2size
))==NULL
)
805 fprintf(stderr
,"Cannot allocate memory for decompressing coordinates.\n");
809 /* Dont bother with compression for three atoms or less */
812 return xdrfile_read_float(ptr
,size3
,xfp
)/3;
813 /* return number of coords, not floats */
815 /* Compression-time if we got here. Read precision first */
816 xdrfile_read_float(precision
,1,xfp
);
818 /* avoid repeated pointer dereferencing. */
821 /* buf2[0-2] are special and do not contain actual data */
822 buf2
[0] = buf2
[1] = buf2
[2] = 0;
823 xdrfile_read_int(minint
,3,xfp
);
824 xdrfile_read_int(maxint
,3,xfp
);
826 sizeint
[0] = maxint
[0] - minint
[0]+1;
827 sizeint
[1] = maxint
[1] - minint
[1]+1;
828 sizeint
[2] = maxint
[2] - minint
[2]+1;
830 /* check if one of the sizes is to big to be multiplied */
831 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff)
833 bitsizeint
[0] = sizeofint(sizeint
[0]);
834 bitsizeint
[1] = sizeofint(sizeint
[1]);
835 bitsizeint
[2] = sizeofint(sizeint
[2]);
836 bitsize
= 0; /* flag the use of large sizes */
840 bitsize
= sizeofints(3, sizeint
);
843 if (xdrfile_read_int(&smallidx
,1,xfp
) == 0)
844 return 0; /* not sure what has happened here or why we return... */
846 maxidx
= (LASTIDX
<tmp
) ? LASTIDX
: tmp
;
847 minidx
= maxidx
- 8; /* often this equal smallidx */
849 tmp
= (FIRSTIDX
>tmp
) ? FIRSTIDX
: tmp
;
850 smaller
= magicints
[tmp
] / 2;
851 smallnum
= magicints
[smallidx
] / 2;
852 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
] ;
853 larger
= magicints
[maxidx
];
855 /* buf2[0] holds the length in bytes */
857 if (xdrfile_read_int(buf2
,1,xfp
) == 0)
859 if (xdrfile_read_opaque((char *)&(buf2
[3]),(unsigned int)buf2
[0],xfp
) == 0)
861 buf2
[0] = buf2
[1] = buf2
[2] = 0;
864 inv_precision
= 1.0 / * precision
;
870 thiscoord
= (int *)(lip
) + i
* 3;
874 thiscoord
[0] = decodebits(buf2
, bitsizeint
[0]);
875 thiscoord
[1] = decodebits(buf2
, bitsizeint
[1]);
876 thiscoord
[2] = decodebits(buf2
, bitsizeint
[2]);
880 decodeints(buf2
, 3, bitsize
, sizeint
, thiscoord
);
884 thiscoord
[0] += minint
[0];
885 thiscoord
[1] += minint
[1];
886 thiscoord
[2] += minint
[2];
888 prevcoord
[0] = thiscoord
[0];
889 prevcoord
[1] = thiscoord
[1];
890 prevcoord
[2] = thiscoord
[2];
892 flag
= decodebits(buf2
, 1);
896 run
= decodebits(buf2
, 5);
897 is_smaller
= run
% 3;
904 for (k
= 0; k
< run
; k
+=3)
906 decodeints(buf2
, 3, smallidx
, sizesmall
, thiscoord
);
908 thiscoord
[0] += prevcoord
[0] - smallnum
;
909 thiscoord
[1] += prevcoord
[1] - smallnum
;
910 thiscoord
[2] += prevcoord
[2] - smallnum
;
912 /* interchange first with second atom for better
913 * compression of water molecules
915 tmp
= thiscoord
[0]; thiscoord
[0] = prevcoord
[0];
917 tmp
= thiscoord
[1]; thiscoord
[1] = prevcoord
[1];
919 tmp
= thiscoord
[2]; thiscoord
[2] = prevcoord
[2];
921 *lfp
++ = prevcoord
[0] * inv_precision
;
922 *lfp
++ = prevcoord
[1] * inv_precision
;
923 *lfp
++ = prevcoord
[2] * inv_precision
;
925 prevcoord
[0] = thiscoord
[0];
926 prevcoord
[1] = thiscoord
[1];
927 prevcoord
[2] = thiscoord
[2];
929 *lfp
++ = thiscoord
[0] * inv_precision
;
930 *lfp
++ = thiscoord
[1] * inv_precision
;
931 *lfp
++ = thiscoord
[2] * inv_precision
;
936 *lfp
++ = thiscoord
[0] * inv_precision
;
937 *lfp
++ = thiscoord
[1] * inv_precision
;
938 *lfp
++ = thiscoord
[2] * inv_precision
;
940 smallidx
+= is_smaller
;
945 if (smallidx
> FIRSTIDX
)
947 smaller
= magicints
[smallidx
- 1] /2;
954 else if (is_smaller
> 0)
957 smallnum
= magicints
[smallidx
] / 2;
959 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
] ;
965 xdrfile_compress_coord_float(float *ptr
,
970 int minint
[3], maxint
[3], mindiff
, *lip
, diff
;
971 int lint1
, lint2
, lint3
, oldlint1
, oldlint2
, oldlint3
, smallidx
;
973 unsigned sizeint
[3], sizesmall
[3], bitsizeint
[3], size3
, *luip
;
975 int smallnum
, smaller
, larger
, i
, j
, is_small
, is_smaller
, run
, prevrun
;
977 int tmp
, tmpsum
, *thiscoord
, prevcoord
[3];
978 unsigned int tmpcoord
[30];
980 unsigned int bitsize
;
990 if(size3
>xfp
->buf1size
)
992 if((xfp
->buf1
=(int *)malloc(sizeof(int)*size3
))==NULL
)
994 fprintf(stderr
,"Cannot allocate memory for compressing coordinates.\n");
998 xfp
->buf2size
=size3
*1.2;
999 if((xfp
->buf2
=(int *)malloc(sizeof(int)*xfp
->buf2size
))==NULL
)
1001 fprintf(stderr
,"Cannot allocate memory for compressing coordinates.\n");
1005 if(xdrfile_write_int(&size
,1,xfp
)==0)
1006 return -1; /* return if we could not write size */
1007 /* Dont bother with compression for three atoms or less */
1010 return xdrfile_write_float(ptr
,size3
,xfp
)/3;
1011 /* return number of coords, not floats */
1013 /* Compression-time if we got here. Write precision first */
1016 xdrfile_write_float(&precision
,1,xfp
);
1017 /* avoid repeated pointer dereferencing. */
1020 /* buf2[0-2] are special and do not contain actual data */
1021 buf2
[0] = buf2
[1] = buf2
[2] = 0;
1022 minint
[0] = minint
[1] = minint
[2] = INT_MAX
;
1023 maxint
[0] = maxint
[1] = maxint
[2] = INT_MIN
;
1028 oldlint1
= oldlint2
= oldlint3
= 0;
1029 while(lfp
< ptr
+ size3
)
1031 /* find nearest integer */
1033 lf
= *lfp
* precision
+ 0.5;
1035 lf
= *lfp
* precision
- 0.5;
1036 if (fabs(lf
) > INT_MAX
-2)
1038 /* scaling would cause overflow */
1039 fprintf(stderr
,"Internal overflow compressing coordinates.\n");
1043 if (lint1
< minint
[0]) minint
[0] = lint1
;
1044 if (lint1
> maxint
[0]) maxint
[0] = lint1
;
1048 lf
= *lfp
* precision
+ 0.5;
1050 lf
= *lfp
* precision
- 0.5;
1051 if (fabs(lf
) > INT_MAX
-2)
1053 /* scaling would cause overflow */
1054 fprintf(stderr
,"Internal overflow compressing coordinates.\n");
1058 if (lint2
< minint
[1]) minint
[1] = lint2
;
1059 if (lint2
> maxint
[1]) maxint
[1] = lint2
;
1063 lf
= *lfp
* precision
+ 0.5;
1065 lf
= *lfp
* precision
- 0.5;
1066 if (fabs(lf
) > INT_MAX
-2)
1071 if (lint3
< minint
[2]) minint
[2] = lint3
;
1072 if (lint3
> maxint
[2]) maxint
[2] = lint3
;
1075 diff
= abs(oldlint1
-lint1
)+abs(oldlint2
-lint2
)+abs(oldlint3
-lint3
);
1076 if (diff
< mindiff
&& lfp
> ptr
+ 3)
1082 xdrfile_write_int(minint
,3,xfp
);
1083 xdrfile_write_int(maxint
,3,xfp
);
1085 if ((float)maxint
[0] - (float)minint
[0] >= INT_MAX
-2 ||
1086 (float)maxint
[1] - (float)minint
[1] >= INT_MAX
-2 ||
1087 (float)maxint
[2] - (float)minint
[2] >= INT_MAX
-2) {
1088 /* turning value in unsigned by subtracting minint
1089 * would cause overflow
1091 fprintf(stderr
,"Internal overflow compressing coordinates.\n");
1094 sizeint
[0] = maxint
[0] - minint
[0]+1;
1095 sizeint
[1] = maxint
[1] - minint
[1]+1;
1096 sizeint
[2] = maxint
[2] - minint
[2]+1;
1098 /* check if one of the sizes is to big to be multiplied */
1099 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff)
1101 bitsizeint
[0] = sizeofint(sizeint
[0]);
1102 bitsizeint
[1] = sizeofint(sizeint
[1]);
1103 bitsizeint
[2] = sizeofint(sizeint
[2]);
1104 bitsize
= 0; /* flag the use of large sizes */
1108 bitsize
= sizeofints(3, sizeint
);
1111 luip
= (unsigned int *) buf1
;
1112 smallidx
= FIRSTIDX
;
1113 while (smallidx
< LASTIDX
&& magicints
[smallidx
] < mindiff
)
1117 xdrfile_write_int(&smallidx
,1,xfp
);
1119 maxidx
= (LASTIDX
<tmp
) ? LASTIDX
: tmp
;
1120 minidx
= maxidx
- 8; /* often this equal smallidx */
1122 tmp
= (FIRSTIDX
>tmp
) ? FIRSTIDX
: tmp
;
1123 smaller
= magicints
[tmp
] / 2;
1124 smallnum
= magicints
[smallidx
] / 2;
1125 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
1126 larger
= magicints
[maxidx
] / 2;
1127 prevcoord
[0] = prevcoord
[1] = prevcoord
[2] = 0;
1132 thiscoord
= (int *)(luip
) + i
* 3;
1133 if (smallidx
< maxidx
&& i
>= 1 &&
1134 abs(thiscoord
[0] - prevcoord
[0]) < larger
&&
1135 abs(thiscoord
[1] - prevcoord
[1]) < larger
&&
1136 abs(thiscoord
[2] - prevcoord
[2]) < larger
) {
1139 else if (smallidx
> minidx
)
1149 if (abs(thiscoord
[0] - thiscoord
[3]) < smallnum
&&
1150 abs(thiscoord
[1] - thiscoord
[4]) < smallnum
&&
1151 abs(thiscoord
[2] - thiscoord
[5]) < smallnum
)
1153 /* interchange first with second atom for better
1154 * compression of water molecules
1156 tmp
= thiscoord
[0]; thiscoord
[0] = thiscoord
[3];
1158 tmp
= thiscoord
[1]; thiscoord
[1] = thiscoord
[4];
1160 tmp
= thiscoord
[2]; thiscoord
[2] = thiscoord
[5];
1165 tmpcoord
[0] = thiscoord
[0] - minint
[0];
1166 tmpcoord
[1] = thiscoord
[1] - minint
[1];
1167 tmpcoord
[2] = thiscoord
[2] - minint
[2];
1170 encodebits(buf2
, bitsizeint
[0], tmpcoord
[0]);
1171 encodebits(buf2
, bitsizeint
[1], tmpcoord
[1]);
1172 encodebits(buf2
, bitsizeint
[2], tmpcoord
[2]);
1176 encodeints(buf2
, 3, bitsize
, sizeint
, tmpcoord
);
1178 prevcoord
[0] = thiscoord
[0];
1179 prevcoord
[1] = thiscoord
[1];
1180 prevcoord
[2] = thiscoord
[2];
1181 thiscoord
= thiscoord
+ 3;
1185 if (is_small
== 0 && is_smaller
== -1)
1187 while (is_small
&& run
< 8*3)
1192 tmp
=thiscoord
[j
] - prevcoord
[j
];
1195 if (is_smaller
== -1 && tmpsum
>= smaller
* smaller
)
1200 tmpcoord
[run
++] = thiscoord
[0] - prevcoord
[0] + smallnum
;
1201 tmpcoord
[run
++] = thiscoord
[1] - prevcoord
[1] + smallnum
;
1202 tmpcoord
[run
++] = thiscoord
[2] - prevcoord
[2] + smallnum
;
1204 prevcoord
[0] = thiscoord
[0];
1205 prevcoord
[1] = thiscoord
[1];
1206 prevcoord
[2] = thiscoord
[2];
1209 thiscoord
= thiscoord
+ 3;
1212 abs(thiscoord
[0] - prevcoord
[0]) < smallnum
&&
1213 abs(thiscoord
[1] - prevcoord
[1]) < smallnum
&&
1214 abs(thiscoord
[2] - prevcoord
[2]) < smallnum
)
1219 if (run
!= prevrun
|| is_smaller
!= 0)
1222 encodebits(buf2
, 1, 1); /* flag the change in run-length */
1223 encodebits(buf2
, 5, run
+is_smaller
+1);
1227 encodebits(buf2
, 1, 0); /* flag the fact that runlength did not change */
1229 for (k
=0; k
< run
; k
+=3)
1231 encodeints(buf2
, 3, smallidx
, sizesmall
, &tmpcoord
[k
]);
1233 if (is_smaller
!= 0)
1235 smallidx
+= is_smaller
;
1239 smaller
= magicints
[smallidx
-1] / 2;
1244 smallnum
= magicints
[smallidx
] / 2;
1246 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
1249 if (buf2
[1] != 0) buf2
[0]++;
1250 xdrfile_write_int(buf2
,1,xfp
); /* buf2[0] holds the length in bytes */
1251 tmp
=xdrfile_write_opaque((char *)&(buf2
[3]),(unsigned int)buf2
[0],xfp
);
1252 if(tmp
==(unsigned int)buf2
[0])
1260 xdrfile_decompress_coord_double(double *ptr
,
1265 int minint
[3], maxint
[3], *lip
;
1266 int smallidx
, minidx
, maxidx
;
1267 unsigned sizeint
[3], sizesmall
[3], bitsizeint
[3], size3
;
1268 int k
, *buf1
, *buf2
, lsize
, flag
;
1269 int smallnum
, smaller
, larger
, i
, is_smaller
, run
;
1270 double *lfp
, inv_precision
;
1271 float float_prec
, tmpdata
[30];
1272 int tmp
, *thiscoord
, prevcoord
[3];
1273 unsigned int bitsize
;
1279 if(xfp
==NULL
|| ptr
==NULL
)
1281 tmp
=xdrfile_read_int(&lsize
,1,xfp
);
1283 return -1; /* return if we could not read size */
1286 fprintf(stderr
, "Requested to decompress %d coords, file contains %d\n",
1292 if(size3
>xfp
->buf1size
)
1294 if((xfp
->buf1
=(int *)malloc(sizeof(int)*size3
))==NULL
)
1296 fprintf(stderr
,"Cannot allocate memory for decompression coordinates.\n");
1299 xfp
->buf1size
=size3
;
1300 xfp
->buf2size
=size3
*1.2;
1301 if((xfp
->buf2
=(int *)malloc(sizeof(int)*xfp
->buf2size
))==NULL
)
1303 fprintf(stderr
,"Cannot allocate memory for decompressing coordinates.\n");
1307 /* Dont bother with compression for three atoms or less */
1310 tmp
=xdrfile_read_float(tmpdata
,size3
,xfp
);
1314 /* return number of coords, not floats */
1316 /* Compression-time if we got here. Read precision first */
1317 xdrfile_read_float(&float_prec
,1,xfp
);
1318 *precision
=float_prec
;
1319 /* avoid repeated pointer dereferencing. */
1322 /* buf2[0-2] are special and do not contain actual data */
1323 buf2
[0] = buf2
[1] = buf2
[2] = 0;
1324 xdrfile_read_int(minint
,3,xfp
);
1325 xdrfile_read_int(maxint
,3,xfp
);
1327 sizeint
[0] = maxint
[0] - minint
[0]+1;
1328 sizeint
[1] = maxint
[1] - minint
[1]+1;
1329 sizeint
[2] = maxint
[2] - minint
[2]+1;
1331 /* check if one of the sizes is to big to be multiplied */
1332 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff)
1334 bitsizeint
[0] = sizeofint(sizeint
[0]);
1335 bitsizeint
[1] = sizeofint(sizeint
[1]);
1336 bitsizeint
[2] = sizeofint(sizeint
[2]);
1337 bitsize
= 0; /* flag the use of large sizes */
1341 bitsize
= sizeofints(3, sizeint
);
1344 if (xdrfile_read_int(&smallidx
,1,xfp
) == 0)
1347 maxidx
= (LASTIDX
<tmp
) ? LASTIDX
: tmp
;
1348 minidx
= maxidx
- 8; /* often this equal smallidx */
1350 tmp
= (FIRSTIDX
>tmp
) ? FIRSTIDX
: tmp
;
1351 smaller
= magicints
[tmp
] / 2;
1352 smallnum
= magicints
[smallidx
] / 2;
1353 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
] ;
1354 larger
= magicints
[maxidx
];
1356 /* buf2[0] holds the length in bytes */
1358 if (xdrfile_read_int(buf2
,1,xfp
) == 0)
1360 if (xdrfile_read_opaque((char *)&(buf2
[3]),(unsigned int)buf2
[0],xfp
) == 0)
1362 buf2
[0] = buf2
[1] = buf2
[2] = 0;
1365 inv_precision
= 1.0 / * precision
;
1371 thiscoord
= (int *)(lip
) + i
* 3;
1375 thiscoord
[0] = decodebits(buf2
, bitsizeint
[0]);
1376 thiscoord
[1] = decodebits(buf2
, bitsizeint
[1]);
1377 thiscoord
[2] = decodebits(buf2
, bitsizeint
[2]);
1379 decodeints(buf2
, 3, bitsize
, sizeint
, thiscoord
);
1383 thiscoord
[0] += minint
[0];
1384 thiscoord
[1] += minint
[1];
1385 thiscoord
[2] += minint
[2];
1387 prevcoord
[0] = thiscoord
[0];
1388 prevcoord
[1] = thiscoord
[1];
1389 prevcoord
[2] = thiscoord
[2];
1391 flag
= decodebits(buf2
, 1);
1395 run
= decodebits(buf2
, 5);
1396 is_smaller
= run
% 3;
1403 for (k
= 0; k
< run
; k
+=3)
1405 decodeints(buf2
, 3, smallidx
, sizesmall
, thiscoord
);
1407 thiscoord
[0] += prevcoord
[0] - smallnum
;
1408 thiscoord
[1] += prevcoord
[1] - smallnum
;
1409 thiscoord
[2] += prevcoord
[2] - smallnum
;
1412 /* interchange first with second atom for better
1413 * compression of water molecules
1415 tmp
= thiscoord
[0]; thiscoord
[0] = prevcoord
[0];
1417 tmp
= thiscoord
[1]; thiscoord
[1] = prevcoord
[1];
1419 tmp
= thiscoord
[2]; thiscoord
[2] = prevcoord
[2];
1421 *lfp
++ = prevcoord
[0] * inv_precision
;
1422 *lfp
++ = prevcoord
[1] * inv_precision
;
1423 *lfp
++ = prevcoord
[2] * inv_precision
;
1427 prevcoord
[0] = thiscoord
[0];
1428 prevcoord
[1] = thiscoord
[1];
1429 prevcoord
[2] = thiscoord
[2];
1431 *lfp
++ = thiscoord
[0] * inv_precision
;
1432 *lfp
++ = thiscoord
[1] * inv_precision
;
1433 *lfp
++ = thiscoord
[2] * inv_precision
;
1436 *lfp
++ = thiscoord
[0] * inv_precision
;
1437 *lfp
++ = thiscoord
[1] * inv_precision
;
1438 *lfp
++ = thiscoord
[2] * inv_precision
;
1440 smallidx
+= is_smaller
;
1441 if (is_smaller
< 0) {
1443 if (smallidx
> FIRSTIDX
) {
1444 smaller
= magicints
[smallidx
- 1] /2;
1448 } else if (is_smaller
> 0) {
1450 smallnum
= magicints
[smallidx
] / 2;
1452 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
] ;
1458 xdrfile_compress_coord_double(double *ptr
,
1463 int minint
[3], maxint
[3], mindiff
, *lip
, diff
;
1464 int lint1
, lint2
, lint3
, oldlint1
, oldlint2
, oldlint3
, smallidx
;
1466 unsigned sizeint
[3], sizesmall
[3], bitsizeint
[3], size3
, *luip
;
1467 int k
, *buf1
, *buf2
;
1468 int smallnum
, smaller
, larger
, i
, j
, is_small
, is_smaller
, run
, prevrun
;
1470 float float_prec
, lf
,tmpdata
[30];
1471 int tmp
, tmpsum
, *thiscoord
, prevcoord
[3];
1472 unsigned int tmpcoord
[30];
1474 unsigned int bitsize
;
1483 if(size3
>xfp
->buf1size
) {
1484 if((xfp
->buf1
=(int *)malloc(sizeof(int)*size3
))==NULL
) {
1485 fprintf(stderr
,"Cannot allocate memory for compressing coordinates.\n");
1488 xfp
->buf1size
=size3
;
1489 xfp
->buf2size
=size3
*1.2;
1490 if((xfp
->buf2
=(int *)malloc(sizeof(int)*xfp
->buf2size
))==NULL
) {
1491 fprintf(stderr
,"Cannot allocate memory for compressing coordinates.\n");
1495 if(xdrfile_write_int(&size
,1,xfp
)==0)
1496 return -1; /* return if we could not write size */
1497 /* Dont bother with compression for three atoms or less */
1501 return xdrfile_write_float(tmpdata
,size3
,xfp
)/3;
1502 /* return number of coords, not floats */
1504 /* Compression-time if we got here. Write precision first */
1507 float_prec
=precision
;
1508 xdrfile_write_float(&float_prec
,1,xfp
);
1509 /* avoid repeated pointer dereferencing. */
1512 /* buf2[0-2] are special and do not contain actual data */
1513 buf2
[0] = buf2
[1] = buf2
[2] = 0;
1514 minint
[0] = minint
[1] = minint
[2] = INT_MAX
;
1515 maxint
[0] = maxint
[1] = maxint
[2] = INT_MIN
;
1520 oldlint1
= oldlint2
= oldlint3
= 0;
1521 while(lfp
< ptr
+ size3
) {
1522 /* find nearest integer */
1524 lf
= (float)*lfp
* float_prec
+ 0.5;
1526 lf
= (float)*lfp
* float_prec
- 0.5;
1527 if (fabs(lf
) > INT_MAX
-2) {
1528 /* scaling would cause overflow */
1529 fprintf(stderr
,"Internal overflow compressing coordinates.\n");
1533 if (lint1
< minint
[0]) minint
[0] = lint1
;
1534 if (lint1
> maxint
[0]) maxint
[0] = lint1
;
1538 lf
= (float)*lfp
* float_prec
+ 0.5;
1540 lf
= (float)*lfp
* float_prec
- 0.5;
1541 if (fabs(lf
) > INT_MAX
-2) {
1542 /* scaling would cause overflow */
1543 fprintf(stderr
,"Internal overflow compressing coordinates.\n");
1547 if (lint2
< minint
[1]) minint
[1] = lint2
;
1548 if (lint2
> maxint
[1]) maxint
[1] = lint2
;
1552 lf
= (float)*lfp
* float_prec
+ 0.5;
1554 lf
= (float)*lfp
* float_prec
- 0.5;
1555 if (fabs(lf
) > INT_MAX
-2) {
1559 if (lint3
< minint
[2]) minint
[2] = lint3
;
1560 if (lint3
> maxint
[2]) maxint
[2] = lint3
;
1563 diff
= abs(oldlint1
-lint1
)+abs(oldlint2
-lint2
)+abs(oldlint3
-lint3
);
1564 if (diff
< mindiff
&& lfp
> ptr
+ 3)
1570 xdrfile_write_int(minint
,3,xfp
);
1571 xdrfile_write_int(maxint
,3,xfp
);
1573 if ((float)maxint
[0] - (float)minint
[0] >= INT_MAX
-2 ||
1574 (float)maxint
[1] - (float)minint
[1] >= INT_MAX
-2 ||
1575 (float)maxint
[2] - (float)minint
[2] >= INT_MAX
-2) {
1576 /* turning value in unsigned by subtracting minint
1577 * would cause overflow
1579 fprintf(stderr
,"Internal overflow compressing coordinates.\n");
1582 sizeint
[0] = maxint
[0] - minint
[0]+1;
1583 sizeint
[1] = maxint
[1] - minint
[1]+1;
1584 sizeint
[2] = maxint
[2] - minint
[2]+1;
1586 /* check if one of the sizes is to big to be multiplied */
1587 if ((sizeint
[0] | sizeint
[1] | sizeint
[2] ) > 0xffffff) {
1588 bitsizeint
[0] = sizeofint(sizeint
[0]);
1589 bitsizeint
[1] = sizeofint(sizeint
[1]);
1590 bitsizeint
[2] = sizeofint(sizeint
[2]);
1591 bitsize
= 0; /* flag the use of large sizes */
1593 bitsize
= sizeofints(3, sizeint
);
1596 luip
= (unsigned int *) buf1
;
1597 smallidx
= FIRSTIDX
;
1598 while (smallidx
< LASTIDX
&& magicints
[smallidx
] < mindiff
) {
1601 xdrfile_write_int(&smallidx
,1,xfp
);
1603 maxidx
= (LASTIDX
<tmp
) ? LASTIDX
: tmp
;
1604 minidx
= maxidx
- 8; /* often this equal smallidx */
1606 tmp
= (FIRSTIDX
>tmp
) ? FIRSTIDX
: tmp
;
1607 smaller
= magicints
[tmp
] / 2;
1608 smallnum
= magicints
[smallidx
] / 2;
1609 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
1610 larger
= magicints
[maxidx
] / 2;
1611 prevcoord
[0] = prevcoord
[1] = prevcoord
[2] = 0;
1615 thiscoord
= (int *)(luip
) + i
* 3;
1616 if (smallidx
< maxidx
&& i
>= 1 &&
1617 abs(thiscoord
[0] - prevcoord
[0]) < larger
&&
1618 abs(thiscoord
[1] - prevcoord
[1]) < larger
&&
1619 abs(thiscoord
[2] - prevcoord
[2]) < larger
) {
1621 } else if (smallidx
> minidx
) {
1627 if (abs(thiscoord
[0] - thiscoord
[3]) < smallnum
&&
1628 abs(thiscoord
[1] - thiscoord
[4]) < smallnum
&&
1629 abs(thiscoord
[2] - thiscoord
[5]) < smallnum
) {
1630 /* interchange first with second atom for better
1631 * compression of water molecules
1633 tmp
= thiscoord
[0]; thiscoord
[0] = thiscoord
[3];
1635 tmp
= thiscoord
[1]; thiscoord
[1] = thiscoord
[4];
1637 tmp
= thiscoord
[2]; thiscoord
[2] = thiscoord
[5];
1642 tmpcoord
[0] = thiscoord
[0] - minint
[0];
1643 tmpcoord
[1] = thiscoord
[1] - minint
[1];
1644 tmpcoord
[2] = thiscoord
[2] - minint
[2];
1646 encodebits(buf2
, bitsizeint
[0], tmpcoord
[0]);
1647 encodebits(buf2
, bitsizeint
[1], tmpcoord
[1]);
1648 encodebits(buf2
, bitsizeint
[2], tmpcoord
[2]);
1650 encodeints(buf2
, 3, bitsize
, sizeint
, tmpcoord
);
1652 prevcoord
[0] = thiscoord
[0];
1653 prevcoord
[1] = thiscoord
[1];
1654 prevcoord
[2] = thiscoord
[2];
1655 thiscoord
= thiscoord
+ 3;
1659 if (is_small
== 0 && is_smaller
== -1)
1661 while (is_small
&& run
< 8*3) {
1664 tmp
=thiscoord
[j
] - prevcoord
[j
];
1667 if (is_smaller
== -1 && tmpsum
>= smaller
* smaller
) {
1671 tmpcoord
[run
++] = thiscoord
[0] - prevcoord
[0] + smallnum
;
1672 tmpcoord
[run
++] = thiscoord
[1] - prevcoord
[1] + smallnum
;
1673 tmpcoord
[run
++] = thiscoord
[2] - prevcoord
[2] + smallnum
;
1675 prevcoord
[0] = thiscoord
[0];
1676 prevcoord
[1] = thiscoord
[1];
1677 prevcoord
[2] = thiscoord
[2];
1680 thiscoord
= thiscoord
+ 3;
1683 abs(thiscoord
[0] - prevcoord
[0]) < smallnum
&&
1684 abs(thiscoord
[1] - prevcoord
[1]) < smallnum
&&
1685 abs(thiscoord
[2] - prevcoord
[2]) < smallnum
) {
1689 if (run
!= prevrun
|| is_smaller
!= 0) {
1691 encodebits(buf2
, 1, 1); /* flag the change in run-length */
1692 encodebits(buf2
, 5, run
+is_smaller
+1);
1694 encodebits(buf2
, 1, 0); /* flag the fact that runlength did not change */
1696 for (k
=0; k
< run
; k
+=3) {
1697 encodeints(buf2
, 3, smallidx
, sizesmall
, &tmpcoord
[k
]);
1699 if (is_smaller
!= 0) {
1700 smallidx
+= is_smaller
;
1701 if (is_smaller
< 0) {
1703 smaller
= magicints
[smallidx
-1] / 2;
1706 smallnum
= magicints
[smallidx
] / 2;
1708 sizesmall
[0] = sizesmall
[1] = sizesmall
[2] = magicints
[smallidx
];
1711 if (buf2
[1] != 0) buf2
[0]++;
1712 xdrfile_write_int(buf2
,1,xfp
); /* buf2[0] holds the length in bytes */
1713 tmp
=xdrfile_write_opaque((char *)&(buf2
[3]),(unsigned int)buf2
[0],xfp
);
1714 if(tmp
==(unsigned int)buf2
[0])
1721 /* Dont try do document Fortran interface, since
1722 * Doxygen barfs at the F77_FUNC macro
1726 /*************************************************************
1727 * Fortran77 interface for reading/writing portable data *
1728 * The routine are not threadsafe when called from Fortran *
1729 * (as they are when called from C) unless you compile with *
1730 * this file with posix thread support. *
1731 * Note that these are not multithread-safe. *
1732 *************************************************************/
1733 /*hms Moved to xdrfile.h to make them available across source files */
1734 /* #define MAX_FORTRAN_XDR 1024 */
1735 /* static XDRFILE *f77xdr[MAX_FORTRAN_XDR]; /\* array of file handles *\/ */
1736 XDRFILE
*f77xdr
[MAX_FORTRAN_XDR
]; /* array of file handles */
1737 static int f77init
= 1; /* zero array first time */
1739 /* internal to this file: C<-->Fortran string conversion */
1740 static int ftocstr(char *dest
, int dest_len
, char *src
, int src_len
);
1741 static int ctofstr(char *dest
, int dest_len
, char *src
);
1745 F77_FUNC(xdropen
,XDROPEN
)(int *fid
, char *filename
, char *mode
,
1746 int fn_len
, int mode_len
)
1748 char cfilename
[512];
1752 /* zero array at first invocation */
1754 for(i
=0;i
<MAX_FORTRAN_XDR
;i
++)
1760 /* f77xdr is always smaller or equal to MAX_FORTRAN_XDR */
1761 while(i
<MAX_FORTRAN_XDR
&& f77xdr
[i
]!=NULL
)
1763 if(i
==MAX_FORTRAN_XDR
) {
1765 } else if (ftocstr(cfilename
, sizeof(cfilename
), filename
, fn_len
)) {
1767 } else if (ftocstr(cmode
, sizeof(cmode
), mode
,mode_len
)) {
1770 f77xdr
[i
]=xdrfile_open(cfilename
,cmode
);
1771 /* return the index in the array as a fortran file handle */
1777 F77_FUNC(xdrclose
,XDRCLOSE
)(int *fid
)
1779 /* first close it */
1780 xdrfile_close(f77xdr
[*fid
]);
1781 /* the remove it from file handle list */
1787 F77_FUNC(xdrrint
,XDRRINT
)(int *fid
, int *data
, int *ndata
, int *ret
)
1789 *ret
= xdrfile_read_int(data
,*ndata
,f77xdr
[*fid
]);
1793 F77_FUNC(xdrwint
,XDRWINT
)(int *fid
, int *data
, int *ndata
, int *ret
)
1795 *ret
= xdrfile_write_int(data
,*ndata
,f77xdr
[*fid
]);
1799 F77_FUNC(xdrruint
,XDRRUINT
)(int *fid
, unsigned int *data
, int *ndata
, int *ret
)
1801 *ret
= xdrfile_read_uint(data
,*ndata
,f77xdr
[*fid
]);
1805 F77_FUNC(xdrwuint
,XDRWUINT
)(int *fid
, unsigned int *data
, int *ndata
, int *ret
)
1807 *ret
= xdrfile_write_uint(data
,*ndata
,f77xdr
[*fid
]);
1811 F77_FUNC(xdrrchar
,XDRRCHAR
)(int *fid
, char *ip
, int *ndata
, int *ret
)
1813 *ret
= xdrfile_read_char(ip
,*ndata
,f77xdr
[*fid
]);
1817 F77_FUNC(xdrwchar
,XDRWCHAR
)(int *fid
, char *ip
, int *ndata
, int *ret
)
1819 *ret
= xdrfile_write_char(ip
,*ndata
,f77xdr
[*fid
]);
1823 F77_FUNC(xdrruchar
,XDRRUCHAR
)(int *fid
, unsigned char *ip
, int *ndata
, int *ret
)
1825 *ret
= xdrfile_read_uchar(ip
,*ndata
,f77xdr
[*fid
]);
1829 F77_FUNC(xdrwuchar
,XDRWUCHAR
)(int *fid
, unsigned char *ip
, int *ndata
, int *ret
)
1831 *ret
= xdrfile_write_uchar(ip
,*ndata
,f77xdr
[*fid
]);
1835 F77_FUNC(xdrrshort
,XDRRSHORT
)(int *fid
, short *ip
, int *ndata
, int *ret
)
1837 *ret
= xdrfile_read_short(ip
,*ndata
,f77xdr
[*fid
]);
1841 F77_FUNC(xdrwshort
,XDRWSHORT
)(int *fid
, short *ip
, int *ndata
, int *ret
)
1843 *ret
= xdrfile_write_short(ip
,*ndata
,f77xdr
[*fid
]);
1847 F77_FUNC(xdrrushort
,XDRRUSHORT
)(int *fid
, unsigned short *ip
, int *ndata
, int *ret
)
1849 *ret
= xdrfile_read_ushort(ip
,*ndata
,f77xdr
[*fid
]);
1853 F77_FUNC(xdrwushort
,XDRWUSHORT
)(int *fid
, unsigned short *ip
, int *ndata
, int *ret
)
1855 *ret
= xdrfile_write_ushort(ip
,*ndata
,f77xdr
[*fid
]);
1859 F77_FUNC(xdrrsingle
,XDRRSINGLE
)(int *fid
, float *data
, int *ndata
, int *ret
)
1861 *ret
= xdrfile_read_float(data
,*ndata
,f77xdr
[*fid
]);
1865 F77_FUNC(xdrwsingle
,XDRWSINGLE
)(int *fid
, float *data
, int *ndata
, int *ret
)
1867 *ret
= xdrfile_write_float(data
,*ndata
,f77xdr
[*fid
]);
1871 F77_FUNC(xdrrdouble
,XDRRDOUBLE
)(int *fid
, double *data
, int *ndata
, int *ret
)
1873 *ret
= xdrfile_read_double(data
,*ndata
,f77xdr
[*fid
]);
1877 F77_FUNC(xdrwdouble
,XDRWDOUBLE
)(int *fid
, double *data
, int *ndata
, int *ret
)
1879 *ret
= xdrfile_write_double(data
,*ndata
,f77xdr
[*fid
]);
1882 static int ftocstr(char *dest
, int destlen
, char *src
, int srclen
)
1887 while ( --p
>= src
&& *p
== ' ' );
1888 srclen
= p
- src
+ 1;
1891 if (srclen
> destlen
)
1900 static int ctofstr(char *dest
, int destlen
, char *src
)
1902 while (destlen
&& *src
) {
1913 F77_FUNC(xdrrstring
,XDRRSTRING
)(int *fid
, char *str
, int *ret
, int len
)
1917 if((cstr
=(char*)malloc((len
+1)*sizeof(char)))==NULL
) {
1921 if (ftocstr(cstr
, len
+1, str
, len
)) {
1927 *ret
= xdrfile_read_string(cstr
, len
+1,f77xdr
[*fid
]);
1928 ctofstr( str
, len
, cstr
);
1933 F77_FUNC(xdrwstring
,XDRWSTRING
)(int *fid
, char *str
, int *ret
, int len
)
1937 if((cstr
=(char*)malloc((len
+1)*sizeof(char)))==NULL
) {
1941 if (ftocstr(cstr
, len
+1, str
, len
)) {
1947 *ret
= xdrfile_write_string(cstr
, f77xdr
[*fid
]);
1948 ctofstr( str
, len
, cstr
);
1953 F77_FUNC(xdrropaque
,XDRROPAQUE
)(int *fid
, char *data
, int *ndata
, int *ret
)
1955 *ret
= xdrfile_read_opaque(data
,*ndata
,f77xdr
[*fid
]);
1959 F77_FUNC(xdrwopaque
,XDRWOPAQUE
)(int *fid
, char *data
, int *ndata
, int *ret
)
1961 *ret
= xdrfile_write_opaque(data
,*ndata
,f77xdr
[*fid
]);
1965 /* Write single-precision compressed 3d coordinates */
1967 F77_FUNC(xdrccs
,XDRCCS
)(int *fid
, float *data
, int *ncoord
,
1968 float *precision
, int *ret
)
1970 *ret
= xdrfile_compress_coord_float(data
,*ncoord
,*precision
,f77xdr
[*fid
]);
1974 /* Read single-precision compressed 3d coordinates */
1976 F77_FUNC(xdrdcs
,XDRDCS
)(int *fid
, float *data
, int *ncoord
,
1977 float *precision
, int *ret
)
1979 *ret
= xdrfile_decompress_coord_float(data
,ncoord
,precision
,f77xdr
[*fid
]);
1983 /* Write compressed 3d coordinates from double precision data */
1985 F77_FUNC(xdrccd
,XDRCCD
)(int *fid
, double *data
, int *ncoord
,
1986 double *precision
, int *ret
)
1988 *ret
= xdrfile_compress_coord_double(data
,*ncoord
,*precision
,f77xdr
[*fid
]);
1991 /* Read compressed 3d coordinates into double precision data */
1993 F77_FUNC(xddcd
,XDRDCD
)(int *fid
, double *data
, int *ncoord
,
1994 double *precision
, int *ret
)
1996 *ret
= xdrfile_decompress_coord_double(data
,ncoord
,precision
,f77xdr
[*fid
]);
2005 #endif /* DOXYGEN */
2007 /*************************************************************
2008 * End of higher-level routines - dont change things below! *
2009 *************************************************************/
2030 /*************************************************************
2031 * The rest of this file contains our own implementation *
2032 * of the XDR calls in case you are compiling without them. *
2033 * You do NOT want to change things here since it would make *
2034 * things incompatible with the standard RPC/XDR routines. *
2035 *************************************************************/
2036 #ifndef HAVE_RPC_XDR_H
2039 * What follows is a modified version of the Sun XDR code. For reference
2040 * we include their copyright and license:
2042 * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
2043 * unrestricted use provided that this legend is included on all tape
2044 * media and as a part of the software program in whole or part. Users
2045 * may copy or modify Sun RPC without charge, but are not authorized
2046 * to license or distribute it to anyone else except as part of a product or
2047 * program developed by the user.
2049 * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
2050 * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
2051 * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
2053 * Sun RPC is provided with no support and without any obligation on the
2054 * part of Sun Microsystems, Inc. to assist in its use, correction,
2055 * modification or enhancement.
2057 * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
2058 * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
2059 * OR ANY PART THEREOF.
2061 * In no event will Sun Microsystems, Inc. be liable for any lost revenue
2062 * or profits or other special, indirect and consequential damages, even if
2063 * Sun has been advised of the possibility of such damages.
2065 * Sun Microsystems, Inc.
2066 * 2550 Garcia Avenue
2067 * Mountain View, California 94043
2070 /* INT_MAX is defined in limits.h according to ANSI C */
2071 #if (INT_MAX > 2147483647)
2072 # error Error: Cannot use builtin XDR support when size of int
2073 # error is larger than 4 bytes. Use your system XDR libraries
2074 # error instead, or modify the source code in xdrfile.c
2075 #endif /* Check for 4 byte int type */
2081 typedef int (*xdrproc_t
) (XDR
*, void *,...);
2083 #define xdr_getlong(xdrs, longp) \
2084 (*(xdrs)->x_ops->x_getlong)(xdrs, longp)
2085 #define xdr_putlong(xdrs, longp) \
2086 (*(xdrs)->x_ops->x_putlong)(xdrs, longp)
2087 #define xdr_getbytes(xdrs, addr, len) \
2088 (*(xdrs)->x_ops->x_getbytes)(xdrs, addr, len)
2089 #define xdr_putbytes(xdrs, addr, len) \
2090 (*(xdrs)->x_ops->x_putbytes)(xdrs, addr, len)
2092 #define BYTES_PER_XDR_UNIT 4
2093 static char xdr_zero
[BYTES_PER_XDR_UNIT
] = {0, 0, 0, 0};
2096 xdr_swapbytes(int32_t x
)
2099 char *px
=(char *)&x
;
2100 char *py
=(char *)&y
;
2109 xdr_htonl(int32_t x
)
2112 if( *((char *)&s
)==(char)0x34) {
2113 /* smallendian,swap bytes */
2114 return xdr_swapbytes(x
);
2116 /* bigendian, do nothing */
2125 if( *((char *)&s
)==(char)0x34) {
2126 /* smallendian, swap bytes */
2127 return xdr_swapbytes(x
);
2129 /* bigendian, do nothing */
2135 xdr_int (XDR
*xdrs
, int *ip
)
2142 i32
= (int32_t) *ip
;
2143 return xdr_putlong (xdrs
, &i32
);
2146 if (!xdr_getlong (xdrs
, &i32
))
2158 xdr_u_int (XDR
*xdrs
, unsigned int *up
)
2165 ui32
= (uint32_t) * up
;
2166 return xdr_putlong (xdrs
, (int32_t *)&ui32
);
2169 if (!xdr_getlong (xdrs
, (int32_t *)&ui32
))
2173 *up
= (uint32_t) ui32
;
2181 xdr_short (XDR
*xdrs
, short *sp
)
2188 i32
= (int32_t) *sp
;
2189 return xdr_putlong (xdrs
, &i32
);
2192 if (!xdr_getlong (xdrs
, &i32
))
2206 xdr_u_short (XDR
*xdrs
, unsigned short *sp
)
2213 ui32
= (uint32_t) *sp
;
2214 return xdr_putlong (xdrs
, (int32_t *)&ui32
);
2217 if (!xdr_getlong (xdrs
, (int32_t *)&ui32
))
2221 *sp
= (unsigned short) ui32
;
2231 xdr_char (XDR
*xdrs
, char *cp
)
2236 if (!xdr_int (xdrs
, &i
))
2245 xdr_u_char (XDR
*xdrs
, unsigned char *cp
)
2250 if (!xdr_u_int (xdrs
, &u
))
2260 * Allows the specification of a fixed size sequence of opaque bytes.
2261 * cp points to the opaque object and cnt gives the byte length.
2264 xdr_opaque (XDR
*xdrs
, char *cp
, unsigned int cnt
)
2267 static char crud
[BYTES_PER_XDR_UNIT
];
2270 * if no data we are done
2276 * round byte count to full xdr units
2278 rndup
= cnt
% BYTES_PER_XDR_UNIT
;
2280 rndup
= BYTES_PER_XDR_UNIT
- rndup
;
2285 if (!xdr_getbytes (xdrs
, cp
, cnt
))
2291 return xdr_getbytes (xdrs
, (char *)crud
, rndup
);
2294 if (!xdr_putbytes (xdrs
, cp
, cnt
))
2300 return xdr_putbytes (xdrs
, xdr_zero
, rndup
);
2305 #undef BYTES_PER_XDR_UNIT
2311 * XDR null terminated ASCII strings
2314 xdr_string (XDR
*xdrs
, char **cpp
, unsigned int maxsize
)
2316 char *sp
= *cpp
; /* sp is the actual string pointer */
2318 unsigned int nodesize
;
2321 * first deal with the length since xdr strings are counted-strings
2328 return 1; /* already free */
2330 /* fall through... */
2339 if (!xdr_u_int (xdrs
, &size
))
2347 nodesize
= size
+ 1;
2350 * now deal with the actual bytes
2360 *cpp
= sp
= (char *) malloc (nodesize
);
2363 (void) fputs ("xdr_string: out of memory\n", stderr
);
2370 return xdr_opaque (xdrs
, sp
, size
);
2382 /* Floating-point stuff */
2385 xdr_float(XDR
*xdrs
, float *fp
)
2387 switch (xdrs
->x_op
) {
2390 if (sizeof(float) == sizeof(int32_t))
2391 return (xdr_putlong(xdrs
, (int32_t *)fp
));
2392 else if (sizeof(float) == sizeof(int)) {
2393 int32_t tmp
= *(int *)fp
;
2394 return (xdr_putlong(xdrs
, &tmp
));
2399 if (sizeof(float) == sizeof(int32_t))
2400 return (xdr_getlong(xdrs
, (int32_t *)fp
));
2401 else if (sizeof(float) == sizeof(int)) {
2403 if (xdr_getlong(xdrs
, &tmp
)) {
2418 xdr_double(XDR
*xdrs
, double *dp
)
2420 /* Gromacs detects floating-point stuff at compile time, which is faster */
2422 # ifndef FLOAT_FORMAT_IEEE754
2423 # error non-IEEE floating point system, or you defined GROMACS yourself...
2426 # ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
2430 # endif /* Big endian word order */
2432 /* Outside Gromacs we rely on dynamic detection of FP order. */
2433 int LSW
; /* Least significant fp word */
2435 double x
=0.987654321; /* Just a number */
2436 unsigned char ix
= *((char *)&x
);
2438 /* Possible representations in IEEE double precision:
2439 * (S=small endian, B=big endian)
2441 * Byte order, Word order, Hex
2442 * S S b8 56 0e 3c dd 9a ef 3f
2443 * B S 3c 0e 56 b8 3f ef 9a dd
2444 * S B dd 9a ef 3f b8 56 0e 3c
2445 * B B 3f ef 9a dd 3c 0e 56 b8
2447 if(ix
==0xdd || ix
==0x3f)
2448 LSW
=1; /* Big endian word order */
2449 else if(ix
==0xb8 || ix
==0x3c)
2450 LSW
=0; /* Small endian word order */
2451 else { /* Catch strange errors */
2452 fprintf(stderr
,"Cannot detect floating-point word order.\n"
2453 "Do you have a non-IEEE system?\n"
2454 "Use system XDR libraries or fix xdr_double().\n");
2457 #endif /* end of dynamic detection of fp word order */
2459 switch (xdrs
->x_op
) {
2462 if (2*sizeof(int32_t) == sizeof(double)) {
2463 int32_t *lp
= (int32_t *)dp
;
2464 return (xdr_putlong(xdrs
, lp
+!LSW
) &&
2465 xdr_putlong(xdrs
, lp
+LSW
));
2466 } else if (2*sizeof(int) == sizeof(double)) {
2467 int *ip
= (int *)dp
;
2471 return (xdr_putlong(xdrs
, tmp
) &&
2472 xdr_putlong(xdrs
, tmp
+1));
2477 if (2*sizeof(int32_t) == sizeof(double)) {
2478 int32_t *lp
= (int32_t *)dp
;
2479 return (xdr_getlong(xdrs
, lp
+!LSW
) &&
2480 xdr_getlong(xdrs
, lp
+LSW
));
2481 } else if (2*sizeof(int) == sizeof(double)) {
2482 int *ip
= (int *)dp
;
2484 if (xdr_getlong(xdrs
, tmp
+!LSW
) &&
2485 xdr_getlong(xdrs
, tmp
+LSW
)) {
2500 static int xdrstdio_getlong (XDR
*, int32_t *);
2501 static int xdrstdio_putlong (XDR
*, int32_t *);
2502 static int xdrstdio_getbytes (XDR
*, char *, unsigned int);
2503 static int xdrstdio_putbytes (XDR
*, char *, unsigned int);
2504 static unsigned int xdrstdio_getpos (XDR
*);
2505 static int xdrstdio_setpos (XDR
*, unsigned int);
2506 static void xdrstdio_destroy (XDR
*);
2509 * Ops vector for stdio type XDR
2511 static const struct xdr_ops xdrstdio_ops
=
2513 xdrstdio_getlong
, /* deserialize a long int */
2514 xdrstdio_putlong
, /* serialize a long int */
2515 xdrstdio_getbytes
, /* deserialize counted bytes */
2516 xdrstdio_putbytes
, /* serialize counted bytes */
2517 xdrstdio_getpos
, /* get offset in the stream */
2518 xdrstdio_setpos
, /* set offset in the stream */
2519 xdrstdio_destroy
, /* destroy stream */
2523 * Initialize a stdio xdr stream.
2524 * Sets the xdr stream handle xdrs for use on the stream file.
2525 * Operation flag is set to op.
2528 xdrstdio_create (XDR
*xdrs
, FILE *file
, enum xdr_op op
)
2532 xdrs
->x_ops
= (struct xdr_ops
*) &xdrstdio_ops
;
2533 xdrs
->x_private
= (char *) file
;
2537 * Destroy a stdio xdr stream.
2538 * Cleans up the xdr stream handle xdrs previously set up by xdrstdio_create.
2541 xdrstdio_destroy (XDR
*xdrs
)
2543 (void) fflush ((FILE *) xdrs
->x_private
);
2544 /* xx should we close the file ?? */
2548 xdrstdio_getlong (XDR
*xdrs
, int32_t *lp
)
2552 if (fread ((char *) & mycopy
, 4, 1, (FILE *) xdrs
->x_private
) != 1)
2554 *lp
= (int32_t) xdr_ntohl (mycopy
);
2559 xdrstdio_putlong (XDR
*xdrs
, int32_t *lp
)
2561 int32_t mycopy
= xdr_htonl (*lp
);
2563 if (fwrite ((char *) lp
, 4, 1, (FILE *) xdrs
->x_private
) != 1)
2569 xdrstdio_getbytes (XDR
*xdrs
, char *addr
, unsigned int len
)
2571 if ((len
!= 0) && (fread (addr
, (int) len
, 1,
2572 (FILE *) xdrs
->x_private
) != 1))
2578 xdrstdio_putbytes (XDR
*xdrs
, char *addr
, unsigned int len
)
2580 if ((len
!= 0) && (fwrite (addr
, (int) len
, 1,
2581 (FILE *) xdrs
->x_private
) != 1))
2586 /* 32 bit fileseek operations */
2588 xdrstdio_getpos (XDR
*xdrs
)
2590 return (unsigned int) ftell ((FILE *) xdrs
->x_private
);
2594 xdrstdio_setpos (XDR
*xdrs
, unsigned int pos
)
2596 return fseek ((FILE *) xdrs
->x_private
, pos
, 0) < 0 ? 0 : 1;
2601 #endif /* HAVE_RPC_XDR_H not defined */