?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Update TNG build system
[gromacs.git] / src / external / tng_io / src / tests / compression / testsuite.c
blobe266c110b8b43995e56ef63ab34f165dfe163eff
1 /* tng compression routines */
2
3 /* Only modify testsuite.c
4 *Then* run testsuite.sh to perform the test.
5 */
6
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <string.h>
10 #include <math.h>
11 #include "compression/tng_compress.h"
12 #include TESTPARAM
13
14 #ifdef TEST_FLOAT
15 #define REAL float
16 #else
17 #define REAL double
18 #endif
19
20 #ifndef TNG_COMPRESS_FILES_DIR
21 #define TNG_COMPRESS_FILES_DIR ""
22 #endif
23
24 #define FUDGE 1.1 /* 10% off target precision is acceptable */
25
26 static void keepinbox(int *val)
27 {
28 while (val[0]>INTMAX1)
29 val[0]-=(INTMAX1-INTMIN1+1);
30 while (val[0]<INTMIN1)
31 val[0]+=(INTMAX1-INTMIN1+1);
32 while (val[1]>INTMAX2)
33 val[1]-=(INTMAX2-INTMIN2+1);
34 while (val[1]<INTMIN2)
35 val[1]+=(INTMAX2-INTMIN2+1);
36 while (val[2]>INTMAX3)
37 val[2]-=(INTMAX3-INTMIN3+1);
38 while (val[2]<INTMIN3)
39 val[2]+=(INTMAX3-INTMIN3+1);
40 }
41
42 static int intsintable[128]={
43 0 , 3215 , 6423 , 9615 , 12785 , 15923 , 19023 , 22078 ,
44 25079 , 28019 , 30892 , 33691 , 36409 , 39039 , 41574 , 44010 ,
45 46340 , 48558 , 50659 , 52638 , 54490 , 56211 , 57796 , 59242 ,
46 60546 , 61704 , 62713 , 63570 , 64275 , 64825 , 65219 , 65456 ,
47 65535 , 65456 , 65219 , 64825 , 64275 , 63570 , 62713 , 61704 ,
48 60546 , 59242 , 57796 , 56211 , 54490 , 52638 , 50659 , 48558 ,
49 46340 , 44010 , 41574 , 39039 , 36409 , 33691 , 30892 , 28019 ,
50 25079 , 22078 , 19023 , 15923 , 12785 , 9615 , 6423 , 3215 ,
51 0 , -3215 , -6423 , -9615 , -12785 , -15923 , -19023 , -22078 ,
52 -25079 , -28019 , -30892 , -33691 , -36409 , -39039 , -41574 , -44010 ,
53 -46340 , -48558 , -50659 , -52638 , -54490 , -56211 , -57796 , -59242 ,
54 -60546 , -61704 , -62713 , -63570 , -64275 , -64825 , -65219 , -65456 ,
55 -65535 , -65456 , -65219 , -64825 , -64275 , -63570 , -62713 , -61704 ,
56 -60546 , -59242 , -57796 , -56211 , -54490 , -52638 , -50659 , -48558 ,
57 -46340 , -44010 , -41574 , -39039 , -36409 , -33691 , -30892 , -28019 ,
58 -25079 , -22078 , -19023 , -15923 , -12785 , -9615 , -6423 , -3215 ,
59 };
60
61 static int intsin(int i)
62 {
63 int sign=1;
64 if (i<0)
65 {
66 i=0;
67 sign=-1;
68 }
69 return sign*intsintable[i%128];
70 }
71
72 static int intcos(int i)
73 {
74 if (i<0)
75 i=0;
76 return intsin(i+32);
77 }
78
79 static void molecule(int *target,
80 int *base,
81 int length,
82 int scale, int *direction,
83 int flip,
84 int iframe)
85 {
86 int i;
87 for (i=0; i<length; i++)
88 {
89 int ifl=i;
90 if ((i==0) && (flip) && (length>1))
91 ifl=1;
92 else if ((i==1) && (flip) && (length>1))
93 ifl=0;
94 target[ifl*3]=base[0]+(intsin((i+iframe)*direction[0])*scale)/256;
95 target[ifl*3+1]=base[1]+(intcos((i+iframe)*direction[1])*scale)/256;
96 target[ifl*3+2]=base[2]+(intcos((i+iframe)*direction[2])*scale)/256;
97 keepinbox(target+ifl*3);
98 }
99 }
100
101 #ifndef FRAMESCALE
102 #define FRAMESCALE 1
103 #endif
104
105 static void genibox(int *intbox, int iframe)
106 {
107 int molecule_length=1;
108 int molpos[3];
109 int direction[3]={1,1,1};
110 int scale=1;
111 int flip=0;
112 int i=0;
113 molpos[0]=intsin(iframe*FRAMESCALE)/32;
114 molpos[1]=1+intcos(iframe*FRAMESCALE)/32;
115 molpos[2]=2+intsin(iframe*FRAMESCALE)/16;
116 keepinbox(molpos);
117 while (i<NATOMS)
118 {
119 int this_mol_length=molecule_length;
120 int dir;
121 #ifdef REGULAR
122 this_mol_length=4;
123 flip=0;
124 scale=1;
125 #endif
126 if (i+this_mol_length>NATOMS)
127 this_mol_length=NATOMS-i;
128 /* We must test the large rle as well. This requires special
129 sequencies to get triggered. So insert these from time to
130 time */
131 #ifndef REGULAR
132 if ((i%10)==0)
133 {
134 int j;
135 intbox[i*3]=molpos[0];
136 intbox[i*3+1]=molpos[1];
137 intbox[i*3+2]=molpos[2];
138 for (j=1; j<this_mol_length; j++)
139 {
140 intbox[(i+j)*3]=intbox[(i+j-1)*3]+(INTMAX1-INTMIN1+1)/5;
141 intbox[(i+j)*3+1]=intbox[(i+j-1)*3+1]+(INTMAX2-INTMIN2+1)/5;
142 intbox[(i+j)*3+2]=intbox[(i+j-1)*3+2]+(INTMAX3-INTMIN3+1)/5;
143 keepinbox(intbox+(i+j)*3);
144 }
145 }
146 else
147 #endif
148 molecule(intbox+i*3,molpos,this_mol_length,scale,direction,flip,iframe*FRAMESCALE);
149 i+=this_mol_length;
150 dir=1;
151 if (intsin(i*3)<0)
152 dir=-1;
153 molpos[0]+=dir*(INTMAX1-INTMIN1+1)/20;
154 dir=1;
155 if (intsin(i*5)<0)
156 dir=-1;
157 molpos[1]+=dir*(INTMAX2-INTMIN2+1)/20;
158 dir=1;
159 if (intsin(i*7)<0)
160 dir=-1;
161 molpos[2]+=dir*(INTMAX3-INTMIN3+1)/20;
162 keepinbox(molpos);
163
164 direction[0]=((direction[0]+1)%7)+1;
165 direction[1]=((direction[1]+1)%3)+1;
166 direction[2]=((direction[2]+1)%6)+1;
167
168 scale++;
169 if (scale>5)
170 scale=1;
171
172 molecule_length++;
173 if (molecule_length>30)
174 molecule_length=1;
175 if (i%9)
176 flip=1-flip;
177 }
178 }
179
180 static void genivelbox(int *intvelbox, int iframe)
181 {
182 int i;
183 for (i=0; i<NATOMS; i++)
184 {
185 #ifdef VELINTMUL
186 intvelbox[i*3]=((intsin((i+iframe*FRAMESCALE)*3))/10)*VELINTMUL+i;
187 intvelbox[i*3+1]=1+((intcos((i+iframe*FRAMESCALE)*5))/10)*VELINTMUL+i;
188 intvelbox[i*3+2]=2+((intsin((i+iframe*FRAMESCALE)*7)+intcos((i+iframe*FRAMESCALE)*9))/20)*VELINTMUL+i;
189 #else
190 intvelbox[i*3]=((intsin((i+iframe*FRAMESCALE)*3))/10);
191 intvelbox[i*3+1]=1+((intcos((i+iframe*FRAMESCALE)*5))/10);
192 intvelbox[i*3+2]=2+((intsin((i+iframe*FRAMESCALE)*7)+intcos((i+iframe*FRAMESCALE)*9))/20);
193 #endif
194 }
195 }
196
197 #ifndef STRIDE1
198 #define STRIDE1 3
199 #endif
200
201 #ifndef STRIDE2
202 #define STRIDE2 3
203 #endif
204
205 #ifndef GENPRECISION
206 #define GENPRECISION PRECISION
207 #endif
208
209 #ifndef GENVELPRECISION
210 #define GENVELPRECISION VELPRECISION
211 #endif
212
213 static void realbox(int *intbox, REAL *realbox, int stride)
214 {
215 int i,j;
216 for (i=0; i<NATOMS; i++)
217 {
218 for (j=0; j<3; j++)
219 realbox[i*stride+j]=(REAL)(intbox[i*3+j]*GENPRECISION*SCALE);
220 for (j=3; j<stride; j++)
221 realbox[i*stride+j]=0.;
222 }
223 }
224
225 static void realvelbox(int *intbox, REAL *realbox, int stride)
226 {
227 int i,j;
228 for (i=0; i<NATOMS; i++)
229 {
230 for (j=0; j<3; j++)
231 realbox[i*stride+j]=(REAL)(intbox[i*3+j]*GENVELPRECISION*SCALE);
232 for (j=3; j<stride; j++)
233 realbox[i*stride+j]=0.;
234 }
235 }
236
237 static int equalarr(REAL *arr1, REAL *arr2, REAL prec, int len, int itemlen, int stride1, int stride2)
238 {
239 REAL maxdiff=0.;
240 int i,j;
241 for (i=0; i<len; i++)
242 {
243 for (j=0; j<itemlen; j++)
244 if (fabs(arr1[i*stride1+j]-arr2[i*stride2+j])>maxdiff)
245 maxdiff=(REAL)fabs(arr1[i*stride1+j]-arr2[i*stride2+j]);
246 }
247 #if 0
248 for (i=0; i<len; i++)
249 {
250 for (j=0; j<itemlen; j++)
251 printf("%d %d: %g %g\n",i,j,arr1[i*stride1+j],arr2[i*stride2+j]);
252 }
253 #endif
254 #if 0
255 fprintf(stderr,"Error is %g. Acceptable error is %g.\n",maxdiff,prec*0.5*FUDGE);
256 #endif
257 if (maxdiff>prec*0.5*FUDGE)
258 {
259 return 0;
260 }
261 else
262 return 1;
263 }
264
265 struct tng_file
266 {
267 FILE *f;
268 int natoms;
269 int chunky;
270 REAL precision;
271 REAL velprecision;
272 int initial_coding;
273 int initial_coding_parameter;
274 int coding;
275 int coding_parameter;
276 int initial_velcoding;
277 int initial_velcoding_parameter;
278 int velcoding;
279 int velcoding_parameter;
280 int speed;
281 int nframes;
282 int nframes_delivered;
283 int writevel;
284 REAL *pos;
285 REAL *vel;
286 int *ipos;
287 int *ivel;
288 unsigned int prec_hi, prec_lo;
289 unsigned int velprec_hi, velprec_lo;
290 };
291
292 static size_t fwrite_int_le(int *x,FILE *f)
293 {
294 unsigned char c[4];
295 unsigned int i=(unsigned int)*x;
296 c[0]=(unsigned char)(i&0xFFU);
297 c[1]=(unsigned char)((i>>8)&0xFFU);
298 c[2]=(unsigned char)((i>>16)&0xFFU);
299 c[3]=(unsigned char)((i>>24)&0xFFU);
300 return fwrite(c,1,4,f);
301 }
302
303 static size_t fread_int_le(int *x,FILE *f)
304 {
305 unsigned char c[4];
306 unsigned int i;
307 size_t n=fread(c,1,4,f);
308 if (n)
309 {
310 i=(((unsigned int)c[3])<<24)|(((unsigned int)c[2])<<16)|(((unsigned int)c[1])<<8)|((unsigned int)c[0]);
311 *x=(int)i;
312 }
313 return n;
314 }
315
316 static struct tng_file *open_tng_file_write(char *filename,
317 int natoms,int chunky,
318 REAL precision,
319 int writevel,
320 REAL velprecision,
321 int initial_coding,
322 int initial_coding_parameter,
323 int coding,
324 int coding_parameter,
325 int initial_velcoding,
326 int initial_velcoding_parameter,
327 int velcoding,
328 int velcoding_parameter,
329 int speed)
330 {
331 struct tng_file *tng_file=malloc(sizeof *tng_file);
332 tng_file->pos=NULL;
333 tng_file->vel=NULL;
334 tng_file->ipos=NULL;
335 tng_file->ivel=NULL;
336 tng_file->nframes=0;
337 tng_file->chunky=chunky;
338 tng_file->precision=precision;
339 tng_file->natoms=natoms;
340 tng_file->writevel=writevel;
341 tng_file->velprecision=velprecision;
342 tng_file->initial_coding=initial_coding;
343 tng_file->initial_coding_parameter=initial_coding_parameter;
344 tng_file->coding=coding;
345 tng_file->coding_parameter=coding_parameter;
346 tng_file->initial_velcoding=initial_velcoding;
347 tng_file->initial_velcoding_parameter=initial_velcoding_parameter;
348 tng_file->velcoding=velcoding;
349 tng_file->velcoding_parameter=velcoding_parameter;
350 tng_file->speed=speed;
351 tng_file->pos=malloc(natoms*chunky*3*sizeof *tng_file->pos);
352 tng_file->f=fopen(filename,"wb");
353 if (writevel)
354 tng_file->vel=malloc(natoms*chunky*3*sizeof *tng_file->vel);
355 fwrite_int_le(&natoms,tng_file->f);
356 return tng_file;
357 }
358
359 static struct tng_file *open_tng_file_write_int(char *filename,
360 int natoms,int chunky,
361 int writevel,
362 int initial_coding,
363 int initial_coding_parameter,
364 int coding,
365 int coding_parameter,
366 int initial_velcoding,
367 int initial_velcoding_parameter,
368 int velcoding,
369 int velcoding_parameter,
370 int speed)
371 {
372 struct tng_file *tng_file=malloc(sizeof *tng_file);
373 tng_file->pos=NULL;
374 tng_file->vel=NULL;
375 tng_file->ipos=NULL;
376 tng_file->ivel=NULL;
377 tng_file->nframes=0;
378 tng_file->chunky=chunky;
379 tng_file->natoms=natoms;
380 tng_file->writevel=writevel;
381 tng_file->initial_coding=initial_coding;
382 tng_file->initial_coding_parameter=initial_coding_parameter;
383 tng_file->coding=coding;
384 tng_file->coding_parameter=coding_parameter;
385 tng_file->initial_velcoding=initial_velcoding;
386 tng_file->initial_velcoding_parameter=initial_velcoding_parameter;
387 tng_file->velcoding=velcoding;
388 tng_file->velcoding_parameter=velcoding_parameter;
389 tng_file->speed=speed;
390 tng_file->ipos=malloc(natoms*chunky*3*sizeof *tng_file->ipos);
391 tng_file->f=fopen(filename,"wb");
392 if (writevel)
393 tng_file->ivel=malloc(natoms*chunky*3*sizeof *tng_file->ivel);
394 fwrite_int_le(&natoms,tng_file->f);
395 return tng_file;
396 }
397
398 static void flush_tng_frames(struct tng_file *tng_file,
399 unsigned long prec_hi, unsigned long prec_lo,
400 unsigned long velprec_hi, unsigned long velprec_lo)
401 {
402 int algo[4];
403 char *buf;
404 int nitems;
405
406 /* Make sure these variables are used to avoid compilation warnings */
407 (void)prec_hi;
408 (void)prec_lo;
409 (void)velprec_hi;
410 (void)velprec_lo;
411
412 fwrite_int_le(&tng_file->nframes,tng_file->f);
413 algo[0]=tng_file->initial_coding;
414 algo[1]=tng_file->initial_coding_parameter;
415 algo[2]=tng_file->coding;
416 algo[3]=tng_file->coding_parameter;
417 #ifdef RECOMPRESS
418 buf=tng_compress_pos_int(tng_file->ipos,
419 tng_file->natoms,
420 tng_file->nframes,
421 prec_hi,prec_lo,
422 tng_file->speed,algo,&nitems);
423 #else /* RECOMPRESS */
424 #ifdef TEST_FLOAT
425 buf=tng_compress_pos_float(tng_file->pos,
426 tng_file->natoms,
427 tng_file->nframes,
428 tng_file->precision,
429 tng_file->speed,algo,&nitems);
430 #else /* TEST_FLOAT */
431 buf=tng_compress_pos(tng_file->pos,
432 tng_file->natoms,
433 tng_file->nframes,
434 tng_file->precision,
435 tng_file->speed,algo,&nitems);
436 #endif /* TEST_FLOAT */
437 #endif /* RECOMPRESS */
438 tng_file->initial_coding=algo[0];
439 tng_file->initial_coding_parameter=algo[1];
440 tng_file->coding=algo[2];
441 tng_file->coding_parameter=algo[3];
442 fwrite_int_le(&nitems,tng_file->f);
443 fwrite(buf,1,nitems,tng_file->f);
444 free(buf);
445 if (tng_file->writevel)
446 {
447 algo[0]=tng_file->initial_velcoding;
448 algo[1]=tng_file->initial_velcoding_parameter;
449 algo[2]=tng_file->velcoding;
450 algo[3]=tng_file->velcoding_parameter;
451 #ifdef RECOMPRESS
452 buf=tng_compress_vel_int(tng_file->ivel,
453 tng_file->natoms,
454 tng_file->nframes,
455 velprec_hi,velprec_lo,
456 tng_file->speed,algo,&nitems);
457 #else /* RECOMPRESS */
458 #ifdef TEST_FLOAT
459 buf=tng_compress_vel_float(tng_file->vel,
460 tng_file->natoms,
461 tng_file->nframes,
462 tng_file->velprecision,
463 tng_file->speed,algo,&nitems);
464 #else /* TEST_FLOAT */
465 buf=tng_compress_vel(tng_file->vel,
466 tng_file->natoms,
467 tng_file->nframes,
468 tng_file->velprecision,
469 tng_file->speed,algo,&nitems);
470 #endif /* TEST_FLOAT */
471 #endif /* RECOMPRESS */
472 tng_file->initial_velcoding=algo[0];
473 tng_file->initial_velcoding_parameter=algo[1];
474 tng_file->velcoding=algo[2];
475 tng_file->velcoding_parameter=algo[3];
476 fwrite_int_le(&nitems,tng_file->f);
477 fwrite(buf,1,nitems,tng_file->f);
478 free(buf);
479 }
480 tng_file->nframes=0;
481 }
482
483 static void write_tng_file(struct tng_file *tng_file,
484 REAL *pos,REAL *vel)
485 {
486 memcpy(tng_file->pos+tng_file->nframes*tng_file->natoms*3,pos,tng_file->natoms*3*sizeof *tng_file->pos);
487 if (tng_file->writevel)
488 memcpy(tng_file->vel+tng_file->nframes*tng_file->natoms*3,vel,tng_file->natoms*3*sizeof *tng_file->vel);
489 tng_file->nframes++;
490 if (tng_file->nframes==tng_file->chunky)
491 flush_tng_frames(tng_file,0,0,0,0);
492 }
493
494 static void write_tng_file_int(struct tng_file *tng_file,
495 int *ipos,int *ivel,
496 unsigned long prec_hi, unsigned long prec_lo,
497 unsigned long velprec_hi, unsigned long velprec_lo)
498 {
499 memcpy(tng_file->ipos+tng_file->nframes*tng_file->natoms*3,ipos,tng_file->natoms*3*sizeof *tng_file->ipos);
500 if (tng_file->writevel)
501 memcpy(tng_file->ivel+tng_file->nframes*tng_file->natoms*3,ivel,tng_file->natoms*3*sizeof *tng_file->ivel);
502 tng_file->nframes++;
503 if (tng_file->nframes==tng_file->chunky)
504 flush_tng_frames(tng_file,prec_hi,prec_lo,velprec_hi,velprec_lo);
505 tng_file->prec_hi=prec_hi;
506 tng_file->prec_lo=prec_lo;
507 tng_file->velprec_hi=velprec_hi;
508 tng_file->velprec_lo=velprec_lo;
509 }
510
511 static void close_tng_file_write(struct tng_file *tng_file)
512 {
513 if (tng_file->nframes)
514 flush_tng_frames(tng_file,tng_file->prec_hi,tng_file->prec_lo,tng_file->velprec_hi,tng_file->velprec_lo);
515 fclose(tng_file->f);
516 free(tng_file->pos);
517 free(tng_file->vel);
518 free(tng_file->ipos);
519 free(tng_file->ivel);
520 free(tng_file);
521 }
522
523 static struct tng_file *open_tng_file_read(char *filename, int writevel)
524 {
525 struct tng_file *tng_file=malloc(sizeof *tng_file);
526 tng_file->pos=NULL;
527 tng_file->vel=NULL;
528 tng_file->ipos=NULL;
529 tng_file->ivel=NULL;
530 tng_file->f=fopen(filename,"rb");
531 tng_file->nframes=0;
532 tng_file->nframes_delivered=0;
533 tng_file->writevel=writevel;
534 if (tng_file->f)
535 fread_int_le(&tng_file->natoms,tng_file->f);
536 else
537 {
538 free(tng_file);
539 tng_file=NULL;
540 }
541 return tng_file;
542 }
543
544 static struct tng_file *open_tng_file_read_int(char *filename, int writevel)
545 {
546 struct tng_file *tng_file=malloc(sizeof *tng_file);
547 tng_file->pos=NULL;
548 tng_file->vel=NULL;
549 tng_file->ipos=NULL;
550 tng_file->ivel=NULL;
551 tng_file->f=fopen(filename,"rb");
552 tng_file->nframes=0;
553 tng_file->nframes_delivered=0;
554 tng_file->writevel=writevel;
555 if (tng_file->f)
556 fread_int_le(&tng_file->natoms,tng_file->f);
557 else
558 {
559 free(tng_file);
560 tng_file=NULL;
561 }
562 return tng_file;
563 }
564
565 static int read_tng_file(struct tng_file *tng_file,
566 REAL *pos,
567 REAL *vel)
568 {
569 if (tng_file->nframes==tng_file->nframes_delivered)
570 {
571 int nitems;
572 char *buf;
573 free(tng_file->pos);
574 free(tng_file->vel);
575 if (!fread_int_le(&tng_file->nframes,tng_file->f))
576 return 1;
577 if (!fread_int_le(&nitems,tng_file->f))
578 return 1;
579 buf=malloc(nitems);
580 if (!fread(buf,1,nitems,tng_file->f))
581 {
582 free(buf);
583 return 1;
584 }
585 tng_file->pos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->pos);
586 if (tng_file->writevel)
587 tng_file->vel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->vel);
588 #if 0
589 {
590 int natoms, nframes, algo[4];
591 double precision;
592 int ivel;
593 char *initial_coding, *coding;
594 tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo);
595 initial_coding=tng_compress_initial_pos_algo(algo);
596 coding=tng_compress_pos_algo(algo);
597 printf("ivel=%d natoms=%d nframes=%d precision=%g initial pos=%s pos=%s\n",ivel,natoms,nframes,precision,initial_coding,coding);
598 }
599 #endif
600 #ifdef TEST_FLOAT
601 tng_compress_uncompress_float(buf,tng_file->pos);
602 #else
603 tng_compress_uncompress(buf,tng_file->pos);
604 #endif
605 free(buf);
606 if (tng_file->writevel)
607 {
608 if (!fread_int_le(&nitems,tng_file->f))
609 return 1;
610 buf=malloc(nitems);
611 if (!fread(buf,1,nitems,tng_file->f))
612 {
613 free(buf);
614 return 1;
615 }
616 #if 0
617 {
618 int natoms, nframes, algo[4];
619 double precision;
620 int ivel;
621 char *initial_coding, *coding;
622 tng_compress_inquire(buf,&ivel,&natoms,&nframes,&precision,algo);
623 initial_coding=tng_compress_initial_vel_algo(algo);
624 coding=tng_compress_vel_algo(algo);
625 printf("ivel=%d natoms=%d nframes=%d precision=%g initial vel=%s vel=%s\n",ivel,natoms,nframes,precision,initial_coding,coding);
626 }
627 #endif
628 #ifdef TEST_FLOAT
629 tng_compress_uncompress_float(buf,tng_file->vel);
630 #else
631 tng_compress_uncompress(buf,tng_file->vel);
632 #endif
633 free(buf);
634 }
635 tng_file->nframes_delivered=0;
636 }
637 memcpy(pos,tng_file->pos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *pos);
638 if (tng_file->writevel)
639 memcpy(vel,tng_file->vel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *vel);
640 tng_file->nframes_delivered++;
641 return 0;
642 }
643
644 static int read_tng_file_int(struct tng_file *tng_file,
645 int *ipos,
646 int *ivel,
647 unsigned long *prec_hi, unsigned long *prec_lo,
648 unsigned long *velprec_hi, unsigned long *velprec_lo)
649 {
650 if (tng_file->nframes==tng_file->nframes_delivered)
651 {
652 int nitems;
653 char *buf;
654 free(tng_file->ipos);
655 free(tng_file->ivel);
656 if (!fread_int_le(&tng_file->nframes,tng_file->f))
657 return 1;
658 if (!fread_int_le(&nitems,tng_file->f))
659 return 1;
660 buf=malloc(nitems);
661 if (!fread(buf,1,nitems,tng_file->f))
662 {
663 free(buf);
664 return 1;
665 }
666 tng_file->ipos=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->ipos);
667 if (tng_file->writevel)
668 tng_file->ivel=malloc(tng_file->natoms*tng_file->nframes*3*sizeof *tng_file->ivel);
669 tng_compress_uncompress_int(buf,tng_file->ipos,prec_hi,prec_lo);
670 free(buf);
671 if (tng_file->writevel)
672 {
673 if (!fread_int_le(&nitems,tng_file->f))
674 return 1;
675 buf=malloc(nitems);
676 if (!fread(buf,1,nitems,tng_file->f))
677 {
678 free(buf);
679 return 1;
680 }
681 tng_compress_uncompress_int(buf,tng_file->ivel,velprec_hi,velprec_lo);
682 free(buf);
683 }
684 tng_file->nframes_delivered=0;
685 }
686 memcpy(ipos,tng_file->ipos+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *ipos);
687 if (tng_file->writevel)
688 memcpy(ivel,tng_file->ivel+tng_file->nframes_delivered*tng_file->natoms*3,tng_file->natoms*3*sizeof *ivel);
689 tng_file->nframes_delivered++;
690 return 0;
691 }
692
693 static void close_tng_file_read(struct tng_file *tng_file)
694 {
695 free(tng_file->vel);
696 free(tng_file->pos);
697 free(tng_file->ivel);
698 free(tng_file->ipos);
699 fclose(tng_file->f);
700 free(tng_file);
701 }
702
703
704
705 #ifndef EXPECTED_FILESIZE
706 #define EXPECTED_FILESIZE 1
707 #endif
708
709 #ifndef INITIALVELCODING
710 #define INITIALVELCODING -1
711 #endif
712 #ifndef INITIALVELCODINGPARAMETER
713 #define INITIALVELCODINGPARAMETER -1
714 #endif
715
716 #ifndef SPEED
717 #define SPEED 5
718 #endif
719
720 /* Return value 1 means file error.
721 Return value 4 means coding error in coordinates.
722 Return value 5 means coding error in velocities.
723 Return value 9 means filesize seems too off.
724
725 Return value 100+ means test specific error.
726 */
727 static int algotest()
728 {
729 int i;
730 int *intbox=malloc(NATOMS*3*sizeof *intbox);
731 int *intvelbox=malloc(NATOMS*3*sizeof *intvelbox);
732 #ifdef RECOMPRESS
733 unsigned long pos_prec_hi,pos_prec_lo;
734 unsigned long vel_prec_hi,vel_prec_lo;
735 #endif
736 REAL *box1=malloc(NATOMS*STRIDE1*sizeof *box1);
737 REAL *velbox1=malloc(NATOMS*STRIDE1*sizeof *velbox1);
738 int startframe=0;
739 int endframe=NFRAMES;
740 #ifdef GEN
741 FILE *file;
742 REAL filesize;
743 REAL *box2=0;
744 REAL *velbox2=0;
745 #else
746 int readreturn;
747 REAL *box2=malloc(NATOMS*STRIDE2*sizeof *box2);
748 REAL *velbox2=malloc(NATOMS*STRIDE2*sizeof *velbox2);
749 #endif
750 #ifdef RECOMPRESS
751 void *dumpfile=open_tng_file_write_int(TNG_COMPRESS_FILES_DIR FILENAME,NATOMS,CHUNKY,
752 WRITEVEL,
753 INITIALCODING,
754 INITIALCODINGPARAMETER,CODING,CODINGPARAMETER,
755 INITIALVELCODING,INITIALVELCODINGPARAMETER,
756 VELCODING,VELCODINGPARAMETER,SPEED);
757 void *dumpfile_recompress=open_tng_file_read_int(TNG_COMPRESS_FILES_DIR RECOMPRESS,WRITEVEL);
758 if (!dumpfile_recompress)
759 {
760 free(intbox);
761 free(intvelbox);
762 free(box1);
763 free(velbox1);
764 if(box2)
765 free(box2);
766 if(velbox2)
767 free(velbox2);
768 return 1;
769 }
770 #else /* RECOMPRESS */
771 #ifdef GEN
772 void *dumpfile=open_tng_file_write(TNG_COMPRESS_FILES_DIR FILENAME,NATOMS,CHUNKY,
773 PRECISION,WRITEVEL,VELPRECISION,
774 INITIALCODING,
775 INITIALCODINGPARAMETER,CODING,CODINGPARAMETER,
776 INITIALVELCODING,INITIALVELCODINGPARAMETER,
777 VELCODING,VELCODINGPARAMETER,SPEED);
778 #else
779 void *dumpfile=open_tng_file_read(TNG_COMPRESS_FILES_DIR FILENAME,WRITEVEL);
780 #endif
781 #endif /* RECOMPRESS */
782 if (!dumpfile)
783 {
784 free(intbox);
785 free(intvelbox);
786 free(box1);
787 free(velbox1);
788 if(box2)
789 free(box2);
790 if(velbox2)
791 free(velbox2);
792 return 1;
793 }
794 for (i=startframe; i<endframe; i++)
795 {
796 #ifdef RECOMPRESS
797 unsigned long prec_hi, prec_lo;
798 unsigned long velprec_hi, velprec_lo;
799 if (read_tng_file_int(dumpfile_recompress,intbox,intvelbox,&prec_hi,&prec_lo,&velprec_hi,&velprec_lo))
800 return 1;
801 write_tng_file_int(dumpfile,intbox,intvelbox,prec_hi,prec_lo,velprec_hi,velprec_lo);
802 #else /* RECOMPRESS */
803 genibox(intbox,i);
804 realbox(intbox,box1,STRIDE1);
805 #if WRITEVEL
806 genivelbox(intvelbox,i);
807 realvelbox(intvelbox,velbox1,STRIDE1);
808 #endif
809 #ifdef GEN
810 write_tng_file(dumpfile,box1,velbox1);
811 #else /* GEN */
812 #ifdef INTTOFLOAT
813 {
814 unsigned long prec_hi, prec_lo;
815 unsigned long velprec_hi, velprec_lo;
816 readreturn=read_tng_file_int(dumpfile,intbox,intvelbox,&prec_hi,&prec_lo,&velprec_hi,&velprec_lo);
817 if (!readreturn)
818 {
819 tng_compress_int_to_float(intbox,prec_hi,prec_lo,NATOMS,1,box2);
820 #if WRITEVEL
821 tng_compress_int_to_float(intvelbox,velprec_hi,velprec_lo,NATOMS,1,velbox2);
822 #endif
823 }
824 }
825 #else /* INTTOFLOAT */
826 #ifdef INTTODOUBLE
827 {
828 unsigned long prec_hi, prec_lo;
829 unsigned long velprec_hi, velprec_lo;
830 readreturn=read_tng_file_int(dumpfile,intbox,intvelbox,&prec_hi,&prec_lo,&velprec_hi,&velprec_lo);
831 if (!readreturn)
832 {
833 tng_compress_int_to_double(intbox,prec_hi,prec_lo,NATOMS,1,box2);
834 #if WRITEVEL
835 tng_compress_int_to_double(intvelbox,velprec_hi,velprec_lo,NATOMS,1,velbox2);
836 #endif
837 }
838 }
839 #else /* INTTODOUBLE */
840 readreturn=read_tng_file(dumpfile,box2,velbox2);
841 #endif /* INTTODOUBLE */
842 #endif /* INTTOFLOAT */
843 if (readreturn==1) /* general read error */
844 {
845 free(intbox);
846 free(intvelbox);
847 free(box1);
848 free(velbox1);
849 if(box2)
850 free(box2);
851 if(velbox2)
852 free(velbox2);
853 return 1;
854 }
855 #endif /* GEN */
856 #ifndef GEN
857 /* Check for equality of boxes. */
858 if (!equalarr(box1,box2,(REAL)PRECISION,NATOMS,3,STRIDE1,STRIDE2))
859 {
860 free(intbox);
861 free(intvelbox);
862 free(box1);
863 free(velbox1);
864 if(box2)
865 free(box2);
866 if(velbox2)
867 free(velbox2);
868 return 4;
869 }
870 #if WRITEVEL
871 if (!equalarr(velbox1,velbox2,(REAL)VELPRECISION,NATOMS,3,STRIDE1,STRIDE2))
872 {
873 free(intbox);
874 free(intvelbox);
875 free(box1);
876 free(velbox1);
877 if(box2)
878 free(box2);
879 if(velbox2)
880 free(velbox2);
881 return 5;
882 }
883 #endif
884 #endif /* GEN */
885 #endif /* RECOMPRESS */
886 }
887 #ifdef GEN
888 close_tng_file_write(dumpfile);
889 #else
890 close_tng_file_read(dumpfile);
891 #endif
892 #ifdef RECOMPRESS
893 close_tng_file_read(dumpfile_recompress);
894 #endif
895 #ifdef GEN
896 /* Check against expected filesize for this test. */
897 if (!(file=fopen(TNG_COMPRESS_FILES_DIR FILENAME,"rb")))
898 {
899 fprintf(stderr,"ERROR: Cannot open file "TNG_COMPRESS_FILES_DIR FILENAME"\n");
900 exit(EXIT_FAILURE);
901 }
902 filesize=0;
903 while(1)
904 {
905 char b;
906 if (!fread(&b,1,1,file))
907 break;
908 filesize++;
909 }
910 fclose(file);
911 if (filesize>0)
912 {
913 if ((fabs(filesize-EXPECTED_FILESIZE)/EXPECTED_FILESIZE)>0.05)
914 {
915 free(intvelbox);
916 free(box1);
917 free(velbox1);
918 if(box2)
919 free(box2);
920 if(velbox2)
921 free(velbox2);
922 return 9;
923 }
924 }
925 #endif
926 free(intvelbox);
927 free(box1);
928 free(velbox1);
929 if(box2)
930 free(box2);
931 if(velbox2)
932 free(velbox2);
933 return 0;
934 }
935
936 int main()
937 {
938 int testval;
939 if (sizeof(int)<4)
940 {
941 fprintf(stderr,"ERROR: sizeof(int) is too small: %d<4\n",(int)sizeof(int));
942 exit(EXIT_FAILURE);
943 }
944 #ifdef GEN
945 printf("Tng compress testsuite generating (writing) test: %s\n",TESTNAME);
946 #else
947 printf("Tng compress testsuite running (reading) test: %s\n",TESTNAME);
948 #endif
949 testval=algotest();
950 if (testval==0)
951 printf("Passed.\n");
952 else if (testval==1)
953 {
954 printf("ERROR: File error.\n");
955 exit(EXIT_FAILURE);
956 }
957 else if (testval==4)
958 {
959 printf("ERROR: Read coding error in coordinates.\n");
960 exit(EXIT_FAILURE);
961 }
962 else if (testval==5)
963 {
964 printf("ERROR: Read coding error in velocities.\n");
965 exit(EXIT_FAILURE);
966 }
967 else if (testval==9)
968 {
969 printf("ERROR: Generated filesize differs too much.\n");
970 exit(EXIT_FAILURE);
971 }
972 else
973 {
974 printf("ERROR: Unknown error.\n");
975 exit(EXIT_FAILURE);
976 }
977 return 0;
978 }
The ultimate molecular dynamics simulation package