Bug Summary

File:gromacs/fileio/libxdrf.c
Location:line 872, column 9
Description:Value stored to 'larger' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <limits.h>
42#include <math.h>
43#include <stdio.h>
44#include <stdlib.h>
45#include <string.h>
46
47#include "xdrf.h"
48#include "xdr_datatype.h"
49#include "gromacs/utility/futil.h"
50
51/* This is just for clarity - it can never be anything but 4! */
52#define XDR_INT_SIZE4 4
53
54/* same order as the definition of xdr_datatype */
55const char *xdr_datatype_names[] =
56{
57 "int",
58 "float",
59 "double",
60 "large int",
61 "char",
62 "string"
63};
64
65
66/*___________________________________________________________________________
67 |
68 | what follows are the C routine to read/write compressed coordinates together
69 | with some routines to assist in this task (those are marked
70 | static and cannot be called from user programs)
71 */
72#define MAXABS2147483647 -2 INT_MAX2147483647-2
73
74#ifndef MIN
75#define MIN(x, y)(((x)<(y))?(x):(y)) ((x) < (y) ? (x) : (y))
76#endif
77#ifndef MAX
78#define MAX(x, y)(((x)>(y))?(x):(y)) ((x) > (y) ? (x) : (y))
79#endif
80#ifndef SQR
81#define SQR(x)((x)*(x)) ((x)*(x))
82#endif
83static const int magicints[] = {
84 0, 0, 0, 0, 0, 0, 0, 0, 0,
85 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
86 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
87 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
88 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
89 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
90 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
91 8388607, 10568983, 13316085, 16777216
92};
93
94#define FIRSTIDX9 9
95/* note that magicints[FIRSTIDX-1] == 0 */
96#define LASTIDX(sizeof(magicints) / sizeof(*magicints)) (sizeof(magicints) / sizeof(*magicints))
97
98
99/*____________________________________________________________________________
100 |
101 | sendbits - encode num into buf using the specified number of bits
102 |
103 | This routines appends the value of num to the bits already present in
104 | the array buf. You need to give it the number of bits to use and you
105 | better make sure that this number of bits is enough to hold the value
106 | Also num must be positive.
107 |
108 */
109
110static void sendbits(int buf[], int num_of_bits, int num)
111{
112
113 unsigned int cnt, lastbyte;
114 int lastbits;
115 unsigned char * cbuf;
116
117 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
118 cnt = (unsigned int) buf[0];
119 lastbits = buf[1];
120 lastbyte = (unsigned int) buf[2];
121 while (num_of_bits >= 8)
122 {
123 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
124 cbuf[cnt++] = lastbyte >> lastbits;
125 num_of_bits -= 8;
126 }
127 if (num_of_bits > 0)
128 {
129 lastbyte = (lastbyte << num_of_bits) | num;
130 lastbits += num_of_bits;
131 if (lastbits >= 8)
132 {
133 lastbits -= 8;
134 cbuf[cnt++] = lastbyte >> lastbits;
135 }
136 }
137 buf[0] = cnt;
138 buf[1] = lastbits;
139 buf[2] = lastbyte;
140 if (lastbits > 0)
141 {
142 cbuf[cnt] = lastbyte << (8 - lastbits);
143 }
144}
145
146/*_________________________________________________________________________
147 |
148 | sizeofint - calculate bitsize of an integer
149 |
150 | return the number of bits needed to store an integer with given max size
151 |
152 */
153
154static int sizeofint(const int size)
155{
156 int num = 1;
157 int num_of_bits = 0;
158
159 while (size >= num && num_of_bits < 32)
160 {
161 num_of_bits++;
162 num <<= 1;
163 }
164 return num_of_bits;
165}
166
167/*___________________________________________________________________________
168 |
169 | sizeofints - calculate 'bitsize' of compressed ints
170 |
171 | given the number of small unsigned integers and the maximum value
172 | return the number of bits needed to read or write them with the
173 | routines receiveints and sendints. You need this parameter when
174 | calling these routines. Note that for many calls I can use
175 | the variable 'smallidx' which is exactly the number of bits, and
176 | So I don't need to call 'sizeofints for those calls.
177 */
178
179static int sizeofints( const int num_of_ints, unsigned int sizes[])
180{
181 int i, num;
182 int bytes[32];
183 unsigned int num_of_bytes, num_of_bits, bytecnt, tmp;
184 num_of_bytes = 1;
185 bytes[0] = 1;
186 num_of_bits = 0;
187 for (i = 0; i < num_of_ints; i++)
188 {
189 tmp = 0;
190 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
191 {
192 tmp = bytes[bytecnt] * sizes[i] + tmp;
193 bytes[bytecnt] = tmp & 0xff;
194 tmp >>= 8;
195 }
196 while (tmp != 0)
197 {
198 bytes[bytecnt++] = tmp & 0xff;
199 tmp >>= 8;
200 }
201 num_of_bytes = bytecnt;
202 }
203 num = 1;
204 num_of_bytes--;
205 while (bytes[num_of_bytes] >= num)
206 {
207 num_of_bits++;
208 num *= 2;
209 }
210 return num_of_bits + num_of_bytes * 8;
211
212}
213
214/*____________________________________________________________________________
215 |
216 | sendints - send a small set of small integers in compressed format
217 |
218 | this routine is used internally by xdr3dfcoord, to send a set of
219 | small integers to the buffer.
220 | Multiplication with fixed (specified maximum ) sizes is used to get
221 | to one big, multibyte integer. Allthough the routine could be
222 | modified to handle sizes bigger than 16777216, or more than just
223 | a few integers, this is not done, because the gain in compression
224 | isn't worth the effort. Note that overflowing the multiplication
225 | or the byte buffer (32 bytes) is unchecked and causes bad results.
226 |
227 */
228
229static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
230 unsigned int sizes[], unsigned int nums[])
231{
232
233 int i, num_of_bytes, bytecnt;
234 unsigned int bytes[32], tmp;
235
236 tmp = nums[0];
237 num_of_bytes = 0;
238 do
239 {
240 bytes[num_of_bytes++] = tmp & 0xff;
241 tmp >>= 8;
242 }
243 while (tmp != 0);
244
245 for (i = 1; i < num_of_ints; i++)
246 {
247 if (nums[i] >= sizes[i])
248 {
249 fprintf(stderrstderr, "major breakdown in sendints num %u doesn't "
250 "match size %u\n", nums[i], sizes[i]);
251 exit(1);
252 }
253 /* use one step multiply */
254 tmp = nums[i];
255 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++)
256 {
257 tmp = bytes[bytecnt] * sizes[i] + tmp;
258 bytes[bytecnt] = tmp & 0xff;
259 tmp >>= 8;
260 }
261 while (tmp != 0)
262 {
263 bytes[bytecnt++] = tmp & 0xff;
264 tmp >>= 8;
265 }
266 num_of_bytes = bytecnt;
267 }
268 if (num_of_bits >= num_of_bytes * 8)
269 {
270 for (i = 0; i < num_of_bytes; i++)
271 {
272 sendbits(buf, 8, bytes[i]);
273 }
274 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
275 }
276 else
277 {
278 for (i = 0; i < num_of_bytes-1; i++)
279 {
280 sendbits(buf, 8, bytes[i]);
281 }
282 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
283 }
284}
285
286
287/*___________________________________________________________________________
288 |
289 | receivebits - decode number from buf using specified number of bits
290 |
291 | extract the number of bits from the array buf and construct an integer
292 | from it. Return that value.
293 |
294 */
295
296static int receivebits(int buf[], int num_of_bits)
297{
298
299 int cnt, num, lastbits;
300 unsigned int lastbyte;
301 unsigned char * cbuf;
302 int mask = (1 << num_of_bits) -1;
303
304 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
305 cnt = buf[0];
306 lastbits = (unsigned int) buf[1];
307 lastbyte = (unsigned int) buf[2];
308
309 num = 0;
310 while (num_of_bits >= 8)
311 {
312 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
313 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
314 num_of_bits -= 8;
315 }
316 if (num_of_bits > 0)
317 {
318 if (lastbits < num_of_bits)
319 {
320 lastbits += 8;
321 lastbyte = (lastbyte << 8) | cbuf[cnt++];
322 }
323 lastbits -= num_of_bits;
324 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
325 }
326 num &= mask;
327 buf[0] = cnt;
328 buf[1] = lastbits;
329 buf[2] = lastbyte;
330 return num;
331}
332
333/*____________________________________________________________________________
334 |
335 | receiveints - decode 'small' integers from the buf array
336 |
337 | this routine is the inverse from sendints() and decodes the small integers
338 | written to buf by calculating the remainder and doing divisions with
339 | the given sizes[]. You need to specify the total number of bits to be
340 | used from buf in num_of_bits.
341 |
342 */
343
344static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
345 unsigned int sizes[], int nums[])
346{
347 int bytes[32];
348 int i, j, num_of_bytes, p, num;
349
350 bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
351 num_of_bytes = 0;
352 while (num_of_bits > 8)
353 {
354 bytes[num_of_bytes++] = receivebits(buf, 8);
355 num_of_bits -= 8;
356 }
357 if (num_of_bits > 0)
358 {
359 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
360 }
361 for (i = num_of_ints-1; i > 0; i--)
362 {
363 num = 0;
364 for (j = num_of_bytes-1; j >= 0; j--)
365 {
366 num = (num << 8) | bytes[j];
367 p = num / sizes[i];
368 bytes[j] = p;
369 num = num - p * sizes[i];
370 }
371 nums[i] = num;
372 }
373 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
374}
375
376/*____________________________________________________________________________
377 |
378 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
379 |
380 | this routine reads or writes (depending on how you opened the file with
381 | xdropen() ) a large number of 3d coordinates (stored in *fp).
382 | The number of coordinates triplets to write is given by *size. On
383 | read this number may be zero, in which case it reads as many as were written
384 | or it may specify the number if triplets to read (which should match the
385 | number written).
386 | Compression is achieved by first converting all floating numbers to integer
387 | using multiplication by *precision and rounding to the nearest integer.
388 | Then the minimum and maximum value are calculated to determine the range.
389 | The limited range of integers so found, is used to compress the coordinates.
390 | In addition the differences between succesive coordinates is calculated.
391 | If the difference happens to be 'small' then only the difference is saved,
392 | compressing the data even more. The notion of 'small' is changed dynamically
393 | and is enlarged or reduced whenever needed or possible.
394 | Extra compression is achieved in the case of GROMOS and coordinates of
395 | water molecules. GROMOS first writes out the Oxygen position, followed by
396 | the two hydrogens. In order to make the differences smaller (and thereby
397 | compression the data better) the order is changed into first one hydrogen
398 | then the oxygen, followed by the other hydrogen. This is rather special, but
399 | it shouldn't harm in the general case.
400 |
401 */
402
403int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
404{
405 int *ip = NULL((void*)0);
406 int *buf = NULL((void*)0);
407 gmx_bool bRead;
408
409 /* preallocate a small buffer and ip on the stack - if we need more
410 we can always malloc(). This is faster for small values of size: */
411 unsigned prealloc_size = 3*16;
412 int prealloc_ip[3*16], prealloc_buf[3*20];
413 int we_should_free = 0;
414
415 int minint[3], maxint[3], mindiff, *lip, diff;
416 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
417 int minidx, maxidx;
418 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
419 int flag, k;
420 int smallnum, smaller, larger, i, is_small, is_smaller, run, prevrun;
421 float *lfp, lf;
422 int tmp, *thiscoord, prevcoord[3];
423 unsigned int tmpcoord[30];
424
425 int bufsize, xdrid, lsize;
426 unsigned int bitsize;
427 float inv_precision;
428 int errval = 1;
429 int rc;
430
431 bRead = (xdrs->x_op == XDR_DECODE);
432 bitsizeint[0] = bitsizeint[1] = bitsizeint[2] = 0;
433 prevcoord[0] = prevcoord[1] = prevcoord[2] = 0;
434
435 if (!bRead)
436 {
437 /* xdrs is open for writing */
438
439 if (xdr_int(xdrs, size) == 0)
440 {
441 return 0;
442 }
443 size3 = *size * 3;
444 /* when the number of coordinates is small, don't try to compress; just
445 * write them as floats using xdr_vector
446 */
447 if (*size <= 9)
448 {
449 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
450 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
451 }
452
453 if (xdr_float(xdrs, precision) == 0)
454 {
455 return 0;
456 }
457
458 if (size3 <= prealloc_size)
459 {
460 ip = prealloc_ip;
461 buf = prealloc_buf;
462 }
463 else
464 {
465 we_should_free = 1;
466 bufsize = size3 * 1.2;
467 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
468 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
469 if (ip == NULL((void*)0) || buf == NULL((void*)0))
470 {
471 fprintf(stderrstderr, "malloc failed\n");
472 exit(1);
473 }
474 }
475 /* buf[0-2] are special and do not contain actual data */
476 buf[0] = buf[1] = buf[2] = 0;
477 minint[0] = minint[1] = minint[2] = INT_MAX2147483647;
478 maxint[0] = maxint[1] = maxint[2] = INT_MIN(-2147483647 -1);
479 prevrun = -1;
480 lfp = fp;
481 lip = ip;
482 mindiff = INT_MAX2147483647;
483 oldlint1 = oldlint2 = oldlint3 = 0;
484 while (lfp < fp + size3)
485 {
486 /* find nearest integer */
487 if (*lfp >= 0.0)
488 {
489 lf = *lfp * *precision + 0.5;
490 }
491 else
492 {
493 lf = *lfp * *precision - 0.5;
494 }
495 if (fabs(lf) > MAXABS2147483647 -2)
496 {
497 /* scaling would cause overflow */
498 errval = 0;
499 }
500 lint1 = lf;
501 if (lint1 < minint[0])
502 {
503 minint[0] = lint1;
504 }
505 if (lint1 > maxint[0])
506 {
507 maxint[0] = lint1;
508 }
509 *lip++ = lint1;
510 lfp++;
511 if (*lfp >= 0.0)
512 {
513 lf = *lfp * *precision + 0.5;
514 }
515 else
516 {
517 lf = *lfp * *precision - 0.5;
518 }
519 if (fabs(lf) > MAXABS2147483647 -2)
520 {
521 /* scaling would cause overflow */
522 errval = 0;
523 }
524 lint2 = lf;
525 if (lint2 < minint[1])
526 {
527 minint[1] = lint2;
528 }
529 if (lint2 > maxint[1])
530 {
531 maxint[1] = lint2;
532 }
533 *lip++ = lint2;
534 lfp++;
535 if (*lfp >= 0.0)
536 {
537 lf = *lfp * *precision + 0.5;
538 }
539 else
540 {
541 lf = *lfp * *precision - 0.5;
542 }
543 if (fabs(lf) > MAXABS2147483647 -2)
544 {
545 /* scaling would cause overflow */
546 errval = 0;
547 }
548 lint3 = lf;
549 if (lint3 < minint[2])
550 {
551 minint[2] = lint3;
552 }
553 if (lint3 > maxint[2])
554 {
555 maxint[2] = lint3;
556 }
557 *lip++ = lint3;
558 lfp++;
559 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
560 if (diff < mindiff && lfp > fp + 3)
561 {
562 mindiff = diff;
563 }
564 oldlint1 = lint1;
565 oldlint2 = lint2;
566 oldlint3 = lint3;
567 }
568 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
569 (xdr_int(xdrs, &(minint[1])) == 0) ||
570 (xdr_int(xdrs, &(minint[2])) == 0) ||
571 (xdr_int(xdrs, &(maxint[0])) == 0) ||
572 (xdr_int(xdrs, &(maxint[1])) == 0) ||
573 (xdr_int(xdrs, &(maxint[2])) == 0))
574 {
575 if (we_should_free)
576 {
577 free(ip);
578 free(buf);
579 }
580 return 0;
581 }
582
583 if ((float)maxint[0] - (float)minint[0] >= MAXABS2147483647 -2 ||
584 (float)maxint[1] - (float)minint[1] >= MAXABS2147483647 -2 ||
585 (float)maxint[2] - (float)minint[2] >= MAXABS2147483647 -2)
586 {
587 /* turning value in unsigned by subtracting minint
588 * would cause overflow
589 */
590 errval = 0;
591 }
592 sizeint[0] = maxint[0] - minint[0]+1;
593 sizeint[1] = maxint[1] - minint[1]+1;
594 sizeint[2] = maxint[2] - minint[2]+1;
595
596 /* check if one of the sizes is to big to be multiplied */
597 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
598 {
599 bitsizeint[0] = sizeofint(sizeint[0]);
600 bitsizeint[1] = sizeofint(sizeint[1]);
601 bitsizeint[2] = sizeofint(sizeint[2]);
602 bitsize = 0; /* flag the use of large sizes */
603 }
604 else
605 {
606 bitsize = sizeofints(3, sizeint);
607 }
608 lip = ip;
609 luip = (unsigned int *) ip;
610 smallidx = FIRSTIDX9;
611 while (smallidx < LASTIDX(sizeof(magicints) / sizeof(*magicints)) && magicints[smallidx] < mindiff)
612 {
613 smallidx++;
614 }
615 if (xdr_int(xdrs, &smallidx) == 0)
616 {
617 if (we_should_free)
618 {
619 free(ip);
620 free(buf);
621 }
622 return 0;
623 }
624
625 maxidx = MIN(LASTIDX, smallidx + 8)((((sizeof(magicints) / sizeof(*magicints)))<(smallidx + 8
))?((sizeof(magicints) / sizeof(*magicints))):(smallidx + 8));
626 minidx = maxidx - 8; /* often this equal smallidx */
627 smaller = magicints[MAX(FIRSTIDX, smallidx-1)(((9)>(smallidx-1))?(9):(smallidx-1))] / 2;
628 smallnum = magicints[smallidx] / 2;
629 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
630 larger = magicints[maxidx] / 2;
631 i = 0;
632 while (i < *size)
633 {
634 is_small = 0;
635 thiscoord = (int *)(luip) + i * 3;
636 if (smallidx < maxidx && i >= 1 &&
637 abs(thiscoord[0] - prevcoord[0]) < larger &&
638 abs(thiscoord[1] - prevcoord[1]) < larger &&
639 abs(thiscoord[2] - prevcoord[2]) < larger)
640 {
641 is_smaller = 1;
642 }
643 else if (smallidx > minidx)
644 {
645 is_smaller = -1;
646 }
647 else
648 {
649 is_smaller = 0;
650 }
651 if (i + 1 < *size)
652 {
653 if (abs(thiscoord[0] - thiscoord[3]) < smallnum &&
654 abs(thiscoord[1] - thiscoord[4]) < smallnum &&
655 abs(thiscoord[2] - thiscoord[5]) < smallnum)
656 {
657 /* interchange first with second atom for better
658 * compression of water molecules
659 */
660 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
661 thiscoord[3] = tmp;
662 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
663 thiscoord[4] = tmp;
664 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
665 thiscoord[5] = tmp;
666 is_small = 1;
667 }
668
669 }
670 tmpcoord[0] = thiscoord[0] - minint[0];
671 tmpcoord[1] = thiscoord[1] - minint[1];
672 tmpcoord[2] = thiscoord[2] - minint[2];
673 if (bitsize == 0)
674 {
675 sendbits(buf, bitsizeint[0], tmpcoord[0]);
676 sendbits(buf, bitsizeint[1], tmpcoord[1]);
677 sendbits(buf, bitsizeint[2], tmpcoord[2]);
678 }
679 else
680 {
681 sendints(buf, 3, bitsize, sizeint, tmpcoord);
682 }
683 prevcoord[0] = thiscoord[0];
684 prevcoord[1] = thiscoord[1];
685 prevcoord[2] = thiscoord[2];
686 thiscoord = thiscoord + 3;
687 i++;
688
689 run = 0;
690 if (is_small == 0 && is_smaller == -1)
691 {
692 is_smaller = 0;
693 }
694 while (is_small && run < 8*3)
695 {
696 if (is_smaller == -1 && (
697 SQR(thiscoord[0] - prevcoord[0])((thiscoord[0] - prevcoord[0])*(thiscoord[0] - prevcoord[0])) +
698 SQR(thiscoord[1] - prevcoord[1])((thiscoord[1] - prevcoord[1])*(thiscoord[1] - prevcoord[1])) +
699 SQR(thiscoord[2] - prevcoord[2])((thiscoord[2] - prevcoord[2])*(thiscoord[2] - prevcoord[2])) >= smaller * smaller))
700 {
701 is_smaller = 0;
702 }
703
704 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + smallnum;
705 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + smallnum;
706 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + smallnum;
707
708 prevcoord[0] = thiscoord[0];
709 prevcoord[1] = thiscoord[1];
710 prevcoord[2] = thiscoord[2];
711
712 i++;
713 thiscoord = thiscoord + 3;
714 is_small = 0;
715 if (i < *size &&
716 abs(thiscoord[0] - prevcoord[0]) < smallnum &&
717 abs(thiscoord[1] - prevcoord[1]) < smallnum &&
718 abs(thiscoord[2] - prevcoord[2]) < smallnum)
719 {
720 is_small = 1;
721 }
722 }
723 if (run != prevrun || is_smaller != 0)
724 {
725 prevrun = run;
726 sendbits(buf, 1, 1); /* flag the change in run-length */
727 sendbits(buf, 5, run+is_smaller+1);
728 }
729 else
730 {
731 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
732 }
733 for (k = 0; k < run; k += 3)
734 {
735 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
736 }
737 if (is_smaller != 0)
738 {
739 smallidx += is_smaller;
740 if (is_smaller < 0)
741 {
742 smallnum = smaller;
743 smaller = magicints[smallidx-1] / 2;
744 }
745 else
746 {
747 smaller = smallnum;
748 smallnum = magicints[smallidx] / 2;
749 }
750 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
751 }
752 }
753 if (buf[1] != 0)
754 {
755 buf[0]++;
756 }
757 /* buf[0] holds the length in bytes */
758 if (xdr_int(xdrs, &(buf[0])) == 0)
759 {
760 if (we_should_free)
761 {
762 free(ip);
763 free(buf);
764 }
765 return 0;
766 }
767
768
769 rc = errval * (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]));
770 if (we_should_free)
771 {
772 free(ip);
773 free(buf);
774 }
775 return rc;
776
777 }
778 else
779 {
780
781 /* xdrs is open for reading */
782
783 if (xdr_int(xdrs, &lsize) == 0)
784 {
785 return 0;
786 }
787 if (*size != 0 && lsize != *size)
788 {
789 fprintf(stderrstderr, "wrong number of coordinates in xdr3dfcoord; "
790 "%d arg vs %d in file", *size, lsize);
791 }
792 *size = lsize;
793 size3 = *size * 3;
794 if (*size <= 9)
795 {
796 *precision = -1;
797 return (xdr_vector(xdrs, (char *) fp, (unsigned int)size3,
798 (unsigned int)sizeof(*fp), (xdrproc_t)xdr_float));
799 }
800 if (xdr_float(xdrs, precision) == 0)
801 {
802 return 0;
803 }
804
805 if (size3 <= prealloc_size)
806 {
807 ip = prealloc_ip;
808 buf = prealloc_buf;
809 }
810 else
811 {
812 we_should_free = 1;
813 bufsize = size3 * 1.2;
814 ip = (int *)malloc((size_t)(size3 * sizeof(*ip)));
815 buf = (int *)malloc((size_t)(bufsize * sizeof(*buf)));
816 if (ip == NULL((void*)0) || buf == NULL((void*)0))
817 {
818 fprintf(stderrstderr, "malloc failed\n");
819 exit(1);
820 }
821 }
822
823 buf[0] = buf[1] = buf[2] = 0;
824
825 if ( (xdr_int(xdrs, &(minint[0])) == 0) ||
826 (xdr_int(xdrs, &(minint[1])) == 0) ||
827 (xdr_int(xdrs, &(minint[2])) == 0) ||
828 (xdr_int(xdrs, &(maxint[0])) == 0) ||
829 (xdr_int(xdrs, &(maxint[1])) == 0) ||
830 (xdr_int(xdrs, &(maxint[2])) == 0))
831 {
832 if (we_should_free)
833 {
834 free(ip);
835 free(buf);
836 }
837 return 0;
838 }
839
840 sizeint[0] = maxint[0] - minint[0]+1;
841 sizeint[1] = maxint[1] - minint[1]+1;
842 sizeint[2] = maxint[2] - minint[2]+1;
843
844 /* check if one of the sizes is to big to be multiplied */
845 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff)
846 {
847 bitsizeint[0] = sizeofint(sizeint[0]);
848 bitsizeint[1] = sizeofint(sizeint[1]);
849 bitsizeint[2] = sizeofint(sizeint[2]);
850 bitsize = 0; /* flag the use of large sizes */
851 }
852 else
853 {
854 bitsize = sizeofints(3, sizeint);
855 }
856
857 if (xdr_int(xdrs, &smallidx) == 0)
858 {
859 if (we_should_free)
860 {
861 free(ip);
862 free(buf);
863 }
864 return 0;
865 }
866
867 maxidx = MIN(LASTIDX, smallidx + 8)((((sizeof(magicints) / sizeof(*magicints)))<(smallidx + 8
))?((sizeof(magicints) / sizeof(*magicints))):(smallidx + 8));
868 minidx = maxidx - 8; /* often this equal smallidx */
869 smaller = magicints[MAX(FIRSTIDX, smallidx-1)(((9)>(smallidx-1))?(9):(smallidx-1))] / 2;
870 smallnum = magicints[smallidx] / 2;
871 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
872 larger = magicints[maxidx];
Value stored to 'larger' is never read
873
874 /* buf[0] holds the length in bytes */
875
876 if (xdr_int(xdrs, &(buf[0])) == 0)
877 {
878 if (we_should_free)
879 {
880 free(ip);
881 free(buf);
882 }
883 return 0;
884 }
885
886
887 if (xdr_opaque(xdrs, (char *)&(buf[3]), (unsigned int)buf[0]) == 0)
888 {
889 if (we_should_free)
890 {
891 free(ip);
892 free(buf);
893 }
894 return 0;
895 }
896
897
898
899 buf[0] = buf[1] = buf[2] = 0;
900
901 lfp = fp;
902 inv_precision = 1.0 / *precision;
903 run = 0;
904 i = 0;
905 lip = ip;
906 while (i < lsize)
907 {
908 thiscoord = (int *)(lip) + i * 3;
909
910 if (bitsize == 0)
911 {
912 thiscoord[0] = receivebits(buf, bitsizeint[0]);
913 thiscoord[1] = receivebits(buf, bitsizeint[1]);
914 thiscoord[2] = receivebits(buf, bitsizeint[2]);
915 }
916 else
917 {
918 receiveints(buf, 3, bitsize, sizeint, thiscoord);
919 }
920
921 i++;
922 thiscoord[0] += minint[0];
923 thiscoord[1] += minint[1];
924 thiscoord[2] += minint[2];
925
926 prevcoord[0] = thiscoord[0];
927 prevcoord[1] = thiscoord[1];
928 prevcoord[2] = thiscoord[2];
929
930
931 flag = receivebits(buf, 1);
932 is_smaller = 0;
933 if (flag == 1)
934 {
935 run = receivebits(buf, 5);
936 is_smaller = run % 3;
937 run -= is_smaller;
938 is_smaller--;
939 }
940 if (run > 0)
941 {
942 thiscoord += 3;
943 for (k = 0; k < run; k += 3)
944 {
945 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
946 i++;
947 thiscoord[0] += prevcoord[0] - smallnum;
948 thiscoord[1] += prevcoord[1] - smallnum;
949 thiscoord[2] += prevcoord[2] - smallnum;
950 if (k == 0)
951 {
952 /* interchange first with second atom for better
953 * compression of water molecules
954 */
955 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
956 prevcoord[0] = tmp;
957 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
958 prevcoord[1] = tmp;
959 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
960 prevcoord[2] = tmp;
961 *lfp++ = prevcoord[0] * inv_precision;
962 *lfp++ = prevcoord[1] * inv_precision;
963 *lfp++ = prevcoord[2] * inv_precision;
964 }
965 else
966 {
967 prevcoord[0] = thiscoord[0];
968 prevcoord[1] = thiscoord[1];
969 prevcoord[2] = thiscoord[2];
970 }
971 *lfp++ = thiscoord[0] * inv_precision;
972 *lfp++ = thiscoord[1] * inv_precision;
973 *lfp++ = thiscoord[2] * inv_precision;
974 }
975 }
976 else
977 {
978 *lfp++ = thiscoord[0] * inv_precision;
979 *lfp++ = thiscoord[1] * inv_precision;
980 *lfp++ = thiscoord[2] * inv_precision;
981 }
982 smallidx += is_smaller;
983 if (is_smaller < 0)
984 {
985 smallnum = smaller;
986 if (smallidx > FIRSTIDX9)
987 {
988 smaller = magicints[smallidx - 1] /2;
989 }
990 else
991 {
992 smaller = 0;
993 }
994 }
995 else if (is_smaller > 0)
996 {
997 smaller = smallnum;
998 smallnum = magicints[smallidx] / 2;
999 }
1000 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1001 }
1002 }
1003 if (we_should_free)
1004 {
1005 free(ip);
1006 free(buf);
1007 }
1008 return 1;
1009}
1010
1011
1012
1013/******************************************************************
1014
1015 XTC files have a relatively simple structure.
1016 They have a header of 16 bytes and the rest are
1017 the compressed coordinates of the files. Due to the
1018 compression 00 is not present in the coordinates.
1019 The first 4 bytes of the header are the magic number
1020 1995 (0x000007CB). If we find this number we are guaranteed
1021 to be in the header, due to the presence of so many zeros.
1022 The second 4 bytes are the number of atoms in the frame, and is
1023 assumed to be constant. The third 4 bytes are the frame number.
1024 The last 4 bytes are a floating point representation of the time.
1025
1026 ********************************************************************/
1027
1028/* Must match definition in xtcio.c */
1029#ifndef XTC_MAGIC1995
1030#define XTC_MAGIC1995 1995
1031#endif
1032
1033static const int header_size = 16;
1034
1035/* Check if we are at the header start.
1036 At the same time it will also read 1 int */
1037static int xtc_at_header_start(FILE *fp, XDR *xdrs,
1038 int natoms, int * timestep, float * time)
1039{
1040 int i_inp[3];
1041 float f_inp[10];
1042 int i;
1043 gmx_off_t off;
1044
1045
1046 if ((off = gmx_ftell(fp)) < 0)
1047 {
1048 return -1;
1049 }
1050 /* read magic natoms and timestep */
1051 for (i = 0; i < 3; i++)
1052 {
1053 if (!xdr_int(xdrs, &(i_inp[i])))
1054 {
1055 gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0);
1056 return -1;
1057 }
1058 }
1059 /* quick return */
1060 if (i_inp[0] != XTC_MAGIC1995)
1061 {
1062 if (gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0))
1063 {
1064 return -1;
1065 }
1066 return 0;
1067 }
1068 /* read time and box */
1069 for (i = 0; i < 10; i++)
1070 {
1071 if (!xdr_float(xdrs, &(f_inp[i])))
1072 {
1073 gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0);
1074 return -1;
1075 }
1076 }
1077 /* Make a rigourous check to see if we are in the beggining of a header
1078 Hopefully there are no ambiguous cases */
1079 /* This check makes use of the fact that the box matrix has 3 zeroes on the upper
1080 right triangle and that the first element must be nonzero unless the entire matrix is zero
1081 */
1082 if (i_inp[1] == natoms &&
1083 ((f_inp[1] != 0 && f_inp[6] == 0) ||
1084 (f_inp[1] == 0 && f_inp[5] == 0 && f_inp[9] == 0)))
1085 {
1086 if (gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0))
1087 {
1088 return -1;
1089 }
1090 *time = f_inp[0];
1091 *timestep = i_inp[2];
1092 return 1;
1093 }
1094 if (gmx_fseek(fp, off+XDR_INT_SIZE4, SEEK_SET0))
1095 {
1096 return -1;
1097 }
1098 return 0;
1099}
1100
1101static int
1102xtc_get_next_frame_number(FILE *fp, XDR *xdrs, int natoms)
1103{
1104 gmx_off_t off;
1105 int step;
1106 float time;
1107 int ret;
1108
1109 if ((off = gmx_ftell(fp)) < 0)
1110 {
1111 return -1;
1112 }
1113
1114 /* read one int just to make sure we dont read this frame but the next */
1115 xdr_int(xdrs, &step);
1116 while (1)
1117 {
1118 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1119 if (ret == 1)
1120 {
1121 if (gmx_fseek(fp, off, SEEK_SET0))
1122 {
1123 return -1;
1124 }
1125 return step;
1126 }
1127 else if (ret == -1)
1128 {
1129 if (gmx_fseek(fp, off, SEEK_SET0))
1130 {
1131 return -1;
1132 }
1133 }
1134 }
1135 return -1;
1136}
1137
1138
1139static float xtc_get_next_frame_time(FILE *fp, XDR *xdrs, int natoms,
1140 gmx_bool * bOK)
1141{
1142 gmx_off_t off;
1143 float time;
1144 int step;
1145 int ret;
1146 *bOK = 0;
1147
1148 if ((off = gmx_ftell(fp)) < 0)
1149 {
1150 return -1;
1151 }
1152 /* read one int just to make sure we dont read this frame but the next */
1153 xdr_int(xdrs, &step);
1154 while (1)
1155 {
1156 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1157 if (ret == 1)
1158 {
1159 *bOK = 1;
1160 if (gmx_fseek(fp, off, SEEK_SET0))
1161 {
1162 *bOK = 0;
1163 return -1;
1164 }
1165 return time;
1166 }
1167 else if (ret == -1)
1168 {
1169 if (gmx_fseek(fp, off, SEEK_SET0))
1170 {
1171 return -1;
1172 }
1173 return -1;
1174 }
1175 }
1176 return -1;
1177}
1178
1179
1180static float
1181xtc_get_current_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1182{
1183 gmx_off_t off;
1184 int step;
1185 float time;
1186 int ret;
1187 *bOK = 0;
1188
1189 if ((off = gmx_ftell(fp)) < 0)
1190 {
1191 return -1;
1192 }
1193
1194 while (1)
1195 {
1196 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1197 if (ret == 1)
1198 {
1199 *bOK = 1;
1200 if (gmx_fseek(fp, off, SEEK_SET0))
1201 {
1202 *bOK = 0;
1203 return -1;
1204 }
1205 return time;
1206 }
1207 else if (ret == -1)
1208 {
1209 if (gmx_fseek(fp, off, SEEK_SET0))
1210 {
1211 return -1;
1212 }
1213 return -1;
1214 }
1215 else if (ret == 0)
1216 {
1217 /*Go back.*/
1218 if (gmx_fseek(fp, -2*XDR_INT_SIZE4, SEEK_CUR1))
1219 {
1220 return -1;
1221 }
1222 }
1223 }
1224 return -1;
1225}
1226
1227/* Currently not used, just for completeness */
1228static int
1229xtc_get_current_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1230{
1231 gmx_off_t off;
1232 int ret;
1233 int step;
1234 float time;
1235 *bOK = 0;
1236
1237 if ((off = gmx_ftell(fp)) < 0)
1238 {
1239 return -1;
1240 }
1241
1242
1243 while (1)
1244 {
1245 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1246 if (ret == 1)
1247 {
1248 *bOK = 1;
1249 if (gmx_fseek(fp, off, SEEK_SET0))
1250 {
1251 *bOK = 0;
1252 return -1;
1253 }
1254 return step;
1255 }
1256 else if (ret == -1)
1257 {
1258 if (gmx_fseek(fp, off, SEEK_SET0))
1259 {
1260 return -1;
1261 }
1262 return -1;
1263
1264 }
1265 else if (ret == 0)
1266 {
1267 /*Go back.*/
1268 if (gmx_fseek(fp, -2*XDR_INT_SIZE4, SEEK_CUR1))
1269 {
1270 return -1;
1271 }
1272 }
1273 }
1274 return -1;
1275}
1276
1277
1278
1279static gmx_off_t xtc_get_next_frame_start(FILE *fp, XDR *xdrs, int natoms)
1280{
1281 int inp;
1282 gmx_off_t res;
1283 int ret;
1284 int step;
1285 float time;
1286 /* read one int just to make sure we dont read this frame but the next */
1287 xdr_int(xdrs, &step);
1288 while (1)
1289 {
1290 ret = xtc_at_header_start(fp, xdrs, natoms, &step, &time);
1291 if (ret == 1)
1292 {
1293 if ((res = gmx_ftell(fp)) >= 0)
1294 {
1295 return res - XDR_INT_SIZE4;
1296 }
1297 else
1298 {
1299 return res;
1300 }
1301 }
1302 else if (ret == -1)
1303 {
1304 return -1;
1305 }
1306 }
1307 return -1;
1308}
1309
1310
1311static
1312float
1313xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1314{
1315 float res;
1316 float tinit;
1317 gmx_off_t off;
1318
1319 *bOK = 0;
1320 if ((off = gmx_ftell(fp)) < 0)
1321 {
1322 return -1;
1323 }
1324
1325 tinit = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1326
1327 if (!(*bOK))
1328 {
1329 return -1;
1330 }
1331
1332 res = xtc_get_next_frame_time(fp, xdrs, natoms, bOK);
1333
1334 if (!(*bOK))
1335 {
1336 return -1;
1337 }
1338
1339 res -= tinit;
1340 if (0 != gmx_fseek(fp, off, SEEK_SET0))
1341 {
1342 *bOK = 0;
1343 return -1;
1344 }
1345 return res;
1346}
1347
1348
1349int
1350xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms)
1351{
1352 gmx_off_t low = 0;
1353 gmx_off_t high, pos;
1354
1355
1356 /* round to 4 bytes */
1357 int fr;
1358 gmx_off_t offset;
1359 if (gmx_fseek(fp, 0, SEEK_END2))
1360 {
1361 return -1;
1362 }
1363
1364 if ((high = gmx_ftell(fp)) < 0)
1365 {
1366 return -1;
1367 }
1368
1369 /* round to 4 bytes */
1370 high /= XDR_INT_SIZE4;
1371 high *= XDR_INT_SIZE4;
1372 offset = ((high/2)/XDR_INT_SIZE4)*XDR_INT_SIZE4;
1373
1374 if (gmx_fseek(fp, offset, SEEK_SET0))
1375 {
1376 return -1;
1377 }
1378
1379 while (1)
1380 {
1381 fr = xtc_get_next_frame_number(fp, xdrs, natoms);
1382 if (fr < 0)
1383 {
1384 return -1;
1385 }
1386 if (fr != frame && abs(low-high) > header_size)
1387 {
1388 if (fr < frame)
1389 {
1390 low = offset;
1391 }
1392 else
1393 {
1394 high = offset;
1395 }
1396 /* round to 4 bytes */
1397 offset = (((high+low)/2)/4)*4;
1398
1399 if (gmx_fseek(fp, offset, SEEK_SET0))
1400 {
1401 return -1;
1402 }
1403 }
1404 else
1405 {
1406 break;
1407 }
1408 }
1409 if (offset <= header_size)
1410 {
1411 offset = low;
1412 }
1413
1414 if (gmx_fseek(fp, offset, SEEK_SET0))
1415 {
1416 return -1;
1417 }
1418
1419 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1420 {
1421 /* we probably hit an end of file */
1422 return -1;
1423 }
1424
1425 if (gmx_fseek(fp, pos, SEEK_SET0))
1426 {
1427 return -1;
1428 }
1429
1430 return 0;
1431}
1432
1433
1434
1435int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms, gmx_bool bSeekForwardOnly)
1436{
1437 float t;
1438 float dt;
1439 gmx_bool bOK = FALSE0;
1440 gmx_off_t low = 0;
1441 gmx_off_t high, offset, pos;
1442 int res;
1443 int dt_sign = 0;
1444
1445 if (bSeekForwardOnly)
1446 {
1447 low = gmx_ftell(fp);
1448 }
1449 if (gmx_fseek(fp, 0, SEEK_END2))
1450 {
1451 return -1;
1452 }
1453
1454 if ((high = gmx_ftell(fp)) < 0)
1455 {
1456 return -1;
1457 }
1458 /* round to int */
1459 high /= XDR_INT_SIZE4;
1460 high *= XDR_INT_SIZE4;
1461 offset = (((high-low) / 2) / XDR_INT_SIZE4) * XDR_INT_SIZE4;
1462
1463 if (gmx_fseek(fp, offset, SEEK_SET0))
1464 {
1465 return -1;
1466 }
1467
1468
1469 /*
1470 * No need to call xdr_xtc_estimate_dt here - since xdr_xtc_estimate_dt is called first thing in the loop
1471 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1472
1473 if (!bOK)
1474 {
1475 return -1;
1476 }
1477 */
1478
1479 while (1)
1480 {
1481 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1482 if (!bOK)
1483 {
1484 return -1;
1485 }
1486 else
1487 {
1488 if (dt > 0)
1489 {
1490 if (dt_sign == -1)
1491 {
1492 /* Found a place in the trajectory that has positive time step while
1493 other has negative time step */
1494 return -2;
1495 }
1496 dt_sign = 1;
1497 }
1498 else if (dt < 0)
1499 {
1500 if (dt_sign == 1)
1501 {
1502 /* Found a place in the trajectory that has positive time step while
1503 other has negative time step */
1504 return -2;
1505 }
1506 dt_sign = -1;
1507 }
1508 }
1509 t = xtc_get_next_frame_time(fp, xdrs, natoms, &bOK);
1510 if (!bOK)
1511 {
1512 return -1;
1513 }
1514
1515 /* If we are before the target time and the time step is positive or 0, or we have
1516 after the target time and the time step is negative, or the difference between
1517 the current time and the target time is bigger than dt and above all the distance between high
1518 and low is bigger than 1 frame, then do another step of binary search. Otherwise stop and check
1519 if we reached the solution */
1520 if ((((t < time && dt_sign >= 0) || (t > time && dt_sign == -1)) || ((t
1521 - time) >= dt && dt_sign >= 0)
1522 || ((time - t) >= -dt && dt_sign < 0)) && (abs(low - high)
1523 > header_size))
1524 {
1525 if (dt >= 0 && dt_sign != -1)
1526 {
1527 if (t < time)
1528 {
1529 low = offset;
1530 }
1531 else
1532 {
1533 high = offset;
1534 }
1535 }
1536 else if (dt <= 0 && dt_sign == -1)
1537 {
1538 if (t >= time)
1539 {
1540 low = offset;
1541 }
1542 else
1543 {
1544 high = offset;
1545 }
1546 }
1547 else
1548 {
1549 /* We should never reach here */
1550 return -1;
1551 }
1552 /* round to 4 bytes and subtract header*/
1553 offset = (((high + low) / 2) / XDR_INT_SIZE4) * XDR_INT_SIZE4;
1554 if (gmx_fseek(fp, offset, SEEK_SET0))
1555 {
1556 return -1;
1557 }
1558 }
1559 else
1560 {
1561 if (abs(low - high) <= header_size)
1562 {
1563 break;
1564 }
1565 /* re-estimate dt */
1566 if (xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK) != dt)
1567 {
1568 if (bOK)
1569 {
1570 dt = xdr_xtc_estimate_dt(fp, xdrs, natoms, &bOK);
1571 }
1572 }
1573 if (t >= time && t - time < dt)
1574 {
1575 break;
1576 }
1577 }
1578 }
1579
1580 if (offset <= header_size)
1581 {
1582 offset = low;
1583 }
1584
1585 gmx_fseek(fp, offset, SEEK_SET0);
1586
1587 if ((pos = xtc_get_next_frame_start(fp, xdrs, natoms)) < 0)
1588 {
1589 return -1;
1590 }
1591
1592 if (gmx_fseek(fp, pos, SEEK_SET0))
1593 {
1594 return -1;
1595 }
1596 return 0;
1597}
1598
1599float
1600xdr_xtc_get_last_frame_time(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1601{
1602 float time;
1603 gmx_off_t off;
1604 int res;
1605 *bOK = 1;
1606 off = gmx_ftell(fp);
1607 if (off < 0)
1608 {
1609 *bOK = 0;
1610 return -1;
1611 }
1612
1613 if ( (res = gmx_fseek(fp, -3*XDR_INT_SIZE4, SEEK_END2)) != 0)
1614 {
1615 *bOK = 0;
1616 return -1;
1617 }
1618
1619 time = xtc_get_current_frame_time(fp, xdrs, natoms, bOK);
1620 if (!(*bOK))
1621 {
1622 return -1;
1623 }
1624
1625 if ( (res = gmx_fseek(fp, off, SEEK_SET0)) != 0)
1626 {
1627 *bOK = 0;
1628 return -1;
1629 }
1630 return time;
1631}
1632
1633
1634int
1635xdr_xtc_get_last_frame_number(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK)
1636{
1637 int frame;
1638 gmx_off_t off;
1639 int res;
1640 *bOK = 1;
1641
1642 if ((off = gmx_ftell(fp)) < 0)
1643 {
1644 *bOK = 0;
1645 return -1;
1646 }
1647
1648
1649 if (gmx_fseek(fp, -3*XDR_INT_SIZE4, SEEK_END2))
1650 {
1651 *bOK = 0;
1652 return -1;
1653 }
1654
1655 frame = xtc_get_current_frame_number(fp, xdrs, natoms, bOK);
1656 if (!bOK)
1657 {
1658 return -1;
1659 }
1660
1661
1662 if (gmx_fseek(fp, off, SEEK_SET0))
1663 {
1664 *bOK = 0;
1665 return -1;
1666 }
1667
1668 return frame;
1669}