?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Two patch fo extra functionality in python and fortran.
[gromacs/libxdrfile.git] / src / xdrfile.c
blob18cd234626ac0b99b1d08c4c248ee1f3a8227d0f
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*-
2 *
3 * $Id$
4 *
5 * Copyright (c) Erik Lindahl, David van der Spoel 2003,2004.
6 * Coordinate compression (c) by Frans van Hoesel.
7 *
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.
12 */
13
14 /* Get HAVE_RPC_XDR_H, F77_FUNC from config.h if available */
15 #ifdef HAVE_CONFIG_H
16 #include <config.h>
17 #endif
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <math.h>
23 #include <limits.h>
24
25 #define _FILE_OFFSET_BITS 64
26
27 /* get fixed-width types if we are using ANSI C99 */
28 #ifdef HAVE_STDINT_H
29 # include <stdint.h>
30 #elif (defined HAVE_INTTYPES_H)
31 # include <inttypes.h>
32 #endif
33
34
35 #ifdef HAVE_RPC_XDR_H
36 # include <rpc/rpc.h>
37 # include <rpc/xdr.h>
38 #endif
39
40 #include "xdrfile.h"
41
42 /* Default FORTRAN name mangling is: lower case name, append underscore */
43 #ifndef F77_FUNC
44 #define F77_FUNC(name,NAME) name ## _
45 #endif
46
47 char *exdr_message[exdrNR] = {
48 "OK",
49 "Header",
50 "String",
51 "Double",
52 "Integer",
53 "Float",
54 "Unsigned integer",
55 "Compressed 3D coordinate",
56 "Closing file",
57 "Magic number",
58 "Not enough memory",
59 "End of file",
60 "File not found"
61 };
62
63 /*
64 * Declare our own XDR routines statically if no libraries are present.
65 * Actual implementation is at the end of this file.
66 *
67 * We don't want the low-level XDR implementation as part of the Gromacs
68 * documentation, so skip it for doxygen too...
69 */
70 #if (!defined HAVE_RPC_XDR_H && !defined DOXYGEN)
71
72 enum xdr_op
73 {
74 XDR_ENCODE = 0,
75 XDR_DECODE = 1,
76 XDR_FREE = 2
77 };
78
79
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.
85 */
86 #if !(defined INT32_MAX || defined DOXYGEN)
87 # if (INT_MAX == 2147483647)
88 # define int32_t int
89 # define uint32_t unsigned int
90 # define INT32_MAX 2147483647
91 # elif (LONG_MAX == 2147483647)
92 # define int32_t long
93 # define uint32_t unsigned long
94 # define INT32_MAX 2147483647L
95 # else
96 # error ERROR: No 32 bit wide integer type found!
97 # error Use system XDR libraries instead, or update xdrfile.c
98 # endif
99 #endif
100
101 typedef struct XDR XDR;
102
103 struct XDR
104 {
105 enum xdr_op x_op;
106 struct xdr_ops
107 {
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);
116 }
117 *x_ops;
118 char *x_private;
119 };
120
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);
132
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) \
138 do { \
139 if ((xdrs)->x_ops->x_destroy) \
140 (*(xdrs)->x_ops->x_destroy)(xdrs); \
141 } while (0)
142 #endif /* end of our own XDR declarations */
143
144
145
146
147
148 /** Contents of the abstract XDRFILE data structure.
149 *
150 * @internal
151 *
152 * This structure is used to provide an XDR file interface that is
153 * virtual identical to the standard UNIX fopen/fread/fwrite/fclose.
154 */
155 struct XDRFILE
156 {
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 */
164 };
165
166
167
168
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 *************************************************************/
174
175 XDRFILE *
176 xdrfile_open(const char *path, const char *mode)
177 {
178 char newmode[5];
179 enum xdr_op xdrmode;
180 XDRFILE *xfp;
181
182 /* make sure XDR files are opened in binary mode... */
183 if(*mode=='w' || *mode=='W')
184 {
185 sprintf(newmode,"wb+");
186 xdrmode=XDR_ENCODE;
187 } else if(*mode == 'a' || *mode == 'A')
188 {
189 sprintf(newmode,"ab+");
190 xdrmode = XDR_ENCODE;
191 } else if(*mode == 'r' || *mode == 'R')
192 {
193 sprintf(newmode,"rb");
194 xdrmode = XDR_DECODE;
195 } else /* cannot determine mode */
196 return NULL;
197
198 if((xfp=(XDRFILE *)malloc(sizeof(XDRFILE)))==NULL)
199 return NULL;
200 if((xfp->fp=fopen(path,newmode))==NULL)
201 {
202 free(xfp);
203 return NULL;
204 }
205 if((xfp->xdr=(XDR *)malloc(sizeof(XDR)))==NULL)
206 {
207 fclose(xfp->fp);
208 free(xfp);
209 return NULL;
210 }
211 xfp->mode=*mode;
212 xdrstdio_create((XDR *)(xfp->xdr),xfp->fp,xdrmode);
213 xfp->buf1 = xfp->buf2 = NULL;
214 xfp->buf1size = xfp->buf2size = 0;
215 return xfp;
216 }
217
218 int
219 xdrfile_close(XDRFILE *xfp)
220 {
221 int ret=exdrCLOSE;
222 if(xfp)
223 {
224 /* flush and destroy XDR stream */
225 if(xfp->xdr)
226 xdr_destroy((XDR *)(xfp->xdr));
227 free(xfp->xdr);
228 /* close the file */
229 ret=fclose(xfp->fp);
230 if(xfp->buf1size)
231 free(xfp->buf1);
232 if(xfp->buf2size)
233 free(xfp->buf2);
234 free(xfp);
235 }
236 return ret; /* return 0 if ok */
237 }
238
239
240
241 int
242 xdrfile_read_int(int *ptr, int ndata, XDRFILE* xfp)
243 {
244 int i=0;
245
246 /* read write is encoded in the XDR struct */
247 while(i<ndata && xdr_int((XDR *)(xfp->xdr),ptr+i))
248 i++;
249
250 return i;
251 }
252
253 int
254 xdrfile_write_int(int *ptr, int ndata, XDRFILE* xfp)
255 {
256 int i=0;
257
258 /* read write is encoded in the XDR struct */
259 while(i<ndata && xdr_int((XDR *)(xfp->xdr),ptr+i))
260 i++;
261 return i;
262 }
263
264
265 int
266 xdrfile_read_uint(unsigned int *ptr, int ndata, XDRFILE* xfp)
267 {
268 int i=0;
269
270 /* read write is encoded in the XDR struct */
271 while(i<ndata && xdr_u_int((XDR *)(xfp->xdr),ptr+i))
272 i++;
273
274 return i;
275 }
276
277 int
278 xdrfile_write_uint(unsigned int *ptr, int ndata, XDRFILE* xfp)
279 {
280 int i=0;
281
282 /* read write is encoded in the XDR struct */
283 while(i<ndata && xdr_u_int((XDR *)(xfp->xdr),ptr+i))
284 i++;
285 return i;
286 }
287
288 int
289 xdrfile_read_char(char *ptr, int ndata, XDRFILE* xfp)
290 {
291 int i=0;
292
293 /* read write is encoded in the XDR struct */
294 while(i<ndata && xdr_char((XDR *)(xfp->xdr),ptr+i))
295 i++;
296
297 return i;
298 }
299
300 int
301 xdrfile_write_char(char *ptr, int ndata, XDRFILE* xfp)
302 {
303 int i=0;
304
305 /* read write is encoded in the XDR struct */
306 while(i<ndata && xdr_char((XDR *)(xfp->xdr),ptr+i))
307 i++;
308 return i;
309 }
310
311
312 int
313 xdrfile_read_uchar(unsigned char *ptr, int ndata, XDRFILE* xfp)
314 {
315 int i=0;
316
317 /* read write is encoded in the XDR struct */
318 while(i<ndata && xdr_u_char((XDR *)(xfp->xdr),ptr+i))
319 i++;
320
321 return i;
322 }
323
324 int
325 xdrfile_write_uchar(unsigned char *ptr, int ndata, XDRFILE* xfp)
326 {
327 int i=0;
328
329 /* read write is encoded in the XDR struct */
330 while(i<ndata && xdr_u_char((XDR *)(xfp->xdr),ptr+i))
331 i++;
332 return i;
333 }
334
335 int
336 xdrfile_read_short(short *ptr, int ndata, XDRFILE* xfp)
337 {
338 int i=0;
339
340 /* read write is encoded in the XDR struct */
341 while(i<ndata && xdr_short((XDR *)(xfp->xdr),ptr+i))
342 i++;
343
344 return i;
345 }
346
347 int
348 xdrfile_write_short(short *ptr, int ndata, XDRFILE* xfp)
349 {
350 int i=0;
351
352 /* read write is encoded in the XDR struct */
353 while(i<ndata && xdr_short((XDR *)(xfp->xdr),ptr+i))
354 i++;
355 return i;
356 }
357
358
359 int
360 xdrfile_read_ushort(unsigned short *ptr, int ndata, XDRFILE* xfp)
361 {
362 int i=0;
363
364 /* read write is encoded in the XDR struct */
365 while(i<ndata && xdr_u_short((XDR *)(xfp->xdr),ptr+i))
366 i++;
367
368 return i;
369 }
370
371 int
372 xdrfile_write_ushort(unsigned short *ptr, int ndata, XDRFILE* xfp)
373 {
374 int i=0;
375
376 /* read write is encoded in the XDR struct */
377 while(i<ndata && xdr_u_short((XDR *)(xfp->xdr),ptr+i))
378 i++;
379 return i;
380 }
381
382 int
383 xdrfile_read_float(float *ptr, int ndata, XDRFILE* xfp)
384 {
385 int i=0;
386 /* read write is encoded in the XDR struct */
387 while(i<ndata && xdr_float((XDR *)(xfp->xdr),ptr+i))
388 i++;
389 return i;
390 }
391
392 int
393 xdrfile_write_float(float *ptr, int ndata, XDRFILE* xfp)
394 {
395 int i=0;
396 /* read write is encoded in the XDR struct */
397 while(i<ndata && xdr_float((XDR *)(xfp->xdr),ptr+i))
398 i++;
399 return i;
400 }
401
402 int
403 xdrfile_read_double(double *ptr, int ndata, XDRFILE* xfp)
404 {
405 int i=0;
406 /* read write is encoded in the XDR struct */
407 while(i<ndata && xdr_double((XDR *)(xfp->xdr),ptr+i))
408 i++;
409 return i;
410 }
411
412 int
413 xdrfile_write_double(double *ptr, int ndata, XDRFILE* xfp)
414 {
415 int i=0;
416 /* read write is encoded in the XDR struct */
417 while(i<ndata && xdr_double((XDR *)(xfp->xdr),ptr+i))
418 i++;
419 return i;
420 }
421
422 int
423 xdrfile_read_string(char *ptr, int maxlen, XDRFILE* xfp)
424 {
425 int i;
426 if(xdr_string((XDR *)(xfp->xdr),&ptr,maxlen)) {
427 i=0;
428 while(i<maxlen && ptr[i]!=0)
429 i++;
430 if(i==maxlen)
431 return maxlen;
432 else
433 return i+1;
434 } else
435 return 0;
436 }
437
438 int
439 xdrfile_write_string(char *ptr, XDRFILE* xfp)
440 {
441 int len=strlen(ptr)+1;
442
443 if(xdr_string((XDR *)(xfp->xdr),&ptr,len))
444 return len;
445 else
446 return 0;
447 }
448
449
450 int
451 xdrfile_read_opaque(char *ptr, int cnt, XDRFILE* xfp)
452 {
453 if(xdr_opaque((XDR *)(xfp->xdr),ptr,cnt))
454 return cnt;
455 else
456 return 0;
457 }
458
459
460 int
461 xdrfile_write_opaque(char *ptr, int cnt, XDRFILE* xfp)
462 {
463 if(xdr_opaque((XDR *)(xfp->xdr),ptr,cnt))
464 return cnt;
465 else
466 return 0;
467 }
468
469
470 /* Internal support routines for reading/writing compressed coordinates
471 * sizeofint - calculate smallest number of bits necessary
472 * to represent a certain integer.
473 */
474 static int
475 sizeofint(int size) {
476 unsigned int num = 1;
477 int num_of_bits = 0;
478
479 while (size >= num && num_of_bits < 32)
480 {
481 num_of_bits++;
482 num <<= 1;
483 }
484 return num_of_bits;
485 }
486
487
488 /*
489 * sizeofints - calculate 'bitsize' of compressed ints
490 *
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
497 * this routine).
498 */
499 static int
500 sizeofints(int num_of_ints, unsigned int sizes[])
501 {
502 int i, num;
503 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
504 num_of_bytes = 1;
505 bytes[0] = 1;
506 num_of_bits = 0;
507 for (i=0; i < num_of_ints; i++)
508 {
509 tmp = 0;
510 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
511 {
512 tmp = bytes[bytecnt] * sizes[i] + tmp;
513 bytes[bytecnt] = tmp & 0xff;
514 tmp >>= 8;
515 }
516 while (tmp != 0)
517 {
518 bytes[bytecnt++] = tmp & 0xff;
519 tmp >>= 8;
520 }
521 num_of_bytes = bytecnt;
522 }
523 num = 1;
524 num_of_bytes--;
525 while (bytes[num_of_bytes] >= num)
526 {
527 num_of_bits++;
528 num *= 2;
529 }
530 return num_of_bits + num_of_bytes * 8;
531
532 }
533
534
535 /*
536 * encodebits - encode num into buf using the specified number of bits
537 *
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.
542 */
543 static void
544 encodebits(int buf[], int num_of_bits, int num)
545 {
546
547 unsigned int cnt, lastbyte;
548 int lastbits;
549 unsigned char * cbuf;
550
551 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
552 cnt = (unsigned int) buf[0];
553 lastbits = buf[1];
554 lastbyte =(unsigned int) buf[2];
555 while (num_of_bits >= 8)
556 {
557 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
558 cbuf[cnt++] = lastbyte >> lastbits;
559 num_of_bits -= 8;
560 }
561 if (num_of_bits > 0)
562 {
563 lastbyte = (lastbyte << num_of_bits) | num;
564 lastbits += num_of_bits;
565 if (lastbits >= 8)
566 {
567 lastbits -= 8;
568 cbuf[cnt++] = lastbyte >> lastbits;
569 }
570 }
571 buf[0] = cnt;
572 buf[1] = lastbits;
573 buf[2] = lastbyte;
574 if (lastbits>0)
575 {
576 cbuf[cnt] = lastbyte << (8 - lastbits);
577 }
578 }
579
580 /*
581 * encodeints - encode a small set of small integers in compressed format
582 *
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...
593 */
594
595 static void
596 encodeints(int buf[], int num_of_ints, int num_of_bits,
597 unsigned int sizes[], unsigned int nums[])
598 {
599
600 int i;
601 unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
602
603 tmp = nums[0];
604 num_of_bytes = 0;
605 do
606 {
607 bytes[num_of_bytes++] = tmp & 0xff;
608 tmp >>= 8;
609 } while (tmp != 0);
610
611 for (i = 1; i < num_of_ints; i++)
612 {
613 if (nums[i] >= sizes[i])
614 {
615 fprintf(stderr,"major breakdown in encodeints - num %u doesn't "
616 "match size %u\n", nums[i], sizes[i]);
617 abort();
618 }
619 /* use one step multiply */
620 tmp = nums[i];
621 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
622 {
623 tmp = bytes[bytecnt] * sizes[i] + tmp;
624 bytes[bytecnt] = tmp & 0xff;
625 tmp >>= 8;
626 }
627 while (tmp != 0)
628 {
629 bytes[bytecnt++] = tmp & 0xff;
630 tmp >>= 8;
631 }
632 num_of_bytes = bytecnt;
633 }
634 if (num_of_bits >= num_of_bytes * 8)
635 {
636 for (i = 0; i < num_of_bytes; i++)
637 {
638 encodebits(buf, 8, bytes[i]);
639 }
640 encodebits(buf, num_of_bits - num_of_bytes * 8, 0);
641 }
642 else
643 {
644 for (i = 0; i < num_of_bytes-1; i++)
645 {
646 encodebits(buf, 8, bytes[i]);
647 }
648 encodebits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
649 }
650 }
651
652
653 /*
654 * decodebits - decode number from buf using specified number of bits
655 *
656 * extract the number of bits from the array buf and construct an integer
657 * from it. Return that value.
658 *
659 */
660
661 static int
662 decodebits(int buf[], int num_of_bits)
663 {
664
665 int cnt, num;
666 unsigned int lastbits, lastbyte;
667 unsigned char * cbuf;
668 int mask = (1 << num_of_bits) -1;
669
670 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
671 cnt = buf[0];
672 lastbits = (unsigned int) buf[1];
673 lastbyte = (unsigned int) buf[2];
674
675 num = 0;
676 while (num_of_bits >= 8)
677 {
678 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
679 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
680 num_of_bits -=8;
681 }
682 if (num_of_bits > 0)
683 {
684 if (lastbits < num_of_bits)
685 {
686 lastbits += 8;
687 lastbyte = (lastbyte << 8) | cbuf[cnt++];
688 }
689 lastbits -= num_of_bits;
690 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
691 }
692 num &= mask;
693 buf[0] = cnt;
694 buf[1] = lastbits;
695 buf[2] = lastbyte;
696 return num;
697 }
698
699 /*
700 * decodeints - decode 'small' integers from the buf array
701 *
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.
706 *
707 */
708
709 static void
710 decodeints(int buf[], int num_of_ints, int num_of_bits,
711 unsigned int sizes[], int nums[])
712 {
713
714 int bytes[32];
715 int i, j, num_of_bytes, p, num;
716
717 bytes[1] = bytes[2] = bytes[3] = 0;
718 num_of_bytes = 0;
719 while (num_of_bits > 8)
720 {
721 bytes[num_of_bytes++] = decodebits(buf, 8);
722 num_of_bits -= 8;
723 }
724 if (num_of_bits > 0)
725 {
726 bytes[num_of_bytes++] = decodebits(buf, num_of_bits);
727 }
728 for (i = num_of_ints-1; i > 0; i--)
729 {
730 num = 0;
731 for (j = num_of_bytes-1; j >=0; j--)
732 {
733 num = (num << 8) | bytes[j];
734 p = num / sizes[i];
735 bytes[j] = p;
736 num = num - p * sizes[i];
737 }
738 nums[i] = num;
739 }
740 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
741 }
742
743
744 static const int magicints[] =
745 {
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
753 };
754
755 #define FIRSTIDX 9
756 /* note that magicints[FIRSTIDX-1] == 0 */
757 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
758
759 /* Compressed coordinate routines - modified from the original
760 * implementation by Frans v. Hoesel to make them threadsafe.
761 */
762 int
763 xdrfile_decompress_coord_float(float *ptr,
764 int *size,
765 float *precision,
766 XDRFILE* xfp)
767 {
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;
776
777 bitsizeint[0] = 0;
778 bitsizeint[1] = 0;
779 bitsizeint[2] = 0;
780
781 if(xfp==NULL || ptr==NULL)
782 return -1;
783 tmp=xdrfile_read_int(&lsize,1,xfp);
784 if(tmp==0)
785 return -1; /* return if we could not read size */
786 if (*size < lsize)
787 {
788 fprintf(stderr, "Requested to decompress %d coords, file contains %d\n",
789 *size, lsize);
790 return -1;
791 }
792 *size = lsize;
793 size3 = *size * 3;
794 if(size3>xfp->buf1size)
795 {
796 if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL)
797 {
798 fprintf(stderr,"Cannot allocate memory for decompressing coordinates.\n");
799 return -1;
800 }
801 xfp->buf1size=size3;
802 xfp->buf2size=size3*1.2;
803 if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL)
804 {
805 fprintf(stderr,"Cannot allocate memory for decompressing coordinates.\n");
806 return -1;
807 }
808 }
809 /* Dont bother with compression for three atoms or less */
810 if(*size<=9)
811 {
812 return xdrfile_read_float(ptr,size3,xfp)/3;
813 /* return number of coords, not floats */
814 }
815 /* Compression-time if we got here. Read precision first */
816 xdrfile_read_float(precision,1,xfp);
817
818 /* avoid repeated pointer dereferencing. */
819 buf1=xfp->buf1;
820 buf2=xfp->buf2;
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);
825
826 sizeint[0] = maxint[0] - minint[0]+1;
827 sizeint[1] = maxint[1] - minint[1]+1;
828 sizeint[2] = maxint[2] - minint[2]+1;
829
830 /* check if one of the sizes is to big to be multiplied */
831 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
832 {
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 */
837 }
838 else
839 {
840 bitsize = sizeofints(3, sizeint);
841 }
842
843 if (xdrfile_read_int(&smallidx,1,xfp) == 0)
844 return 0; /* not sure what has happened here or why we return... */
845 tmp=smallidx+8;
846 maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
847 minidx = maxidx - 8; /* often this equal smallidx */
848 tmp = smallidx-1;
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];
854
855 /* buf2[0] holds the length in bytes */
856
857 if (xdrfile_read_int(buf2,1,xfp) == 0)
858 return 0;
859 if (xdrfile_read_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp) == 0)
860 return 0;
861 buf2[0] = buf2[1] = buf2[2] = 0;
862
863 lfp = ptr;
864 inv_precision = 1.0 / * precision;
865 run = 0;
866 i = 0;
867 lip = buf1;
868 while ( i < lsize )
869 {
870 thiscoord = (int *)(lip) + i * 3;
871
872 if (bitsize == 0)
873 {
874 thiscoord[0] = decodebits(buf2, bitsizeint[0]);
875 thiscoord[1] = decodebits(buf2, bitsizeint[1]);
876 thiscoord[2] = decodebits(buf2, bitsizeint[2]);
877 }
878 else
879 {
880 decodeints(buf2, 3, bitsize, sizeint, thiscoord);
881 }
882
883 i++;
884 thiscoord[0] += minint[0];
885 thiscoord[1] += minint[1];
886 thiscoord[2] += minint[2];
887
888 prevcoord[0] = thiscoord[0];
889 prevcoord[1] = thiscoord[1];
890 prevcoord[2] = thiscoord[2];
891
892 flag = decodebits(buf2, 1);
893 is_smaller = 0;
894 if (flag == 1)
895 {
896 run = decodebits(buf2, 5);
897 is_smaller = run % 3;
898 run -= is_smaller;
899 is_smaller--;
900 }
901 if (run > 0)
902 {
903 thiscoord += 3;
904 for (k = 0; k < run; k+=3)
905 {
906 decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
907 i++;
908 thiscoord[0] += prevcoord[0] - smallnum;
909 thiscoord[1] += prevcoord[1] - smallnum;
910 thiscoord[2] += prevcoord[2] - smallnum;
911 if (k == 0) {
912 /* interchange first with second atom for better
913 * compression of water molecules
914 */
915 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
916 prevcoord[0] = tmp;
917 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
918 prevcoord[1] = tmp;
919 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
920 prevcoord[2] = tmp;
921 *lfp++ = prevcoord[0] * inv_precision;
922 *lfp++ = prevcoord[1] * inv_precision;
923 *lfp++ = prevcoord[2] * inv_precision;
924 } else {
925 prevcoord[0] = thiscoord[0];
926 prevcoord[1] = thiscoord[1];
927 prevcoord[2] = thiscoord[2];
928 }
929 *lfp++ = thiscoord[0] * inv_precision;
930 *lfp++ = thiscoord[1] * inv_precision;
931 *lfp++ = thiscoord[2] * inv_precision;
932 }
933 }
934 else
935 {
936 *lfp++ = thiscoord[0] * inv_precision;
937 *lfp++ = thiscoord[1] * inv_precision;
938 *lfp++ = thiscoord[2] * inv_precision;
939 }
940 smallidx += is_smaller;
941 if (is_smaller < 0)
942 {
943 smallnum = smaller;
944
945 if (smallidx > FIRSTIDX)
946 {
947 smaller = magicints[smallidx - 1] /2;
948 }
949 else
950 {
951 smaller = 0;
952 }
953 }
954 else if (is_smaller > 0)
955 {
956 smaller = smallnum;
957 smallnum = magicints[smallidx] / 2;
958 }
959 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
960 }
961 return *size;
962 }
963
964 int
965 xdrfile_compress_coord_float(float *ptr,
966 int size,
967 float precision,
968 XDRFILE* xfp)
969 {
970 int minint[3], maxint[3], mindiff, *lip, diff;
971 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
972 int minidx, maxidx;
973 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
974 int k, *buf1, *buf2;
975 int smallnum, smaller, larger, i, j, is_small, is_smaller, run, prevrun;
976 float *lfp, lf;
977 int tmp, tmpsum, *thiscoord, prevcoord[3];
978 unsigned int tmpcoord[30];
979 int errval=1;
980 unsigned int bitsize;
981
982 if(xfp==NULL)
983 return -1;
984 size3=3*size;
985
986 bitsizeint[0] = 0;
987 bitsizeint[1] = 0;
988 bitsizeint[2] = 0;
989
990 if(size3>xfp->buf1size)
991 {
992 if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL)
993 {
994 fprintf(stderr,"Cannot allocate memory for compressing coordinates.\n");
995 return -1;
996 }
997 xfp->buf1size=size3;
998 xfp->buf2size=size3*1.2;
999 if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL)
1000 {
1001 fprintf(stderr,"Cannot allocate memory for compressing coordinates.\n");
1002 return -1;
1003 }
1004 }
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 */
1008 if(size<=9)
1009 {
1010 return xdrfile_write_float(ptr,size3,xfp)/3;
1011 /* return number of coords, not floats */
1012 }
1013 /* Compression-time if we got here. Write precision first */
1014 if (precision <= 0)
1015 precision = 1000;
1016 xdrfile_write_float(&precision,1,xfp);
1017 /* avoid repeated pointer dereferencing. */
1018 buf1=xfp->buf1;
1019 buf2=xfp->buf2;
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;
1024 prevrun = -1;
1025 lfp = ptr;
1026 lip = buf1;
1027 mindiff = INT_MAX;
1028 oldlint1 = oldlint2 = oldlint3 = 0;
1029 while(lfp < ptr + size3 )
1030 {
1031 /* find nearest integer */
1032 if (*lfp >= 0.0)
1033 lf = *lfp * precision + 0.5;
1034 else
1035 lf = *lfp * precision - 0.5;
1036 if (fabs(lf) > INT_MAX-2)
1037 {
1038 /* scaling would cause overflow */
1039 fprintf(stderr,"Internal overflow compressing coordinates.\n");
1040 errval=0;
1041 }
1042 lint1 = lf;
1043 if (lint1 < minint[0]) minint[0] = lint1;
1044 if (lint1 > maxint[0]) maxint[0] = lint1;
1045 *lip++ = lint1;
1046 lfp++;
1047 if (*lfp >= 0.0)
1048 lf = *lfp * precision + 0.5;
1049 else
1050 lf = *lfp * precision - 0.5;
1051 if (fabs(lf) > INT_MAX-2)
1052 {
1053 /* scaling would cause overflow */
1054 fprintf(stderr,"Internal overflow compressing coordinates.\n");
1055 errval=0;
1056 }
1057 lint2 = lf;
1058 if (lint2 < minint[1]) minint[1] = lint2;
1059 if (lint2 > maxint[1]) maxint[1] = lint2;
1060 *lip++ = lint2;
1061 lfp++;
1062 if (*lfp >= 0.0)
1063 lf = *lfp * precision + 0.5;
1064 else
1065 lf = *lfp * precision - 0.5;
1066 if (fabs(lf) > INT_MAX-2)
1067 {
1068 errval=0;
1069 }
1070 lint3 = lf;
1071 if (lint3 < minint[2]) minint[2] = lint3;
1072 if (lint3 > maxint[2]) maxint[2] = lint3;
1073 *lip++ = lint3;
1074 lfp++;
1075 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
1076 if (diff < mindiff && lfp > ptr + 3)
1077 mindiff = diff;
1078 oldlint1 = lint1;
1079 oldlint2 = lint2;
1080 oldlint3 = lint3;
1081 }
1082 xdrfile_write_int(minint,3,xfp);
1083 xdrfile_write_int(maxint,3,xfp);
1084
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
1090 */
1091 fprintf(stderr,"Internal overflow compressing coordinates.\n");
1092 errval=0;
1093 }
1094 sizeint[0] = maxint[0] - minint[0]+1;
1095 sizeint[1] = maxint[1] - minint[1]+1;
1096 sizeint[2] = maxint[2] - minint[2]+1;
1097
1098 /* check if one of the sizes is to big to be multiplied */
1099 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
1100 {
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 */
1105 }
1106 else
1107 {
1108 bitsize = sizeofints(3, sizeint);
1109 }
1110 lip = buf1;
1111 luip = (unsigned int *) buf1;
1112 smallidx = FIRSTIDX;
1113 while (smallidx < LASTIDX && magicints[smallidx] < mindiff)
1114 {
1115 smallidx++;
1116 }
1117 xdrfile_write_int(&smallidx,1,xfp);
1118 tmp=smallidx+8;
1119 maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
1120 minidx = maxidx - 8; /* often this equal smallidx */
1121 tmp=smallidx-1;
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;
1128 i = 0;
1129 while (i < size)
1130 {
1131 is_small = 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) {
1137 is_smaller = 1;
1138 }
1139 else if (smallidx > minidx)
1140 {
1141 is_smaller = -1;
1142 }
1143 else
1144 {
1145 is_smaller = 0;
1146 }
1147 if (i + 1 < size)
1148 {
1149 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
1150 abs(thiscoord[1] - thiscoord[4]) < smallnum &&
1151 abs(thiscoord[2] - thiscoord[5]) < smallnum)
1152 {
1153 /* interchange first with second atom for better
1154 * compression of water molecules
1155 */
1156 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
1157 thiscoord[3] = tmp;
1158 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
1159 thiscoord[4] = tmp;
1160 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
1161 thiscoord[5] = tmp;
1162 is_small = 1;
1163 }
1164 }
1165 tmpcoord[0] = thiscoord[0] - minint[0];
1166 tmpcoord[1] = thiscoord[1] - minint[1];
1167 tmpcoord[2] = thiscoord[2] - minint[2];
1168 if (bitsize == 0)
1169 {
1170 encodebits(buf2, bitsizeint[0], tmpcoord[0]);
1171 encodebits(buf2, bitsizeint[1], tmpcoord[1]);
1172 encodebits(buf2, bitsizeint[2], tmpcoord[2]);
1173 }
1174 else
1175 {
1176 encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
1177 }
1178 prevcoord[0] = thiscoord[0];
1179 prevcoord[1] = thiscoord[1];
1180 prevcoord[2] = thiscoord[2];
1181 thiscoord = thiscoord + 3;
1182 i++;
1183
1184 run = 0;
1185 if (is_small == 0 && is_smaller == -1)
1186 is_smaller = 0;
1187 while (is_small && run < 8*3)
1188 {
1189 tmpsum=0;
1190 for(j=0;j<3;j++)
1191 {
1192 tmp=thiscoord[j] - prevcoord[j];
1193 tmpsum+=tmp*tmp;
1194 }
1195 if (is_smaller == -1 && tmpsum >= smaller * smaller)
1196 {
1197 is_smaller = 0;
1198 }
1199
1200 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
1201 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
1202 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
1203
1204 prevcoord[0] = thiscoord[0];
1205 prevcoord[1] = thiscoord[1];
1206 prevcoord[2] = thiscoord[2];
1207
1208 i++;
1209 thiscoord = thiscoord + 3;
1210 is_small = 0;
1211 if (i < size &&
1212 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
1213 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
1214 abs(thiscoord[2] - prevcoord[2]) < smallnum)
1215 {
1216 is_small = 1;
1217 }
1218 }
1219 if (run != prevrun || is_smaller != 0)
1220 {
1221 prevrun = run;
1222 encodebits(buf2, 1, 1); /* flag the change in run-length */
1223 encodebits(buf2, 5, run+is_smaller+1);
1224 }
1225 else
1226 {
1227 encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
1228 }
1229 for (k=0; k < run; k+=3)
1230 {
1231 encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
1232 }
1233 if (is_smaller != 0)
1234 {
1235 smallidx += is_smaller;
1236 if (is_smaller < 0)
1237 {
1238 smallnum = smaller;
1239 smaller = magicints[smallidx-1] / 2;
1240 }
1241 else
1242 {
1243 smaller = smallnum;
1244 smallnum = magicints[smallidx] / 2;
1245 }
1246 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1247 }
1248 }
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])
1253 return size;
1254 else
1255 return -1;
1256 }
1257
1258
1259 int
1260 xdrfile_decompress_coord_double(double *ptr,
1261 int *size,
1262 double *precision,
1263 XDRFILE* xfp)
1264 {
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;
1274
1275 bitsizeint[0] = 0;
1276 bitsizeint[1] = 0;
1277 bitsizeint[2] = 0;
1278
1279 if(xfp==NULL || ptr==NULL)
1280 return -1;
1281 tmp=xdrfile_read_int(&lsize,1,xfp);
1282 if(tmp==0)
1283 return -1; /* return if we could not read size */
1284 if (*size < lsize)
1285 {
1286 fprintf(stderr, "Requested to decompress %d coords, file contains %d\n",
1287 *size, lsize);
1288 return -1;
1289 }
1290 *size = lsize;
1291 size3 = *size * 3;
1292 if(size3>xfp->buf1size)
1293 {
1294 if((xfp->buf1=(int *)malloc(sizeof(int)*size3))==NULL)
1295 {
1296 fprintf(stderr,"Cannot allocate memory for decompression coordinates.\n");
1297 return -1;
1298 }
1299 xfp->buf1size=size3;
1300 xfp->buf2size=size3*1.2;
1301 if((xfp->buf2=(int *)malloc(sizeof(int)*xfp->buf2size))==NULL)
1302 {
1303 fprintf(stderr,"Cannot allocate memory for decompressing coordinates.\n");
1304 return -1;
1305 }
1306 }
1307 /* Dont bother with compression for three atoms or less */
1308 if(*size<=9)
1309 {
1310 tmp=xdrfile_read_float(tmpdata,size3,xfp);
1311 for(i=0;i<9*3;i++)
1312 ptr[i]=tmpdata[i];
1313 return tmp/3;
1314 /* return number of coords, not floats */
1315 }
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. */
1320 buf1=xfp->buf1;
1321 buf2=xfp->buf2;
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);
1326
1327 sizeint[0] = maxint[0] - minint[0]+1;
1328 sizeint[1] = maxint[1] - minint[1]+1;
1329 sizeint[2] = maxint[2] - minint[2]+1;
1330
1331 /* check if one of the sizes is to big to be multiplied */
1332 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
1333 {
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 */
1338 }
1339 else
1340 {
1341 bitsize = sizeofints(3, sizeint);
1342 }
1343
1344 if (xdrfile_read_int(&smallidx,1,xfp) == 0)
1345 return 0;
1346 tmp=smallidx+8;
1347 maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
1348 minidx = maxidx - 8; /* often this equal smallidx */
1349 tmp = smallidx-1;
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];
1355
1356 /* buf2[0] holds the length in bytes */
1357
1358 if (xdrfile_read_int(buf2,1,xfp) == 0)
1359 return 0;
1360 if (xdrfile_read_opaque((char *)&(buf2[3]),(unsigned int)buf2[0],xfp) == 0)
1361 return 0;
1362 buf2[0] = buf2[1] = buf2[2] = 0;
1363
1364 lfp = ptr;
1365 inv_precision = 1.0 / * precision;
1366 run = 0;
1367 i = 0;
1368 lip = buf1;
1369 while ( i < lsize )
1370 {
1371 thiscoord = (int *)(lip) + i * 3;
1372
1373 if (bitsize == 0)
1374 {
1375 thiscoord[0] = decodebits(buf2, bitsizeint[0]);
1376 thiscoord[1] = decodebits(buf2, bitsizeint[1]);
1377 thiscoord[2] = decodebits(buf2, bitsizeint[2]);
1378 } else {
1379 decodeints(buf2, 3, bitsize, sizeint, thiscoord);
1380 }
1381
1382 i++;
1383 thiscoord[0] += minint[0];
1384 thiscoord[1] += minint[1];
1385 thiscoord[2] += minint[2];
1386
1387 prevcoord[0] = thiscoord[0];
1388 prevcoord[1] = thiscoord[1];
1389 prevcoord[2] = thiscoord[2];
1390
1391 flag = decodebits(buf2, 1);
1392 is_smaller = 0;
1393 if (flag == 1)
1394 {
1395 run = decodebits(buf2, 5);
1396 is_smaller = run % 3;
1397 run -= is_smaller;
1398 is_smaller--;
1399 }
1400 if (run > 0)
1401 {
1402 thiscoord += 3;
1403 for (k = 0; k < run; k+=3)
1404 {
1405 decodeints(buf2, 3, smallidx, sizesmall, thiscoord);
1406 i++;
1407 thiscoord[0] += prevcoord[0] - smallnum;
1408 thiscoord[1] += prevcoord[1] - smallnum;
1409 thiscoord[2] += prevcoord[2] - smallnum;
1410 if (k == 0)
1411 {
1412 /* interchange first with second atom for better
1413 * compression of water molecules
1414 */
1415 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1416 prevcoord[0] = tmp;
1417 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1418 prevcoord[1] = tmp;
1419 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1420 prevcoord[2] = tmp;
1421 *lfp++ = prevcoord[0] * inv_precision;
1422 *lfp++ = prevcoord[1] * inv_precision;
1423 *lfp++ = prevcoord[2] * inv_precision;
1424 }
1425 else
1426 {
1427 prevcoord[0] = thiscoord[0];
1428 prevcoord[1] = thiscoord[1];
1429 prevcoord[2] = thiscoord[2];
1430 }
1431 *lfp++ = thiscoord[0] * inv_precision;
1432 *lfp++ = thiscoord[1] * inv_precision;
1433 *lfp++ = thiscoord[2] * inv_precision;
1434 }
1435 } else {
1436 *lfp++ = thiscoord[0] * inv_precision;
1437 *lfp++ = thiscoord[1] * inv_precision;
1438 *lfp++ = thiscoord[2] * inv_precision;
1439 }
1440 smallidx += is_smaller;
1441 if (is_smaller < 0) {
1442 smallnum = smaller;
1443 if (smallidx > FIRSTIDX) {
1444 smaller = magicints[smallidx - 1] /2;
1445 } else {
1446 smaller = 0;
1447 }
1448 } else if (is_smaller > 0) {
1449 smaller = smallnum;
1450 smallnum = magicints[smallidx] / 2;
1451 }
1452 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1453 }
1454 return *size;
1455 }
1456
1457 int
1458 xdrfile_compress_coord_double(double *ptr,
1459 int size,
1460 double precision,
1461 XDRFILE* xfp)
1462 {
1463 int minint[3], maxint[3], mindiff, *lip, diff;
1464 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
1465 int minidx, maxidx;
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;
1469 double *lfp;
1470 float float_prec, lf,tmpdata[30];
1471 int tmp, tmpsum, *thiscoord, prevcoord[3];
1472 unsigned int tmpcoord[30];
1473 int errval=1;
1474 unsigned int bitsize;
1475
1476 bitsizeint[0] = 0;
1477 bitsizeint[1] = 0;
1478 bitsizeint[2] = 0;
1479
1480 if(xfp==NULL)
1481 return -1;
1482 size3=3*size;
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");
1486 return -1;
1487 }
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");
1492 return -1;
1493 }
1494 }
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 */
1498 if(size<=9) {
1499 for(i=0;i<9*3;i++)
1500 tmpdata[i]=ptr[i];
1501 return xdrfile_write_float(tmpdata,size3,xfp)/3;
1502 /* return number of coords, not floats */
1503 }
1504 /* Compression-time if we got here. Write precision first */
1505 if (precision <= 0)
1506 precision = 1000;
1507 float_prec=precision;
1508 xdrfile_write_float(&float_prec,1,xfp);
1509 /* avoid repeated pointer dereferencing. */
1510 buf1=xfp->buf1;
1511 buf2=xfp->buf2;
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;
1516 prevrun = -1;
1517 lfp = ptr;
1518 lip = buf1;
1519 mindiff = INT_MAX;
1520 oldlint1 = oldlint2 = oldlint3 = 0;
1521 while(lfp < ptr + size3 ) {
1522 /* find nearest integer */
1523 if (*lfp >= 0.0)
1524 lf = (float)*lfp * float_prec + 0.5;
1525 else
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");
1530 errval=0;
1531 }
1532 lint1 = lf;
1533 if (lint1 < minint[0]) minint[0] = lint1;
1534 if (lint1 > maxint[0]) maxint[0] = lint1;
1535 *lip++ = lint1;
1536 lfp++;
1537 if (*lfp >= 0.0)
1538 lf = (float)*lfp * float_prec + 0.5;
1539 else
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");
1544 errval=0;
1545 }
1546 lint2 = lf;
1547 if (lint2 < minint[1]) minint[1] = lint2;
1548 if (lint2 > maxint[1]) maxint[1] = lint2;
1549 *lip++ = lint2;
1550 lfp++;
1551 if (*lfp >= 0.0)
1552 lf = (float)*lfp * float_prec + 0.5;
1553 else
1554 lf = (float)*lfp * float_prec - 0.5;
1555 if (fabs(lf) > INT_MAX-2) {
1556 errval=0;
1557 }
1558 lint3 = lf;
1559 if (lint3 < minint[2]) minint[2] = lint3;
1560 if (lint3 > maxint[2]) maxint[2] = lint3;
1561 *lip++ = lint3;
1562 lfp++;
1563 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
1564 if (diff < mindiff && lfp > ptr + 3)
1565 mindiff = diff;
1566 oldlint1 = lint1;
1567 oldlint2 = lint2;
1568 oldlint3 = lint3;
1569 }
1570 xdrfile_write_int(minint,3,xfp);
1571 xdrfile_write_int(maxint,3,xfp);
1572
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
1578 */
1579 fprintf(stderr,"Internal overflow compressing coordinates.\n");
1580 errval=0;
1581 }
1582 sizeint[0] = maxint[0] - minint[0]+1;
1583 sizeint[1] = maxint[1] - minint[1]+1;
1584 sizeint[2] = maxint[2] - minint[2]+1;
1585
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 */
1592 } else {
1593 bitsize = sizeofints(3, sizeint);
1594 }
1595 lip = buf1;
1596 luip = (unsigned int *) buf1;
1597 smallidx = FIRSTIDX;
1598 while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
1599 smallidx++;
1600 }
1601 xdrfile_write_int(&smallidx,1,xfp);
1602 tmp=smallidx+8;
1603 maxidx = (LASTIDX<tmp) ? LASTIDX : tmp;
1604 minidx = maxidx - 8; /* often this equal smallidx */
1605 tmp=smallidx-1;
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;
1612 i = 0;
1613 while (i < size) {
1614 is_small = 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) {
1620 is_smaller = 1;
1621 } else if (smallidx > minidx) {
1622 is_smaller = -1;
1623 } else {
1624 is_smaller = 0;
1625 }
1626 if (i + 1 < size) {
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
1632 */
1633 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
1634 thiscoord[3] = tmp;
1635 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
1636 thiscoord[4] = tmp;
1637 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
1638 thiscoord[5] = tmp;
1639 is_small = 1;
1640 }
1641 }
1642 tmpcoord[0] = thiscoord[0] - minint[0];
1643 tmpcoord[1] = thiscoord[1] - minint[1];
1644 tmpcoord[2] = thiscoord[2] - minint[2];
1645 if (bitsize == 0) {
1646 encodebits(buf2, bitsizeint[0], tmpcoord[0]);
1647 encodebits(buf2, bitsizeint[1], tmpcoord[1]);
1648 encodebits(buf2, bitsizeint[2], tmpcoord[2]);
1649 } else {
1650 encodeints(buf2, 3, bitsize, sizeint, tmpcoord);
1651 }
1652 prevcoord[0] = thiscoord[0];
1653 prevcoord[1] = thiscoord[1];
1654 prevcoord[2] = thiscoord[2];
1655 thiscoord = thiscoord + 3;
1656 i++;
1657
1658 run = 0;
1659 if (is_small == 0 && is_smaller == -1)
1660 is_smaller = 0;
1661 while (is_small && run < 8*3) {
1662 tmpsum=0;
1663 for(j=0;j<3;j++) {
1664 tmp=thiscoord[j] - prevcoord[j];
1665 tmpsum+=tmp*tmp;
1666 }
1667 if (is_smaller == -1 && tmpsum >= smaller * smaller) {
1668 is_smaller = 0;
1669 }
1670
1671 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
1672 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
1673 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
1674
1675 prevcoord[0] = thiscoord[0];
1676 prevcoord[1] = thiscoord[1];
1677 prevcoord[2] = thiscoord[2];
1678
1679 i++;
1680 thiscoord = thiscoord + 3;
1681 is_small = 0;
1682 if (i < size &&
1683 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
1684 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
1685 abs(thiscoord[2] - prevcoord[2]) < smallnum) {
1686 is_small = 1;
1687 }
1688 }
1689 if (run != prevrun || is_smaller != 0) {
1690 prevrun = run;
1691 encodebits(buf2, 1, 1); /* flag the change in run-length */
1692 encodebits(buf2, 5, run+is_smaller+1);
1693 } else {
1694 encodebits(buf2, 1, 0); /* flag the fact that runlength did not change */
1695 }
1696 for (k=0; k < run; k+=3) {
1697 encodeints(buf2, 3, smallidx, sizesmall, &tmpcoord[k]);
1698 }
1699 if (is_smaller != 0) {
1700 smallidx += is_smaller;
1701 if (is_smaller < 0) {
1702 smallnum = smaller;
1703 smaller = magicints[smallidx-1] / 2;
1704 } else {
1705 smaller = smallnum;
1706 smallnum = magicints[smallidx] / 2;
1707 }
1708 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1709 }
1710 }
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])
1715 return size;
1716 else
1717 return -1;
1718 }
1719
1720
1721 /* Dont try do document Fortran interface, since
1722 * Doxygen barfs at the F77_FUNC macro
1723 */
1724 #ifndef DOXYGEN
1725
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 */
1738
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);
1742
1743
1744 void
1745 F77_FUNC(xdropen,XDROPEN)(int *fid, char *filename, char *mode,
1746 int fn_len, int mode_len)
1747 {
1748 char cfilename[512];
1749 char cmode[5];
1750 int i;
1751
1752 /* zero array at first invocation */
1753 if(f77init) {
1754 for(i=0;i<MAX_FORTRAN_XDR;i++)
1755 f77xdr[i]=NULL;
1756 f77init=0;
1757 }
1758 i=0;
1759
1760 /* f77xdr is always smaller or equal to MAX_FORTRAN_XDR */
1761 while(i<MAX_FORTRAN_XDR && f77xdr[i]!=NULL)
1762 i++;
1763 if(i==MAX_FORTRAN_XDR) {
1764 *fid = -1;
1765 } else if (ftocstr(cfilename, sizeof(cfilename), filename, fn_len)) {
1766 *fid = -1;
1767 } else if (ftocstr(cmode, sizeof(cmode), mode,mode_len)) {
1768 *fid = -1;
1769 } else {
1770 f77xdr[i]=xdrfile_open(cfilename,cmode);
1771 /* return the index in the array as a fortran file handle */
1772 *fid=i;
1773 }
1774 }
1775
1776 void
1777 F77_FUNC(xdrclose,XDRCLOSE)(int *fid)
1778 {
1779 /* first close it */
1780 xdrfile_close(f77xdr[*fid]);
1781 /* the remove it from file handle list */
1782 f77xdr[*fid]=NULL;
1783 }
1784
1785
1786 void
1787 F77_FUNC(xdrrint,XDRRINT)(int *fid, int *data, int *ndata, int *ret)
1788 {
1789 *ret = xdrfile_read_int(data,*ndata,f77xdr[*fid]);
1790 }
1791
1792 void
1793 F77_FUNC(xdrwint,XDRWINT)(int *fid, int *data, int *ndata, int *ret)
1794 {
1795 *ret = xdrfile_write_int(data,*ndata,f77xdr[*fid]);
1796 }
1797
1798 void
1799 F77_FUNC(xdrruint,XDRRUINT)(int *fid, unsigned int *data, int *ndata, int *ret)
1800 {
1801 *ret = xdrfile_read_uint(data,*ndata,f77xdr[*fid]);
1802 }
1803
1804 void
1805 F77_FUNC(xdrwuint,XDRWUINT)(int *fid, unsigned int *data, int *ndata, int *ret)
1806 {
1807 *ret = xdrfile_write_uint(data,*ndata,f77xdr[*fid]);
1808 }
1809
1810 void
1811 F77_FUNC(xdrrchar,XDRRCHAR)(int *fid, char *ip, int *ndata, int *ret)
1812 {
1813 *ret = xdrfile_read_char(ip,*ndata,f77xdr[*fid]);
1814 }
1815
1816 void
1817 F77_FUNC(xdrwchar,XDRWCHAR)(int *fid, char *ip, int *ndata, int *ret)
1818 {
1819 *ret = xdrfile_write_char(ip,*ndata,f77xdr[*fid]);
1820 }
1821
1822 void
1823 F77_FUNC(xdrruchar,XDRRUCHAR)(int *fid, unsigned char *ip, int *ndata, int *ret)
1824 {
1825 *ret = xdrfile_read_uchar(ip,*ndata,f77xdr[*fid]);
1826 }
1827
1828 void
1829 F77_FUNC(xdrwuchar,XDRWUCHAR)(int *fid, unsigned char *ip, int *ndata, int *ret)
1830 {
1831 *ret = xdrfile_write_uchar(ip,*ndata,f77xdr[*fid]);
1832 }
1833
1834 void
1835 F77_FUNC(xdrrshort,XDRRSHORT)(int *fid, short *ip, int *ndata, int *ret)
1836 {
1837 *ret = xdrfile_read_short(ip,*ndata,f77xdr[*fid]);
1838 }
1839
1840 void
1841 F77_FUNC(xdrwshort,XDRWSHORT)(int *fid, short *ip, int *ndata, int *ret)
1842 {
1843 *ret = xdrfile_write_short(ip,*ndata,f77xdr[*fid]);
1844 }
1845
1846 void
1847 F77_FUNC(xdrrushort,XDRRUSHORT)(int *fid, unsigned short *ip, int *ndata, int *ret)
1848 {
1849 *ret = xdrfile_read_ushort(ip,*ndata,f77xdr[*fid]);
1850 }
1851
1852 void
1853 F77_FUNC(xdrwushort,XDRWUSHORT)(int *fid, unsigned short *ip, int *ndata, int *ret)
1854 {
1855 *ret = xdrfile_write_ushort(ip,*ndata,f77xdr[*fid]);
1856 }
1857
1858 void
1859 F77_FUNC(xdrrsingle,XDRRSINGLE)(int *fid, float *data, int *ndata, int *ret)
1860 {
1861 *ret = xdrfile_read_float(data,*ndata,f77xdr[*fid]);
1862 }
1863
1864 void
1865 F77_FUNC(xdrwsingle,XDRWSINGLE)(int *fid, float *data, int *ndata, int *ret)
1866 {
1867 *ret = xdrfile_write_float(data,*ndata,f77xdr[*fid]);
1868 }
1869
1870 void
1871 F77_FUNC(xdrrdouble,XDRRDOUBLE)(int *fid, double *data, int *ndata, int *ret)
1872 {
1873 *ret = xdrfile_read_double(data,*ndata,f77xdr[*fid]);
1874 }
1875
1876 void
1877 F77_FUNC(xdrwdouble,XDRWDOUBLE)(int *fid, double *data, int *ndata, int *ret)
1878 {
1879 *ret = xdrfile_write_double(data,*ndata,f77xdr[*fid]);
1880 }
1881
1882 static int ftocstr(char *dest, int destlen, char *src, int srclen)
1883 {
1884 char *p;
1885
1886 p = src + srclen;
1887 while ( --p >= src && *p == ' ' );
1888 srclen = p - src + 1;
1889 destlen--;
1890 dest[0] = 0;
1891 if (srclen > destlen)
1892 return 1;
1893 while (srclen--)
1894 (*dest++ = *src++);
1895 *dest = '\0';
1896 return 0;
1897 }
1898
1899
1900 static int ctofstr(char *dest, int destlen, char *src)
1901 {
1902 while (destlen && *src) {
1903 *dest++ = *src++;
1904 destlen--;
1905 }
1906 while (destlen--)
1907 *dest++ = ' ';
1908 return 0;
1909 }
1910
1911
1912 void
1913 F77_FUNC(xdrrstring,XDRRSTRING)(int *fid, char *str, int *ret, int len)
1914 {
1915 char *cstr;
1916
1917 if((cstr=(char*)malloc((len+1)*sizeof(char)))==NULL) {
1918 *ret = 0;
1919 return;
1920 }
1921 if (ftocstr(cstr, len+1, str, len)) {
1922 *ret = 0;
1923 free(cstr);
1924 return;
1925 }
1926
1927 *ret = xdrfile_read_string(cstr, len+1,f77xdr[*fid]);
1928 ctofstr( str, len , cstr);
1929 free(cstr);
1930 }
1931
1932 void
1933 F77_FUNC(xdrwstring,XDRWSTRING)(int *fid, char *str, int *ret, int len)
1934 {
1935 char *cstr;
1936
1937 if((cstr=(char*)malloc((len+1)*sizeof(char)))==NULL) {
1938 *ret = 0;
1939 return;
1940 }
1941 if (ftocstr(cstr, len+1, str, len)) {
1942 *ret = 0;
1943 free(cstr);
1944 return;
1945 }
1946
1947 *ret = xdrfile_write_string(cstr, f77xdr[*fid]);
1948 ctofstr( str, len , cstr);
1949 free(cstr);
1950 }
1951
1952 void
1953 F77_FUNC(xdrropaque,XDRROPAQUE)(int *fid, char *data, int *ndata, int *ret)
1954 {
1955 *ret = xdrfile_read_opaque(data,*ndata,f77xdr[*fid]);
1956 }
1957
1958 void
1959 F77_FUNC(xdrwopaque,XDRWOPAQUE)(int *fid, char *data, int *ndata, int *ret)
1960 {
1961 *ret = xdrfile_write_opaque(data,*ndata,f77xdr[*fid]);
1962 }
1963
1964
1965 /* Write single-precision compressed 3d coordinates */
1966 void
1967 F77_FUNC(xdrccs,XDRCCS)(int *fid, float *data, int *ncoord,
1968 float *precision, int *ret)
1969 {
1970 *ret = xdrfile_compress_coord_float(data,*ncoord,*precision,f77xdr[*fid]);
1971 }
1972
1973
1974 /* Read single-precision compressed 3d coordinates */
1975 void
1976 F77_FUNC(xdrdcs,XDRDCS)(int *fid, float *data, int *ncoord,
1977 float *precision, int *ret)
1978 {
1979 *ret = xdrfile_decompress_coord_float(data,ncoord,precision,f77xdr[*fid]);
1980 }
1981
1982
1983 /* Write compressed 3d coordinates from double precision data */
1984 void
1985 F77_FUNC(xdrccd,XDRCCD)(int *fid, double *data, int *ncoord,
1986 double *precision, int *ret)
1987 {
1988 *ret = xdrfile_compress_coord_double(data,*ncoord,*precision,f77xdr[*fid]);
1989 }
1990
1991 /* Read compressed 3d coordinates into double precision data */
1992 void
1993 F77_FUNC(xddcd,XDRDCD)(int *fid, double *data, int *ncoord,
1994 double *precision, int *ret)
1995 {
1996 *ret = xdrfile_decompress_coord_double(data,ncoord,precision,f77xdr[*fid]);
1997 }
1998
1999
2000
2001
2002
2003
2004
2005 #endif /* DOXYGEN */
2006
2007 /*************************************************************
2008 * End of higher-level routines - dont change things below! *
2009 *************************************************************/
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
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
2037
2038 /*
2039 * What follows is a modified version of the Sun XDR code. For reference
2040 * we include their copyright and license:
2041 *
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.
2048 *
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.
2052 *
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.
2056 *
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.
2060 *
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.
2064 *
2065 * Sun Microsystems, Inc.
2066 * 2550 Garcia Avenue
2067 * Mountain View, California 94043
2068 */
2069
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 */
2076
2077
2078
2079
2080
2081 typedef int (*xdrproc_t) (XDR *, void *,...);
2082
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)
2091
2092 #define BYTES_PER_XDR_UNIT 4
2093 static char xdr_zero[BYTES_PER_XDR_UNIT] = {0, 0, 0, 0};
2094
2095 static int32_t
2096 xdr_swapbytes(int32_t x)
2097 {
2098 int32_t y,i;
2099 char *px=(char *)&x;
2100 char *py=(char *)&y;
2101
2102 for(i=0;i<4;i++)
2103 py[i]=px[3-i];
2104
2105 return y;
2106 }
2107
2108 static int32_t
2109 xdr_htonl(int32_t x)
2110 {
2111 int s=0x1234;
2112 if( *((char *)&s)==(char)0x34) {
2113 /* smallendian,swap bytes */
2114 return xdr_swapbytes(x);
2115 } else {
2116 /* bigendian, do nothing */
2117 return x;
2118 }
2119 }
2120
2121 static int32_t
2122 xdr_ntohl(int x)
2123 {
2124 int s=0x1234;
2125 if( *((char *)&s)==(char)0x34) {
2126 /* smallendian, swap bytes */
2127 return xdr_swapbytes(x);
2128 } else {
2129 /* bigendian, do nothing */
2130 return x;
2131 }
2132 }
2133
2134 static int
2135 xdr_int (XDR *xdrs, int *ip)
2136 {
2137 int32_t i32;
2138
2139 switch (xdrs->x_op)
2140 {
2141 case XDR_ENCODE:
2142 i32 = (int32_t) *ip;
2143 return xdr_putlong (xdrs, &i32);
2144
2145 case XDR_DECODE:
2146 if (!xdr_getlong (xdrs, &i32))
2147 {
2148 return 0;
2149 }
2150 *ip = (int) i32;
2151 case XDR_FREE:
2152 return 1;
2153 }
2154 return 0;
2155 }
2156
2157 static int
2158 xdr_u_int (XDR *xdrs, unsigned int *up)
2159 {
2160 uint32_t ui32;
2161
2162 switch (xdrs->x_op)
2163 {
2164 case XDR_ENCODE:
2165 ui32 = (uint32_t) * up;
2166 return xdr_putlong (xdrs, (int32_t *)&ui32);
2167
2168 case XDR_DECODE:
2169 if (!xdr_getlong (xdrs, (int32_t *)&ui32))
2170 {
2171 return 0;
2172 }
2173 *up = (uint32_t) ui32;
2174 case XDR_FREE:
2175 return 1;
2176 }
2177 return 0;
2178 }
2179
2180 static int
2181 xdr_short (XDR *xdrs, short *sp)
2182 {
2183 int32_t i32;
2184
2185 switch (xdrs->x_op)
2186 {
2187 case XDR_ENCODE:
2188 i32 = (int32_t) *sp;
2189 return xdr_putlong (xdrs, &i32);
2190
2191 case XDR_DECODE:
2192 if (!xdr_getlong (xdrs, &i32))
2193 {
2194 return 0;
2195 }
2196 *sp = (short) i32;
2197 return 1;
2198
2199 case XDR_FREE:
2200 return 1;
2201 }
2202 return 0;
2203 }
2204
2205 static int
2206 xdr_u_short (XDR *xdrs, unsigned short *sp)
2207 {
2208 uint32_t ui32;
2209
2210 switch (xdrs->x_op)
2211 {
2212 case XDR_ENCODE:
2213 ui32 = (uint32_t) *sp;
2214 return xdr_putlong (xdrs, (int32_t *)&ui32);
2215
2216 case XDR_DECODE:
2217 if (!xdr_getlong (xdrs, (int32_t *)&ui32))
2218 {
2219 return 0;
2220 }
2221 *sp = (unsigned short) ui32;
2222 return 1;
2223
2224 case XDR_FREE:
2225 return 1;
2226 }
2227 return 0;
2228 }
2229
2230 static int
2231 xdr_char (XDR *xdrs, char *cp)
2232 {
2233 int i;
2234
2235 i = (*cp);
2236 if (!xdr_int (xdrs, &i))
2237 {
2238 return 0;
2239 }
2240 *cp = i;
2241 return 1;
2242 }
2243
2244 static int
2245 xdr_u_char (XDR *xdrs, unsigned char *cp)
2246 {
2247 unsigned int u;
2248
2249 u = (*cp);
2250 if (!xdr_u_int (xdrs, &u))
2251 {
2252 return 0;
2253 }
2254 *cp = u;
2255 return 1;
2256 }
2257
2258 /*
2259 * XDR opaque data
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.
2262 */
2263 static int
2264 xdr_opaque (XDR *xdrs, char *cp, unsigned int cnt)
2265 {
2266 unsigned int rndup;
2267 static char crud[BYTES_PER_XDR_UNIT];
2268
2269 /*
2270 * if no data we are done
2271 */
2272 if (cnt == 0)
2273 return 1;
2274
2275 /*
2276 * round byte count to full xdr units
2277 */
2278 rndup = cnt % BYTES_PER_XDR_UNIT;
2279 if (rndup > 0)
2280 rndup = BYTES_PER_XDR_UNIT - rndup;
2281
2282 switch (xdrs->x_op)
2283 {
2284 case XDR_DECODE:
2285 if (!xdr_getbytes (xdrs, cp, cnt))
2286 {
2287 return 0;
2288 }
2289 if (rndup == 0)
2290 return 1;
2291 return xdr_getbytes (xdrs, (char *)crud, rndup);
2292
2293 case XDR_ENCODE:
2294 if (!xdr_putbytes (xdrs, cp, cnt))
2295 {
2296 return 0;
2297 }
2298 if (rndup == 0)
2299 return 1;
2300 return xdr_putbytes (xdrs, xdr_zero, rndup);
2301
2302 case XDR_FREE:
2303 return 1;
2304 }
2305 #undef BYTES_PER_XDR_UNIT
2306 return 0;
2307 }
2308
2309
2310 /*
2311 * XDR null terminated ASCII strings
2312 */
2313 static int
2314 xdr_string (XDR *xdrs, char **cpp, unsigned int maxsize)
2315 {
2316 char *sp = *cpp; /* sp is the actual string pointer */
2317 unsigned int size;
2318 unsigned int nodesize;
2319
2320 /*
2321 * first deal with the length since xdr strings are counted-strings
2322 */
2323 switch (xdrs->x_op)
2324 {
2325 case XDR_FREE:
2326 if (sp == NULL)
2327 {
2328 return 1; /* already free */
2329 }
2330 /* fall through... */
2331 case XDR_ENCODE:
2332 if (sp == NULL)
2333 return 0;
2334 size = strlen (sp);
2335 break;
2336 case XDR_DECODE:
2337 break;
2338 }
2339 if (!xdr_u_int (xdrs, &size))
2340 {
2341 return 0;
2342 }
2343 if (size > maxsize)
2344 {
2345 return 0;
2346 }
2347 nodesize = size + 1;
2348
2349 /*
2350 * now deal with the actual bytes
2351 */
2352 switch (xdrs->x_op)
2353 {
2354 case XDR_DECODE:
2355 if (nodesize == 0)
2356 {
2357 return 1;
2358 }
2359 if (sp == NULL)
2360 *cpp = sp = (char *) malloc (nodesize);
2361 if (sp == NULL)
2362 {
2363 (void) fputs ("xdr_string: out of memory\n", stderr);
2364 return 0;
2365 }
2366 sp[size] = 0;
2367 /* fall into ... */
2368
2369 case XDR_ENCODE:
2370 return xdr_opaque (xdrs, sp, size);
2371
2372 case XDR_FREE:
2373 free (sp);
2374 *cpp = NULL;
2375 return 1;
2376 }
2377 return 0;
2378 }
2379
2380
2381
2382 /* Floating-point stuff */
2383
2384 static int
2385 xdr_float(XDR *xdrs, float *fp)
2386 {
2387 switch (xdrs->x_op) {
2388
2389 case XDR_ENCODE:
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));
2395 }
2396 break;
2397
2398 case XDR_DECODE:
2399 if (sizeof(float) == sizeof(int32_t))
2400 return (xdr_getlong(xdrs, (int32_t *)fp));
2401 else if (sizeof(float) == sizeof(int)) {
2402 int32_t tmp;
2403 if (xdr_getlong(xdrs, &tmp)) {
2404 *(int *)fp = tmp;
2405 return (1);
2406 }
2407 }
2408 break;
2409
2410 case XDR_FREE:
2411 return (1);
2412 }
2413 return (0);
2414 }
2415
2416
2417 static int
2418 xdr_double(XDR *xdrs, double *dp)
2419 {
2420 /* Gromacs detects floating-point stuff at compile time, which is faster */
2421 #ifdef GROMACS
2422 # ifndef FLOAT_FORMAT_IEEE754
2423 # error non-IEEE floating point system, or you defined GROMACS yourself...
2424 # endif
2425 int LSW;
2426 # ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
2427 int LSW=1;
2428 # else
2429 int LSW=0;
2430 # endif /* Big endian word order */
2431 #else
2432 /* Outside Gromacs we rely on dynamic detection of FP order. */
2433 int LSW; /* Least significant fp word */
2434
2435 double x=0.987654321; /* Just a number */
2436 unsigned char ix = *((char *)&x);
2437
2438 /* Possible representations in IEEE double precision:
2439 * (S=small endian, B=big endian)
2440 *
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
2446 */
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");
2455 abort();
2456 }
2457 #endif /* end of dynamic detection of fp word order */
2458
2459 switch (xdrs->x_op) {
2460
2461 case XDR_ENCODE:
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;
2468 int32_t tmp[2];
2469 tmp[0] = ip[!LSW];
2470 tmp[1] = ip[LSW];
2471 return (xdr_putlong(xdrs, tmp) &&
2472 xdr_putlong(xdrs, tmp+1));
2473 }
2474 break;
2475
2476 case XDR_DECODE:
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;
2483 int32_t tmp[2];
2484 if (xdr_getlong(xdrs, tmp+!LSW) &&
2485 xdr_getlong(xdrs, tmp+LSW)) {
2486 ip[0] = tmp[0];
2487 ip[1] = tmp[1];
2488 return (1);
2489 }
2490 }
2491 break;
2492
2493 case XDR_FREE:
2494 return (1);
2495 }
2496 return (0);
2497 }
2498
2499
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 *);
2507
2508 /*
2509 * Ops vector for stdio type XDR
2510 */
2511 static const struct xdr_ops xdrstdio_ops =
2512 {
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 */
2520 };
2521
2522 /*
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.
2526 */
2527 static void
2528 xdrstdio_create (XDR *xdrs, FILE *file, enum xdr_op op)
2529 {
2530 xdrs->x_op = op;
2531
2532 xdrs->x_ops = (struct xdr_ops *) &xdrstdio_ops;
2533 xdrs->x_private = (char *) file;
2534 }
2535
2536 /*
2537 * Destroy a stdio xdr stream.
2538 * Cleans up the xdr stream handle xdrs previously set up by xdrstdio_create.
2539 */
2540 static void
2541 xdrstdio_destroy (XDR *xdrs)
2542 {
2543 (void) fflush ((FILE *) xdrs->x_private);
2544 /* xx should we close the file ?? */
2545 }
2546
2547 static int
2548 xdrstdio_getlong (XDR *xdrs, int32_t *lp)
2549 {
2550 int32_t mycopy;
2551
2552 if (fread ((char *) & mycopy, 4, 1, (FILE *) xdrs->x_private) != 1)
2553 return 0;
2554 *lp = (int32_t) xdr_ntohl (mycopy);
2555 return 1;
2556 }
2557
2558 static int
2559 xdrstdio_putlong (XDR *xdrs, int32_t *lp)
2560 {
2561 int32_t mycopy = xdr_htonl (*lp);
2562 lp = &mycopy;
2563 if (fwrite ((char *) lp, 4, 1, (FILE *) xdrs->x_private) != 1)
2564 return 0;
2565 return 1;
2566 }
2567
2568 static int
2569 xdrstdio_getbytes (XDR *xdrs, char *addr, unsigned int len)
2570 {
2571 if ((len != 0) && (fread (addr, (int) len, 1,
2572 (FILE *) xdrs->x_private) != 1))
2573 return 0;
2574 return 1;
2575 }
2576
2577 static int
2578 xdrstdio_putbytes (XDR *xdrs, char *addr, unsigned int len)
2579 {
2580 if ((len != 0) && (fwrite (addr, (int) len, 1,
2581 (FILE *) xdrs->x_private) != 1))
2582 return 0;
2583 return 1;
2584 }
2585
2586 /* 32 bit fileseek operations */
2587 static unsigned int
2588 xdrstdio_getpos (XDR *xdrs)
2589 {
2590 return (unsigned int) ftell ((FILE *) xdrs->x_private);
2591 }
2592
2593 static int
2594 xdrstdio_setpos (XDR *xdrs, unsigned int pos)
2595 {
2596 return fseek ((FILE *) xdrs->x_private, pos, 0) < 0 ? 0 : 1;
2597 }
2598
2599
2600
2601 #endif /* HAVE_RPC_XDR_H not defined */
XTC Library