| blob | 111839832ef986c0c23675be89c39f7323d9202b |
1 /*____________________________________________________________________________
2 |
3 | libxdrf - portable fortran interface to xdr. some xdr routines
4 | are C routines for compressed coordinates
5 |
6 | version 1.1
7 |
8 | This collection of routines is intended to write and read
9 | data in a portable way to a file, so data written on one type
10 | of machine can be read back on a different type.
11 |
12 | all fortran routines use an integer 'xdrid', which is an id to the
13 | current xdr file, and is set by xdrfopen.
14 | most routines have in integer 'ret' which is the return value.
15 | The value of 'ret' is zero on failure, and most of the time one
16 | on succes.
17 |
18 | There are three routines useful for C users:
19 | xdropen(), xdrclose(), xdr3dfcoord().
20 | The first two replace xdrstdio_create and xdr_destroy, and *must* be
21 | used when you plan to use xdr3dfcoord(). (they are also a bit
22 | easier to interface). For writing data other than compressed coordinates
23 | you should use the standard C xdr routines (see xdr man page)
24 |
25 | xdrfopen(xdrid, filename, mode, ret)
26 | character *(*) filename
27 | character *(*) mode
28 |
29 | this will open the file with the given filename (string)
30 | and the given mode, it returns an id in xdrid, which is
31 | to be used in all other calls to xdrf routines.
32 | mode is 'w' to create, or update an file, for all other
33 | values of mode the file is opened for reading
34 |
35 | you need to call xdrfclose to flush the output and close
36 | the file.
37 | Note that you should not use xdrstdio_create, which comes with the
38 | standard xdr library
39 |
40 | xdrfclose(xdrid, ret)
41 | flush the data to the file, and closes the file;
42 | You should not use xdr_destroy (which comes standard with
43 | the xdr libraries.
44 |
45 | xdrfbool(xdrid, bp, ret)
46 | integer pb
47 |
48 | This filter produces values of either 1 or 0
49 |
50 | xdrfchar(xdrid, cp, ret)
51 | character cp
52 |
53 | filter that translate between characters and their xdr representation
54 | Note that the characters in not compressed and occupies 4 bytes.
55 |
56 | xdrfdouble(xdrid, dp, ret)
57 | double dp
58 |
59 | read/write a double.
60 |
61 | xdrffloat(xdrid, fp, ret)
62 | float fp
63 |
64 | read/write a float.
65 |
66 | xdrfint(xdrid, ip, ret)
67 | integer ip
68 |
69 | read/write integer.
70 |
71 | xdrflong(xdrid, lp, ret)
72 | integer lp
73 |
74 | this routine has a possible portablility problem due to 64 bits longs.
75 |
76 | xdrfshort(xdrid, sp, ret)
77 | integer *2 sp
78 |
79 | xdrfstring(xdrid, sp, maxsize, ret)
80 | character *(*)
81 | integer maxsize
82 |
83 | read/write a string, with maximum length given by maxsize
84 |
85 | xdrfwrapstring(xdris, sp, ret)
86 | character *(*)
87 |
88 | read/write a string (it is the same as xdrfstring accept that it finds
89 | the stringlength itself.
90 |
91 | xdrfvector(xdrid, cp, size, xdrfproc, ret)
92 | character *(*)
93 | integer size
94 | external xdrfproc
95 |
96 | read/write an array pointed to by cp, with number of elements
97 | defined by 'size'. the routine 'xdrfproc' is the name
98 | of one of the above routines to read/write data (like xdrfdouble)
99 | In contrast with the c-version you don't need to specify the
100 | byte size of an element.
101 | xdrfstring is not allowed here (it is in the c version)
102 |
103 | xdrf3dfcoord(xdrid, fp, size, precision, ret)
104 | real (*) fp
105 | real precision
106 | integer size
107 |
108 | this is *NOT* a standard xdr routine. I named it this way, because
109 | it invites people to use the other xdr routines.
110 | It is introduced to store specifically 3d coordinates of molecules
111 | (as found in molecular dynamics) and it writes it in a compressed way.
112 | It starts by multiplying all numbers by precision and
113 | rounding the result to integer. effectively converting
114 | all floating point numbers to fixed point.
115 | it uses an algorithm for compression that is optimized for
116 | molecular data, but could be used for other 3d coordinates
117 | as well. There is subtantial overhead involved, so call this
118 | routine only if you have a large number of coordinates to read/write
119 |
120 | ________________________________________________________________________
121 |
122 | Below are the rotuiens to be used by C programmers. Use the 'normal'
123 | xdr routines to write integers, floats, etc (see man xdr)
124 |
125 | int xdropen(XDR *xdrs, const char *filename, const char *type)
126 | This will open the file with the given filename and the
127 | given mode. You should pass it an allocated XDR struct
128 | in xdrs, to be used in all other calls to xdr routines.
129 | Mode is 'w' to create, or update an file, and for all
130 | other values of mode the file is opened for reading.
131 | You need to call xdrclose to flush the output and close
132 | the file.
133 |
134 | Note that you should not use xdrstdio_create, which
135 | comes with the standard xdr library.
136 |
137 | int xdrclose(XDR *xdrs)
138 | Flush the data to the file, and close the file;
139 | You should not use xdr_destroy (which comes standard
140 | with the xdr libraries).
141 |
142 | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
143 | This is \fInot\fR a standard xdr routine. I named it this
144 | way, because it invites people to use the other xdr
145 | routines.
146 |
147 | frans van hoesel hoesel@chem.rug.nl
148 */
150 #include <limits.h>
151 #include <malloc.h>
152 #include <math.h>
153 #include <stdio.h>
154 #include <stdlib.h>
155 #include "xdrf.h"
157 int ftocstr(char *, int, char *, int);
158 int ctofstr(char *, int, char *);
160 #define MAXID 20
161 static FILE *xdrfiles[MAXID];
162 static XDR *xdridptr[MAXID];
163 static char xdrmodes[MAXID];
164 static unsigned int cnt;
166 typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *);
168 void
169 FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret')
170 int *xdrid, *ret;
171 int *pb;
172 {
173 *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb);
174 cnt += sizeof(int);
175 }
177 void
178 FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret')
179 int *xdrid, *ret;
180 char *cp;
181 {
182 *ret = xdr_char(xdridptr[*xdrid], cp);
183 cnt += sizeof(char);
184 }
186 void
187 FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret')
188 int *xdrid, *ret;
189 double *dp;
190 {
191 *ret = xdr_double(xdridptr[*xdrid], dp);
192 cnt += sizeof(double);
193 }
195 void
196 FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret')
197 int *xdrid, *ret;
198 float *fp;
199 {
200 *ret = xdr_float(xdridptr[*xdrid], fp);
201 cnt += sizeof(float);
202 }
204 void
205 FUNCTION(xdrfint) ARGS(`xdrid, ip, ret')
206 int *xdrid, *ret;
207 int *ip;
208 {
209 *ret = xdr_int(xdridptr[*xdrid], ip);
210 cnt += sizeof(int);
211 }
213 void
214 FUNCTION(xdrflong) ARGS(`xdrid, lp, ret')
215 int *xdrid, *ret;
216 long *lp;
217 {
218 *ret = xdr_long(xdridptr[*xdrid], lp);
219 cnt += sizeof(long);
220 }
222 void
223 FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret')
224 int *xdrid, *ret;
225 short *sp;
226 {
227 *ret = xdr_short(xdridptr[*xdrid], sp);
228 cnt += sizeof(sp);
229 }
231 void
232 FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret')
233 int *xdrid, *ret;
234 unsigned char *ucp;
235 {
236 *ret = xdr_u_char(xdridptr[*xdrid], (u_char *)ucp);
237 cnt += sizeof(char);
238 }
240 void
241 FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret')
242 int *xdrid, *ret;
243 unsigned long *ulp;
244 {
245 *ret = xdr_u_long(xdridptr[*xdrid], (u_long *)ulp);
246 cnt += sizeof(unsigned long);
247 }
249 void
250 FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret')
251 int *xdrid, *ret;
252 unsigned short *usp;
253 {
254 *ret = xdr_u_short(xdridptr[*xdrid], (u_short *)usp);
255 cnt += sizeof(unsigned short);
256 }
258 void
259 FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret')
260 int *xdrid, *ret;
261 float *fp;
262 int *size;
263 float *precision;
264 {
265 *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
266 }
268 void
269 FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret')
270 int *xdrid, *ret;
271 STRING_ARG_DECL(sp);
272 int *maxsize;
273 {
274 char *tsp;
276 tsp = (char*) malloc((size_t)(((STRING_LEN(sp)) + 1) * sizeof(char)));
277 if (tsp == NULL) {
278 *ret = -1;
279 return;
280 }
281 if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) {
282 *ret = -1;
283 free(tsp);
284 return;
285 }
286 *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize);
287 ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
288 cnt += *maxsize;
289 free(tsp);
290 }
292 void
293 FUNCTION(xdrfwrapstring) ARGS(`xdrid, STRING_ARG(sp), ret')
294 int *xdrid, *ret;
295 STRING_ARG_DECL(sp);
296 {
297 char *tsp;
298 int maxsize;
299 maxsize = (STRING_LEN(sp)) + 1;
300 tsp = (char*) malloc((size_t)(maxsize * sizeof(char)));
301 if (tsp == NULL) {
302 *ret = -1;
303 return;
304 }
305 if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) {
306 *ret = -1;
307 free(tsp);
308 return;
309 }
310 *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
311 ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
312 cnt += maxsize;
313 free(tsp);
314 }
316 void
317 FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret')
318 int *xdrid, *ret;
319 caddr_t *cp;
320 int *ccnt;
321 {
322 *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
323 cnt += *ccnt;
324 }
326 void
327 FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret')
328 int *xdrid, *ret;
329 int *pos;
330 {
331 *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
332 }
334 void
335 FUNCTION(xdrf) ARGS(`xdrid, pos')
336 int *xdrid, *pos;
337 {
338 *pos = xdr_getpos(xdridptr[*xdrid]);
339 }
341 void
342 FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret')
343 int *xdrid, *ret;
344 char *cp;
345 int *size;
346 FUNCTION(xdrfproc) elproc;
347 {
348 int lcnt;
349 cnt = 0;
350 for (lcnt = 0; lcnt < *size; lcnt++) {
351 elproc(xdrid, (cp+cnt) , ret);
352 }
353 }
356 void
357 FUNCTION(xdrfclose) ARGS(`xdrid, ret')
358 int *xdrid;
359 int *ret;
360 {
361 *ret = xdrclose(xdridptr[*xdrid]);
362 cnt = 0;
363 }
365 void
366 FUNCTION(xdrfopen) ARGS(`xdrid, STRING_ARG(fp), STRING_ARG(mode), ret')
367 int *xdrid;
368 STRING_ARG_DECL(fp);
369 STRING_ARG_DECL(mode);
370 int *ret;
371 {
372 char fname[512];
373 char fmode[3];
375 if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) {
376 *ret = 0;
377 }
378 if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode),
379 STRING_LEN(mode))) {
380 *ret = 0;
381 }
383 *xdrid = xdropen(NULL, fname, fmode);
384 if (*xdrid == 0)
385 *ret = 0;
386 else
387 *ret = 1;
388 }
390 /*___________________________________________________________________________
391 |
392 | what follows are the C routines for opening, closing xdr streams
393 | and the routine to read/write compressed coordinates together
394 | with some routines to assist in this task (those are marked
395 | static and cannot be called from user programs)
396 */
397 #define MAXABS INT_MAX-2
399 #ifndef MIN
400 #define MIN(x,y) ((x) < (y) ? (x):(y))
401 #endif
402 #ifndef MAX
403 #define MAX(x,y) ((x) > (y) ? (x):(y))
404 #endif
405 #ifndef SQR
406 #define SQR(x) ((x)*(x))
407 #endif
408 static int magicints[] = {
409 0, 0, 0, 0, 0, 0, 0, 0, 0,
410 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
411 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
412 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
413 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
414 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
415 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
416 8388607, 10568983, 13316085, 16777216 };
418 #define FIRSTIDX 9
419 /* note that magicints[FIRSTIDX-1] == 0 */
420 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
423 /*__________________________________________________________________________
424 |
425 | xdropen - open xdr file
426 |
427 | This versions differs from xdrstdio_create, because I need to know
428 | the state of the file (read or write) so I can use xdr3dfcoord
429 | in eigther read or write mode, and the file descriptor
430 | so I can close the file (something xdr_destroy doesn't do).
431 |
432 */
434 int xdropen(XDR *xdrs, const char *filename, const char *type) {
435 static int init_done = 0;
436 enum xdr_op lmode;
437 int xdrid;
439 if (init_done == 0) {
440 for (xdrid = 1; xdrid < MAXID; xdrid++) {
441 xdridptr[xdrid] = NULL;
442 }
443 init_done = 1;
444 }
445 xdrid = 1;
446 while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
447 xdrid++;
448 }
449 if (xdrid == MAXID) {
450 return 0;
451 }
452 if (*type == 'w' || *type == 'W') {
453 type = "w+";
454 lmode = XDR_ENCODE;
455 } else if (*type == 'a' || *type == 'A') {
456 type = "a+";
457 lmode = XDR_ENCODE;
458 } else {
459 type = "r";
460 lmode = XDR_DECODE;
461 }
462 xdrfiles[xdrid] = fopen(filename, type);
463 if (xdrfiles[xdrid] == NULL) {
464 xdrs = NULL;
465 return 0;
466 }
467 xdrmodes[xdrid] = *type;
468 /* next test isn't usefull in the case of C language
469 * but is used for the Fortran interface
470 * (C users are expected to pass the address of an already allocated
471 * XDR staructure)
472 */
473 if (xdrs == NULL) {
474 xdridptr[xdrid] = (XDR *) malloc((size_t)sizeof(XDR));
475 xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
476 } else {
477 xdridptr[xdrid] = xdrs;
478 xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
479 }
480 return xdrid;
481 }
483 /*_________________________________________________________________________
484 |
485 | xdrclose - close a xdr file
486 |
487 | This will flush the xdr buffers, and destroy the xdr stream.
488 | It also closes the associated file descriptor (this is *not*
489 | done by xdr_destroy).
490 |
491 */
493 int xdrclose(XDR *xdrs) {
494 int xdrid;
496 if (xdrs == NULL) {
497 fprintf(stderr, "xdrclose: passed a NULL pointer\n");
498 exit(1);
499 }
500 for (xdrid = 1; xdrid < MAXID; xdrid++) {
501 if (xdridptr[xdrid] == xdrs) {
503 xdr_destroy(xdrs);
504 fclose(xdrfiles[xdrid]);
505 xdridptr[xdrid] = NULL;
506 return 1;
507 }
508 }
509 fprintf(stderr, "xdrclose: no such open xdr file\n");
510 exit(1);
512 /* to make some compilers happy: */
513 return 0;
514 }
516 /*____________________________________________________________________________
517 |
518 | sendbits - encode num into buf using the specified number of bits
519 |
520 | This routines appends the value of num to the bits already present in
521 | the array buf. You need to give it the number of bits to use and you
522 | better make sure that this number of bits is enough to hold the value
523 | Also num must be positive.
524 |
525 */
527 static void sendbits(int buf[], int num_of_bits, int num) {
529 unsigned int cnt, lastbyte;
530 int lastbits;
531 unsigned char * cbuf;
533 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
534 cnt = (unsigned int) buf[0];
535 lastbits = buf[1];
536 lastbyte =(unsigned int) buf[2];
537 while (num_of_bits >= 8) {
538 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
539 cbuf[cnt++] = lastbyte >> lastbits;
540 num_of_bits -= 8;
541 }
542 if (num_of_bits > 0) {
543 lastbyte = (lastbyte << num_of_bits) | num;
544 lastbits += num_of_bits;
545 if (lastbits >= 8) {
546 lastbits -= 8;
547 cbuf[cnt++] = lastbyte >> lastbits;
548 }
549 }
550 buf[0] = cnt;
551 buf[1] = lastbits;
552 buf[2] = lastbyte;
553 if (lastbits>0) {
554 cbuf[cnt] = lastbyte << (8 - lastbits);
555 }
556 }
558 /*_________________________________________________________________________
559 |
560 | sizeofint - calculate bitsize of an integer
561 |
562 | return the number of bits needed to store an integer with given max size
563 |
564 */
566 static int sizeofint(const int size) {
567 unsigned int num = 1;
568 int num_of_bits = 0;
570 while (size >= num && num_of_bits < 32) {
571 num_of_bits++;
572 num <<= 1;
573 }
574 return num_of_bits;
575 }
577 /*___________________________________________________________________________
578 |
579 | sizeofints - calculate 'bitsize' of compressed ints
580 |
581 | given the number of small unsigned integers and the maximum value
582 | return the number of bits needed to read or write them with the
583 | routines receiveints and sendints. You need this parameter when
584 | calling these routines. Note that for many calls I can use
585 | the variable 'smallidx' which is exactly the number of bits, and
586 | So I don't need to call 'sizeofints for those calls.
587 */
589 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
590 int i, num;
591 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
592 num_of_bytes = 1;
593 bytes[0] = 1;
594 num_of_bits = 0;
595 for (i=0; i < num_of_ints; i++) {
596 tmp = 0;
597 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
598 tmp = bytes[bytecnt] * sizes[i] + tmp;
599 bytes[bytecnt] = tmp & 0xff;
600 tmp >>= 8;
601 }
602 while (tmp != 0) {
603 bytes[bytecnt++] = tmp & 0xff;
604 tmp >>= 8;
605 }
606 num_of_bytes = bytecnt;
607 }
608 num = 1;
609 num_of_bytes--;
610 while (bytes[num_of_bytes] >= num) {
611 num_of_bits++;
612 num *= 2;
613 }
614 return num_of_bits + num_of_bytes * 8;
616 }
618 /*____________________________________________________________________________
619 |
620 | sendints - send a small set of small integers in compressed format
621 |
622 | this routine is used internally by xdr3dfcoord, to send a set of
623 | small integers to the buffer.
624 | Multiplication with fixed (specified maximum ) sizes is used to get
625 | to one big, multibyte integer. Allthough the routine could be
626 | modified to handle sizes bigger than 16777216, or more than just
627 | a few integers, this is not done, because the gain in compression
628 | isn't worth the effort. Note that overflowing the multiplication
629 | or the byte buffer (32 bytes) is unchecked and causes bad results.
630 |
631 */
633 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
634 unsigned int sizes[], unsigned int nums[]) {
636 int i;
637 unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
639 tmp = nums[0];
640 num_of_bytes = 0;
641 do {
642 bytes[num_of_bytes++] = tmp & 0xff;
643 tmp >>= 8;
644 } while (tmp != 0);
646 for (i = 1; i < num_of_ints; i++) {
647 if (nums[i] >= sizes[i]) {
648 fprintf(stderr,"major breakdown in sendints num %u doesn't "
649 "match size %u\n", nums[i], sizes[i]);
650 exit(1);
651 }
652 /* use one step multiply */
653 tmp = nums[i];
654 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
655 tmp = bytes[bytecnt] * sizes[i] + tmp;
656 bytes[bytecnt] = tmp & 0xff;
657 tmp >>= 8;
658 }
659 while (tmp != 0) {
660 bytes[bytecnt++] = tmp & 0xff;
661 tmp >>= 8;
662 }
663 num_of_bytes = bytecnt;
664 }
665 if (num_of_bits >= num_of_bytes * 8) {
666 for (i = 0; i < num_of_bytes; i++) {
667 sendbits(buf, 8, bytes[i]);
668 }
669 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
670 } else {
671 for (i = 0; i < num_of_bytes-1; i++) {
672 sendbits(buf, 8, bytes[i]);
673 }
674 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
675 }
676 }
679 /*___________________________________________________________________________
680 |
681 | receivebits - decode number from buf using specified number of bits
682 |
683 | extract the number of bits from the array buf and construct an integer
684 | from it. Return that value.
685 |
686 */
688 static int receivebits(int buf[], int num_of_bits) {
690 int cnt, num;
691 unsigned int lastbits, lastbyte;
692 unsigned char * cbuf;
693 int mask = (1 << num_of_bits) -1;
695 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
696 cnt = buf[0];
697 lastbits = (unsigned int) buf[1];
698 lastbyte = (unsigned int) buf[2];
700 num = 0;
701 while (num_of_bits >= 8) {
702 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
703 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
704 num_of_bits -=8;
705 }
706 if (num_of_bits > 0) {
707 if (lastbits < num_of_bits) {
708 lastbits += 8;
709 lastbyte = (lastbyte << 8) | cbuf[cnt++];
710 }
711 lastbits -= num_of_bits;
712 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
713 }
714 num &= mask;
715 buf[0] = cnt;
716 buf[1] = lastbits;
717 buf[2] = lastbyte;
718 return num;
719 }
721 /*____________________________________________________________________________
722 |
723 | receiveints - decode 'small' integers from the buf array
724 |
725 | this routine is the inverse from sendints() and decodes the small integers
726 | written to buf by calculating the remainder and doing divisions with
727 | the given sizes[]. You need to specify the total number of bits to be
728 | used from buf in num_of_bits.
729 |
730 */
732 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
733 unsigned int sizes[], int nums[]) {
734 int bytes[32];
735 int i, j, num_of_bytes, p, num;
737 bytes[1] = bytes[2] = bytes[3] = 0;
738 num_of_bytes = 0;
739 while (num_of_bits > 8) {
740 bytes[num_of_bytes++] = receivebits(buf, 8);
741 num_of_bits -= 8;
742 }
743 if (num_of_bits > 0) {
744 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
745 }
746 for (i = num_of_ints-1; i > 0; i--) {
747 num = 0;
748 for (j = num_of_bytes-1; j >=0; j--) {
749 num = (num << 8) | bytes[j];
750 p = num / sizes[i];
751 bytes[j] = p;
752 num = num - p * sizes[i];
753 }
754 nums[i] = num;
755 }
756 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
757 }
759 /*____________________________________________________________________________
760 |
761 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
762 |
763 | this routine reads or writes (depending on how you opened the file with
764 | xdropen() ) a large number of 3d coordinates (stored in *fp).
765 | The number of coordinates triplets to write is given by *size. On
766 | read this number may be zero, in which case it reads as many as were written
767 | or it may specify the number if triplets to read (which should match the
768 | number written).
769 | Compression is achieved by first converting all floating numbers to integer
770 | using multiplication by *precision and rounding to the nearest integer.
771 | Then the minimum and maximum value are calculated to determine the range.
772 | The limited range of integers so found, is used to compress the coordinates.
773 | In addition the differences between succesive coordinates is calculated.
774 | If the difference happens to be 'small' then only the difference is saved,
775 | compressing the data even more. The notion of 'small' is changed dynamically
776 | and is enlarged or reduced whenever needed or possible.
777 | Extra compression is achieved in the case of GROMOS and coordinates of
778 | water molecules. GROMOS first writes out the Oxygen position, followed by
779 | the two hydrogens. In order to make the differences smaller (and thereby
780 | compression the data better) the order is changed into first one hydrogen
781 | then the oxygen, followed by the other hydrogen. This is rather special, but
782 | it shouldn't harm in the general case.
783 |
784 */
786 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
789 static int *ip = NULL;
790 static int oldsize;
791 static int *buf;
793 int minint[3], maxint[3], mindiff, *lip, diff;
794 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
795 int minidx, maxidx;
796 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
797 int flag, k;
798 int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
799 float *lfp, lf;
800 int tmp, *thiscoord, prevcoord[3];
801 unsigned int tmpcoord[30];
803 int bufsize, xdrid, lsize;
804 unsigned int bitsize;
805 float inv_precision;
806 int errval = 1;
808 /* find out if xdrs is opened for reading or for writing */
809 xdrid = 0;
810 while (xdridptr[xdrid] != xdrs) {
811 xdrid++;
812 if (xdrid >= MAXID) {
813 fprintf(stderr, "xdr error. no open xdr stream\n");
814 exit (1);
815 }
816 }
817 if ((xdrmodes[xdrid] == 'w') || (xdrmodes[xdrid] == 'a')) {
819 /* xdrs is open for writing */
821 if (xdr_int(xdrs, size) == 0)
822 return 0;
823 size3 = *size * 3;
824 /* when the number of coordinates is small, don't try to compress; just
825 * write them as floats using xdr_vector
826 */
827 if (*size <= 9 ) {
828 return (xdr_vector(xdrs, (char *) fp, (u_int)size3,
829 (u_int)sizeof(*fp), (xdrproc_t)xdr_float));
830 }
832 xdr_float(xdrs, precision);
833 if (ip == NULL) {
834 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
835 if (ip == NULL) {
836 fprintf(stderr,"malloc failed\n");
837 exit(1);
838 }
839 bufsize = size3 * 1.2;
840 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
841 if (buf == NULL) {
842 fprintf(stderr,"malloc failed\n");
843 exit(1);
844 }
845 oldsize = *size;
846 } else if (*size > oldsize) {
847 ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
848 if (ip == NULL) {
849 fprintf(stderr,"malloc failed\n");
850 exit(1);
851 }
852 bufsize = size3 * 1.2;
853 buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
854 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
855 if (buf == NULL) {
856 fprintf(stderr,"malloc failed\n");
857 exit(1);
858 }
859 oldsize = *size;
860 }
861 /* buf[0-2] are special and do not contain actual data */
862 buf[0] = buf[1] = buf[2] = 0;
863 minint[0] = minint[1] = minint[2] = INT_MAX;
864 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
865 prevrun = -1;
866 lfp = fp;
867 lip = ip;
868 mindiff = INT_MAX;
869 oldlint1 = oldlint2 = oldlint3 = 0;
870 while(lfp < fp + size3 ) {
871 /* find nearest integer */
872 if (*lfp >= 0.0)
873 lf = *lfp * *precision + 0.5;
874 else
875 lf = *lfp * *precision - 0.5;
876 if (fabs(lf) > MAXABS) {
877 /* scaling would cause overflow */
878 errval = 0;
879 }
880 lint1 = lf;
881 if (lint1 < minint[0]) minint[0] = lint1;
882 if (lint1 > maxint[0]) maxint[0] = lint1;
883 *lip++ = lint1;
884 lfp++;
885 if (*lfp >= 0.0)
886 lf = *lfp * *precision + 0.5;
887 else
888 lf = *lfp * *precision - 0.5;
889 if (fabs(lf) > MAXABS) {
890 /* scaling would cause overflow */
891 errval = 0;
892 }
893 lint2 = lf;
894 if (lint2 < minint[1]) minint[1] = lint2;
895 if (lint2 > maxint[1]) maxint[1] = lint2;
896 *lip++ = lint2;
897 lfp++;
898 if (*lfp >= 0.0)
899 lf = *lfp * *precision + 0.5;
900 else
901 lf = *lfp * *precision - 0.5;
902 if (fabs(lf) > MAXABS) {
903 /* scaling would cause overflow */
904 errval = 0;
905 }
906 lint3 = lf;
907 if (lint3 < minint[2]) minint[2] = lint3;
908 if (lint3 > maxint[2]) maxint[2] = lint3;
909 *lip++ = lint3;
910 lfp++;
911 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
912 if (diff < mindiff && lfp > fp + 3)
913 mindiff = diff;
914 oldlint1 = lint1;
915 oldlint2 = lint2;
916 oldlint3 = lint3;
917 }
918 xdr_int(xdrs, &(minint[0]));
919 xdr_int(xdrs, &(minint[1]));
920 xdr_int(xdrs, &(minint[2]));
922 xdr_int(xdrs, &(maxint[0]));
923 xdr_int(xdrs, &(maxint[1]));
924 xdr_int(xdrs, &(maxint[2]));
926 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
927 (float)maxint[1] - (float)minint[1] >= MAXABS ||
928 (float)maxint[2] - (float)minint[2] >= MAXABS) {
929 /* turning value in unsigned by subtracting minint
930 * would cause overflow
931 */
932 errval = 0;
933 }
934 sizeint[0] = maxint[0] - minint[0]+1;
935 sizeint[1] = maxint[1] - minint[1]+1;
936 sizeint[2] = maxint[2] - minint[2]+1;
938 /* check if one of the sizes is to big to be multiplied */
939 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
940 bitsizeint[0] = sizeofint(sizeint[0]);
941 bitsizeint[1] = sizeofint(sizeint[1]);
942 bitsizeint[2] = sizeofint(sizeint[2]);
943 bitsize = 0; /* flag the use of large sizes */
944 } else {
945 bitsize = sizeofints(3, sizeint);
946 }
947 lip = ip;
948 luip = (unsigned int *) ip;
949 smallidx = FIRSTIDX;
950 while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
951 smallidx++;
952 }
953 xdr_int(xdrs, &smallidx);
954 maxidx = MIN(LASTIDX, smallidx + 8) ;
955 minidx = maxidx - 8; /* often this equal smallidx */
956 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
957 small = magicints[smallidx] / 2;
958 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
959 larger = magicints[maxidx] / 2;
960 i = 0;
961 while (i < *size) {
962 is_small = 0;
963 thiscoord = (int *)(luip) + i * 3;
964 if (smallidx < maxidx && i >= 1 &&
965 abs(thiscoord[0] - prevcoord[0]) < larger &&
966 abs(thiscoord[1] - prevcoord[1]) < larger &&
967 abs(thiscoord[2] - prevcoord[2]) < larger) {
968 is_smaller = 1;
969 } else if (smallidx > minidx) {
970 is_smaller = -1;
971 } else {
972 is_smaller = 0;
973 }
974 if (i + 1 < *size) {
975 if (abs(thiscoord[0] - thiscoord[3]) < small &&
976 abs(thiscoord[1] - thiscoord[4]) < small &&
977 abs(thiscoord[2] - thiscoord[5]) < small) {
978 /* interchange first with second atom for better
979 * compression of water molecules
980 */
981 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
982 thiscoord[3] = tmp;
983 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
984 thiscoord[4] = tmp;
985 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
986 thiscoord[5] = tmp;
987 is_small = 1;
988 }
990 }
991 tmpcoord[0] = thiscoord[0] - minint[0];
992 tmpcoord[1] = thiscoord[1] - minint[1];
993 tmpcoord[2] = thiscoord[2] - minint[2];
994 if (bitsize == 0) {
995 sendbits(buf, bitsizeint[0], tmpcoord[0]);
996 sendbits(buf, bitsizeint[1], tmpcoord[1]);
997 sendbits(buf, bitsizeint[2], tmpcoord[2]);
998 } else {
999 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1000 }
1001 prevcoord[0] = thiscoord[0];
1002 prevcoord[1] = thiscoord[1];
1003 prevcoord[2] = thiscoord[2];
1004 thiscoord = thiscoord + 3;
1005 i++;
1007 run = 0;
1008 if (is_small == 0 && is_smaller == -1)
1009 is_smaller = 0;
1010 while (is_small && run < 8*3) {
1011 if (is_smaller == -1 && (
1012 SQR(thiscoord[0] - prevcoord[0]) +
1013 SQR(thiscoord[1] - prevcoord[1]) +
1014 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1015 is_smaller = 0;
1016 }
1018 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1019 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1020 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1022 prevcoord[0] = thiscoord[0];
1023 prevcoord[1] = thiscoord[1];
1024 prevcoord[2] = thiscoord[2];
1026 i++;
1027 thiscoord = thiscoord + 3;
1028 is_small = 0;
1029 if (i < *size &&
1030 abs(thiscoord[0] - prevcoord[0]) < small &&
1031 abs(thiscoord[1] - prevcoord[1]) < small &&
1032 abs(thiscoord[2] - prevcoord[2]) < small) {
1033 is_small = 1;
1034 }
1035 }
1036 if (run != prevrun || is_smaller != 0) {
1037 prevrun = run;
1038 sendbits(buf, 1, 1); /* flag the change in run-length */
1039 sendbits(buf, 5, run+is_smaller+1);
1040 } else {
1041 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1042 }
1043 for (k=0; k < run; k+=3) {
1044 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1045 }
1046 if (is_smaller != 0) {
1047 smallidx += is_smaller;
1048 if (is_smaller < 0) {
1049 small = smaller;
1050 smaller = magicints[smallidx-1] / 2;
1051 } else {
1052 smaller = small;
1053 small = magicints[smallidx] / 2;
1054 }
1055 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1056 }
1057 }
1058 if (buf[1] != 0) buf[0]++;
1059 xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1060 return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1061 } else {
1063 /* xdrs is open for reading */
1065 if (xdr_int(xdrs, &lsize) == 0)
1066 return 0;
1067 if (*size != 0 && lsize != *size) {
1068 fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1069 "%d arg vs %d in file", *size, lsize);
1070 }
1071 *size = lsize;
1072 size3 = *size * 3;
1073 if (*size <= 9) {
1074 return (xdr_vector(xdrs, (char *) fp, (u_int)size3,
1075 (u_int)sizeof(*fp), (xdrproc_t)xdr_float));
1076 }
1077 xdr_float(xdrs, precision);
1078 if (ip == NULL) {
1079 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
1080 if (ip == NULL) {
1081 fprintf(stderr,"malloc failed\n");
1082 exit(1);
1083 }
1084 bufsize = size3 * 1.2;
1085 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
1086 if (buf == NULL) {
1087 fprintf(stderr,"malloc failed\n");
1088 exit(1);
1089 }
1090 oldsize = *size;
1091 } else if (*size > oldsize) {
1092 ip = (int *)realloc(ip, (size_t)(size3 * sizeof(*ip)));
1093 if (ip == NULL) {
1094 fprintf(stderr,"malloc failed\n");
1095 exit(1);
1096 }
1097 bufsize = size3 * 1.2;
1098 buf = (int *)realloc(buf, (size_t)(bufsize * sizeof(*buf)));
1099 if (buf == NULL) {
1100 fprintf(stderr,"malloc failed\n");
1101 exit(1);
1102 }
1103 oldsize = *size;
1104 }
1105 buf[0] = buf[1] = buf[2] = 0;
1107 xdr_int(xdrs, &(minint[0]));
1108 xdr_int(xdrs, &(minint[1]));
1109 xdr_int(xdrs, &(minint[2]));
1111 xdr_int(xdrs, &(maxint[0]));
1112 xdr_int(xdrs, &(maxint[1]));
1113 xdr_int(xdrs, &(maxint[2]));
1115 sizeint[0] = maxint[0] - minint[0]+1;
1116 sizeint[1] = maxint[1] - minint[1]+1;
1117 sizeint[2] = maxint[2] - minint[2]+1;
1119 /* check if one of the sizes is to big to be multiplied */
1120 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1121 bitsizeint[0] = sizeofint(sizeint[0]);
1122 bitsizeint[1] = sizeofint(sizeint[1]);
1123 bitsizeint[2] = sizeofint(sizeint[2]);
1124 bitsize = 0; /* flag the use of large sizes */
1125 } else {
1126 bitsize = sizeofints(3, sizeint);
1127 }
1129 if (xdr_int(xdrs, &smallidx) == 0)
1130 return 0;
1131 maxidx = MIN(LASTIDX, smallidx + 8) ;
1132 minidx = maxidx - 8; /* often this equal smallidx */
1133 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1134 small = magicints[smallidx] / 2;
1135 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1136 larger = magicints[maxidx];
1138 /* buf[0] holds the length in bytes */
1140 if (xdr_int(xdrs, &(buf[0])) == 0)
1141 return 0;
1142 if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1143 return 0;
1144 buf[0] = buf[1] = buf[2] = 0;
1146 lfp = fp;
1147 inv_precision = 1.0 / * precision;
1148 run = 0;
1149 i = 0;
1150 lip = ip;
1151 while ( i < lsize ) {
1152 thiscoord = (int *)(lip) + i * 3;
1154 if (bitsize == 0) {
1155 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1156 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1157 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1158 } else {
1159 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1160 }
1162 i++;
1163 thiscoord[0] += minint[0];
1164 thiscoord[1] += minint[1];
1165 thiscoord[2] += minint[2];
1167 prevcoord[0] = thiscoord[0];
1168 prevcoord[1] = thiscoord[1];
1169 prevcoord[2] = thiscoord[2];
1172 flag = receivebits(buf, 1);
1173 is_smaller = 0;
1174 if (flag == 1) {
1175 run = receivebits(buf, 5);
1176 is_smaller = run % 3;
1177 run -= is_smaller;
1178 is_smaller--;
1179 }
1180 if (run > 0) {
1181 thiscoord += 3;
1182 for (k = 0; k < run; k+=3) {
1183 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1184 i++;
1185 thiscoord[0] += prevcoord[0] - small;
1186 thiscoord[1] += prevcoord[1] - small;
1187 thiscoord[2] += prevcoord[2] - small;
1188 if (k == 0) {
1189 /* interchange first with second atom for better
1190 * compression of water molecules
1191 */
1192 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1193 prevcoord[0] = tmp;
1194 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1195 prevcoord[1] = tmp;
1196 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1197 prevcoord[2] = tmp;
1198 *lfp++ = prevcoord[0] * inv_precision;
1199 *lfp++ = prevcoord[1] * inv_precision;
1200 *lfp++ = prevcoord[2] * inv_precision;
1201 } else {
1202 prevcoord[0] = thiscoord[0];
1203 prevcoord[1] = thiscoord[1];
1204 prevcoord[2] = thiscoord[2];
1205 }
1206 *lfp++ = thiscoord[0] * inv_precision;
1207 *lfp++ = thiscoord[1] * inv_precision;
1208 *lfp++ = thiscoord[2] * inv_precision;
1209 }
1210 } else {
1211 *lfp++ = thiscoord[0] * inv_precision;
1212 *lfp++ = thiscoord[1] * inv_precision;
1213 *lfp++ = thiscoord[2] * inv_precision;
1214 }
1215 smallidx += is_smaller;
1216 if (is_smaller < 0) {
1217 small = smaller;
1218 if (smallidx > FIRSTIDX) {
1219 smaller = magicints[smallidx - 1] /2;
1220 } else {
1221 smaller = 0;
1222 }
1223 } else if (is_smaller > 0) {
1224 smaller = small;
1225 small = magicints[smallidx] / 2;
1226 }
1227 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1228 }
1229 }
1230 return 1;
1231 }