| File: | external/tng_io/src/compression/xtc3.c |
| Location: | line 761, column 35 |
| Description: | The left expression of the compound assignment is an uninitialized value. The computed value will also be garbage |
| 1 | /* This code is part of the tng compression routines. | |||
| 2 | * | |||
| 3 | * Written by Daniel Spangberg | |||
| 4 | * Copyright (c) 2010, 2013, The GROMACS development team. | |||
| 5 | * | |||
| 6 | * | |||
| 7 | * This program is free software; you can redistribute it and/or | |||
| 8 | * modify it under the terms of the Revised BSD License. | |||
| 9 | */ | |||
| 10 | ||||
| 11 | /* This code is heavily influenced by | |||
| 12 | http://hpcv100.rc.rug.nl/xdrf.html | |||
| 13 | Based on coordinate compression (c) by Frans van Hoesel. | |||
| 14 | and GROMACS xtc files (http://www.gromacs.org) | |||
| 15 | (c) Copyright (c) Erik Lindahl, David van der Spoel | |||
| 16 | */ | |||
| 17 | ||||
| 18 | /* The cost estimates are ripped right out of xtc2.c, so take these | |||
| 19 | with a grain (truckload) of salt. */ | |||
| 20 | ||||
| 21 | #include <limits.h> | |||
| 22 | #include <stdio.h> | |||
| 23 | #include <stdlib.h> | |||
| 24 | #include <string.h> | |||
| 25 | #include <math.h> | |||
| 26 | #include "../../include/compression/warnmalloc.h" | |||
| 27 | #include "../../include/compression/widemuldiv.h" | |||
| 28 | #include "../../include/compression/bwlzh.h" | |||
| 29 | ||||
| 30 | static const double iflipgaincheck=0.89089871814033927; /* 1./(2**(1./6)) */ | |||
| 31 | ||||
| 32 | #define MAX_LARGE_RLE1024 1024 /* Maximum number of large atoms for large RLE. */ | |||
| 33 | #define MAX_SMALL_RLE12 12 /* Maximum number of small atoms in one group. */ | |||
| 34 | ||||
| 35 | #define TRESHOLD_INTRA_INTER_DIRECT1.5 1.5 /* How much larger can the direct | |||
| 36 | frame deltas for the small | |||
| 37 | triplets be and be accepted anyway | |||
| 38 | as better than the intra/inter frame | |||
| 39 | deltas. For better instructions/RLEs. */ | |||
| 40 | ||||
| 41 | #define TRESHOLD_INTER_INTRA5.0 5.0 /* How much larger can the intra | |||
| 42 | frame deltas for the small | |||
| 43 | triplets be and be accepted anyway | |||
| 44 | as better than the inter frame | |||
| 45 | deltas. */ | |||
| 46 | ||||
| 47 | /* Difference in indices used for determining whether to store as | |||
| 48 | large or small. A fun detail in this compression algorithm is that | |||
| 49 | if everything works fine, large can often be smaller than small, or | |||
| 50 | at least not as large as is large in magic.c. This is a key idea of | |||
| 51 | xtc3. */ | |||
| 52 | #define QUITE_LARGE3 3 | |||
| 53 | #define IS_LARGE6 6 | |||
| 54 | ||||
| 55 | #if 0 | |||
| 56 | #define SHOWIT | |||
| 57 | #endif | |||
| 58 | ||||
| 59 | #if 0 | |||
| 60 | #define SHOWIT_LIGHT | |||
| 61 | #endif | |||
| 62 | ||||
| 63 | /* These routines are in xtc2.c */ | |||
| 64 | int Ptngc_magic(unsigned int i); | |||
| 65 | int Ptngc_find_magic_index(unsigned int maxval); | |||
| 66 | ||||
| 67 | static unsigned int positive_int(int item) | |||
| 68 | { | |||
| 69 | int s=0; | |||
| 70 | if (item>0) | |||
| 71 | s=1+(item-1)*2; | |||
| 72 | else if (item<0) | |||
| 73 | s=2+(-item-1)*2; | |||
| 74 | return s; | |||
| 75 | } | |||
| 76 | ||||
| 77 | static int unpositive_int(int val) | |||
| 78 | { | |||
| 79 | int s=(val+1)/2; | |||
| 80 | if ((val%2)==0) | |||
| 81 | s=-s; | |||
| 82 | return s; | |||
| 83 | } | |||
| 84 | ||||
| 85 | ||||
| 86 | /* Sequence instructions */ | |||
| 87 | #define INSTR_DEFAULT0U 0U | |||
| 88 | #define INSTR_SMALL_RUNLENGTH1U 1U | |||
| 89 | #define INSTR_ONLY_LARGE2U 2U | |||
| 90 | #define INSTR_ONLY_SMALL3U 3U | |||
| 91 | #define INSTR_FLIP4U 4U | |||
| 92 | #define INSTR_LARGE_RLE5U 5U | |||
| 93 | #define INSTR_LARGE_DIRECT6U 6U | |||
| 94 | #define INSTR_LARGE_INTRA_DELTA7U 7U | |||
| 95 | #define INSTR_LARGE_INTER_DELTA8U 8U | |||
| 96 | ||||
| 97 | #define MAXINSTR9 9 | |||
| 98 | ||||
| 99 | struct xtc3_context | |||
| 100 | { | |||
| 101 | unsigned int *instructions; | |||
| 102 | int ninstr, ninstr_alloc; | |||
| 103 | unsigned int *rle; | |||
| 104 | int nrle, nrle_alloc; | |||
| 105 | unsigned int *large_direct; | |||
| 106 | int nlargedir, nlargedir_alloc; | |||
| 107 | unsigned int *large_intra_delta; | |||
| 108 | int nlargeintra, nlargeintra_alloc; | |||
| 109 | unsigned int *large_inter_delta; | |||
| 110 | int nlargeinter, nlargeinter_alloc; | |||
| 111 | unsigned int *smallintra; | |||
| 112 | int nsmallintra, nsmallintra_alloc; | |||
| 113 | int minint[3],maxint[3]; | |||
| 114 | int has_large; | |||
| 115 | int has_large_ints[MAX_LARGE_RLE1024*3]; /* Large cache. */ | |||
| 116 | int has_large_type[MAX_LARGE_RLE1024]; /* What kind of type this large | |||
| 117 | int is. */ | |||
| 118 | int current_large_type; | |||
| 119 | }; | |||
| 120 | ||||
| 121 | static void init_xtc3_context(struct xtc3_context *xtc3_context) | |||
| 122 | { | |||
| 123 | xtc3_context->instructions=NULL((void*)0); | |||
| 124 | xtc3_context->ninstr=0; | |||
| 125 | xtc3_context->ninstr_alloc=0; | |||
| 126 | xtc3_context->rle=NULL((void*)0); | |||
| 127 | xtc3_context->nrle=0; | |||
| 128 | xtc3_context->nrle_alloc=0; | |||
| 129 | xtc3_context->large_direct=NULL((void*)0); | |||
| 130 | xtc3_context->nlargedir=0; | |||
| 131 | xtc3_context->nlargedir_alloc=0; | |||
| 132 | xtc3_context->large_intra_delta=NULL((void*)0); | |||
| 133 | xtc3_context->nlargeintra=0; | |||
| 134 | xtc3_context->nlargeintra_alloc=0; | |||
| 135 | xtc3_context->large_inter_delta=NULL((void*)0); | |||
| 136 | xtc3_context->nlargeinter=0; | |||
| 137 | xtc3_context->nlargeinter_alloc=0; | |||
| 138 | xtc3_context->smallintra=NULL((void*)0); | |||
| 139 | xtc3_context->nsmallintra=0; | |||
| 140 | xtc3_context->nsmallintra_alloc=0; | |||
| 141 | xtc3_context->has_large=0; | |||
| 142 | xtc3_context->current_large_type=0; | |||
| 143 | } | |||
| 144 | ||||
| 145 | static void free_xtc3_context(struct xtc3_context *xtc3_context) | |||
| 146 | { | |||
| 147 | free(xtc3_context->instructions); | |||
| 148 | free(xtc3_context->rle); | |||
| 149 | free(xtc3_context->large_direct); | |||
| 150 | free(xtc3_context->large_intra_delta); | |||
| 151 | free(xtc3_context->large_inter_delta); | |||
| 152 | free(xtc3_context->smallintra); | |||
| 153 | } | |||
| 154 | ||||
| 155 | /* Modifies three integer values for better compression of water */ | |||
| 156 | static void swap_ints(int *in, int *out) | |||
| 157 | { | |||
| 158 | out[0]=in[0]+in[1]; | |||
| 159 | out[1]=-in[1]; | |||
| 160 | out[2]=in[1]+in[2]; | |||
| 161 | } | |||
| 162 | ||||
| 163 | static void swap_is_better(int *input, int *minint, int *sum_normal, int *sum_swapped) | |||
| 164 | { | |||
| 165 | int normal_max=0; | |||
| 166 | int swapped_max=0; | |||
| 167 | int i,j; | |||
| 168 | int normal[3]; | |||
| 169 | int swapped[3]; | |||
| 170 | for (i=0; i<3; i++) | |||
| 171 | { | |||
| 172 | normal[0]=input[i]-minint[i]; | |||
| 173 | normal[1]=input[3+i]-input[i]; /* minint[i]-minint[i] cancels out */ | |||
| 174 | normal[2]=input[6+i]-input[3+i]; /* minint[i]-minint[i] cancels out */ | |||
| 175 | swap_ints(normal,swapped); | |||
| 176 | for (j=1; j<3; j++) | |||
| 177 | { | |||
| 178 | if (positive_int(normal[j])>(unsigned int)normal_max) | |||
| 179 | normal_max=positive_int(normal[j]); | |||
| 180 | if (positive_int(swapped[j])>(unsigned int)swapped_max) | |||
| 181 | swapped_max=positive_int(swapped[j]); | |||
| 182 | } | |||
| 183 | } | |||
| 184 | if (normal_max==0) | |||
| 185 | normal_max=1; | |||
| 186 | if (swapped_max==0) | |||
| 187 | swapped_max=1; | |||
| 188 | *sum_normal=normal_max; | |||
| 189 | *sum_swapped=swapped_max; | |||
| 190 | } | |||
| 191 | ||||
| 192 | static void allocate_enough_memory(unsigned int **ptr, int *nele, int *nele_alloc) | |||
| 193 | { | |||
| 194 | (*nele)++; | |||
| 195 | if (*nele>*nele_alloc) | |||
| 196 | { | |||
| 197 | *nele_alloc=*nele + *nele/2; | |||
| 198 | *ptr=warnrealloc(*ptr,*nele_alloc*sizeof **ptr)Ptngc_warnrealloc_x(*ptr,*nele_alloc*sizeof **ptr,"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,198); | |||
| 199 | } | |||
| 200 | } | |||
| 201 | ||||
| 202 | static void insert_value_in_array(unsigned int **ptr, int *nele, int *nele_alloc, | |||
| 203 | unsigned int value, | |||
| 204 | char *arrayname) | |||
| 205 | { | |||
| 206 | #ifndef SHOWIT | |||
| 207 | (void)arrayname; | |||
| 208 | #endif | |||
| 209 | allocate_enough_memory(ptr,nele,nele_alloc); | |||
| 210 | #ifdef SHOWIT | |||
| 211 | fprintf(stderrstderr,"Inserting value %u into array %s @ %d\n",value,arrayname,(*nele)-1); | |||
| 212 | #endif | |||
| 213 | (*ptr)[(*nele)-1]=value; | |||
| 214 | } | |||
| 215 | ||||
| 216 | ||||
| 217 | ||||
| 218 | static void swapdecide(struct xtc3_context *xtc3_context, int *input,int *swapatoms, int *large_index, int *minint) | |||
| 219 | { | |||
| 220 | int didswap=0; | |||
| 221 | int normal,swapped; | |||
| 222 | (void)large_index; | |||
| 223 | swap_is_better(input,minint,&normal,&swapped); | |||
| 224 | /* We have to determine if it is worth to change the behaviour. | |||
| 225 | If diff is positive it means that it is worth something to | |||
| 226 | swap. But it costs 4 bits to do the change. If we assume that | |||
| 227 | we gain 0.17 bit by the swap per value, and the runlength>2 | |||
| 228 | for four molecules in a row, we gain something. So check if we | |||
| 229 | gain at least 0.17 bits to even attempt the swap. | |||
| 230 | */ | |||
| 231 | #ifdef SHOWIT | |||
| 232 | fprintf(stderrstderr,"Trying Flip: %g %g\n",(double)swapped/normal, (double)normal/swapped); | |||
| 233 | #endif | |||
| 234 | if (((swapped<normal) && (fabs((double)swapped/normal)<iflipgaincheck)) || | |||
| 235 | ((normal<swapped) && (fabs((double)normal/swapped)<iflipgaincheck))) | |||
| 236 | { | |||
| 237 | if (swapped<normal) | |||
| 238 | { | |||
| 239 | if (!*swapatoms) | |||
| 240 | { | |||
| 241 | *swapatoms=1; | |||
| 242 | didswap=1; | |||
| 243 | } | |||
| 244 | } | |||
| 245 | else | |||
| 246 | { | |||
| 247 | if (*swapatoms) | |||
| 248 | { | |||
| 249 | *swapatoms=0; | |||
| 250 | didswap=1; | |||
| 251 | } | |||
| 252 | } | |||
| 253 | } | |||
| 254 | if (didswap) | |||
| 255 | { | |||
| 256 | #ifdef SHOWIT | |||
| 257 | fprintf(stderrstderr,"Flip. Swapatoms is now %d\n",*swapatoms); | |||
| 258 | #endif | |||
| 259 | insert_value_in_array(&xtc3_context->instructions, | |||
| 260 | &xtc3_context->ninstr, | |||
| 261 | &xtc3_context->ninstr_alloc, | |||
| 262 | INSTR_FLIP4U,"instr"); | |||
| 263 | } | |||
| 264 | } | |||
| 265 | ||||
| 266 | /* It is "large" if we have to increase the small index quite a | |||
| 267 | bit. Not so much to be rejected by the not very large check | |||
| 268 | later. */ | |||
| 269 | static int is_quite_large(int *input, int small_index, int max_large_index) | |||
| 270 | { | |||
| 271 | int is=0; | |||
| 272 | int i; | |||
| 273 | if (small_index+QUITE_LARGE3>=max_large_index) | |||
| 274 | is=1; | |||
| 275 | else | |||
| 276 | { | |||
| 277 | for (i=0; i<3; i++) | |||
| 278 | if (positive_int(input[i])>(unsigned int)Ptngc_magic(small_index+QUITE_LARGE3)) | |||
| 279 | { | |||
| 280 | is=1; | |||
| 281 | break; | |||
| 282 | } | |||
| 283 | } | |||
| 284 | return is; | |||
| 285 | } | |||
| 286 | ||||
| 287 | #ifdef SHOWIT | |||
| 288 | int nbits_sum; | |||
| 289 | int nvalues_sum; | |||
| 290 | #endif | |||
| 291 | ||||
| 292 | static void insert_batch(int *input_ptr, int ntriplets_left, int *prevcoord, int *encode_ints, int startenc, int *nenc) | |||
| 293 | { | |||
| 294 | int nencode=startenc*3; | |||
| 295 | int tmp_prevcoord[3]; | |||
| 296 | ||||
| 297 | tmp_prevcoord[0]=prevcoord[0]; | |||
| 298 | tmp_prevcoord[1]=prevcoord[1]; | |||
| 299 | tmp_prevcoord[2]=prevcoord[2]; | |||
| 300 | ||||
| 301 | if (startenc) | |||
| 302 | { | |||
| 303 | int i; | |||
| 304 | for (i=0; i<startenc; i++) | |||
| 305 | { | |||
| 306 | tmp_prevcoord[0]+=encode_ints[i*3]; | |||
| 307 | tmp_prevcoord[1]+=encode_ints[i*3+1]; | |||
| 308 | tmp_prevcoord[2]+=encode_ints[i*3+2]; | |||
| 309 | #ifdef SHOWIT | |||
| 310 | fprintf(stderrstderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",i*3, | |||
| 311 | tmp_prevcoord[0],tmp_prevcoord[1],tmp_prevcoord[2], | |||
| 312 | encode_ints[i*3], | |||
| 313 | encode_ints[i*3+1], | |||
| 314 | encode_ints[i*3+2], | |||
| 315 | positive_int(encode_ints[i*3]), | |||
| 316 | positive_int(encode_ints[i*3+1]), | |||
| 317 | positive_int(encode_ints[i*3+2])); | |||
| 318 | #endif | |||
| 319 | } | |||
| 320 | } | |||
| 321 | ||||
| 322 | #ifdef SHOWIT | |||
| 323 | fprintf(stderrstderr,"New batch\n"); | |||
| 324 | #endif | |||
| 325 | while ((nencode<3+MAX_SMALL_RLE12*3) && (nencode<ntriplets_left*3)) | |||
| 326 | { | |||
| 327 | encode_ints[nencode]=input_ptr[nencode]-tmp_prevcoord[0]; | |||
| 328 | encode_ints[nencode+1]=input_ptr[nencode+1]-tmp_prevcoord[1]; | |||
| 329 | encode_ints[nencode+2]=input_ptr[nencode+2]-tmp_prevcoord[2]; | |||
| 330 | #ifdef SHOWIT | |||
| 331 | fprintf(stderrstderr,"%6d: %6d %6d %6d\t\t%6d %6d %6d\t\t%6d %6d %6d\n",nencode, | |||
| 332 | input_ptr[nencode], | |||
| 333 | input_ptr[nencode+1], | |||
| 334 | input_ptr[nencode+2], | |||
| 335 | encode_ints[nencode], | |||
| 336 | encode_ints[nencode+1], | |||
| 337 | encode_ints[nencode+2], | |||
| 338 | positive_int(encode_ints[nencode]), | |||
| 339 | positive_int(encode_ints[nencode+1]), | |||
| 340 | positive_int(encode_ints[nencode+2])); | |||
| 341 | #endif | |||
| 342 | tmp_prevcoord[0]=input_ptr[nencode]; | |||
| 343 | tmp_prevcoord[1]=input_ptr[nencode+1]; | |||
| 344 | tmp_prevcoord[2]=input_ptr[nencode+2]; | |||
| 345 | nencode+=3; | |||
| 346 | } | |||
| 347 | *nenc=nencode; | |||
| 348 | } | |||
| 349 | ||||
| 350 | static void large_instruction_change(struct xtc3_context *xtc3_context, int i) | |||
| 351 | { | |||
| 352 | /* If the first large is of a different kind than the currently used we must | |||
| 353 | emit an "instruction" to change the large type. */ | |||
| 354 | if (xtc3_context->has_large_type[i]!=xtc3_context->current_large_type) | |||
| 355 | { | |||
| 356 | unsigned int instr; | |||
| 357 | xtc3_context->current_large_type=xtc3_context->has_large_type[i]; | |||
| 358 | if (xtc3_context->current_large_type==0) | |||
| 359 | instr=INSTR_LARGE_DIRECT6U; | |||
| 360 | else if (xtc3_context->current_large_type==1) | |||
| 361 | instr=INSTR_LARGE_INTRA_DELTA7U; | |||
| 362 | else | |||
| 363 | instr=INSTR_LARGE_INTER_DELTA8U; | |||
| 364 | insert_value_in_array(&xtc3_context->instructions, | |||
| 365 | &xtc3_context->ninstr, | |||
| 366 | &xtc3_context->ninstr_alloc, | |||
| 367 | instr,"instr"); | |||
| 368 | } | |||
| 369 | } | |||
| 370 | ||||
| 371 | static void write_three_large(struct xtc3_context *xtc3_context, | |||
| 372 | int i) | |||
| 373 | { | |||
| 374 | int m; | |||
| 375 | if (xtc3_context->current_large_type==0) | |||
| 376 | { | |||
| 377 | for (m=0; m<3; m++) | |||
| 378 | insert_value_in_array(&xtc3_context->large_direct, | |||
| 379 | &xtc3_context->nlargedir, | |||
| 380 | &xtc3_context->nlargedir_alloc, | |||
| 381 | xtc3_context->has_large_ints[i*3+m],"large direct"); | |||
| 382 | } | |||
| 383 | else if (xtc3_context->current_large_type==1) | |||
| 384 | { | |||
| 385 | for (m=0; m<3; m++) | |||
| 386 | insert_value_in_array(&xtc3_context->large_intra_delta, | |||
| 387 | &xtc3_context->nlargeintra, | |||
| 388 | &xtc3_context->nlargeintra_alloc, | |||
| 389 | xtc3_context->has_large_ints[i*3+m],"large intra"); | |||
| 390 | } | |||
| 391 | else | |||
| 392 | { | |||
| 393 | for (m=0; m<3; m++) | |||
| 394 | insert_value_in_array(&xtc3_context->large_inter_delta, | |||
| 395 | &xtc3_context->nlargeinter, | |||
| 396 | &xtc3_context->nlargeinter_alloc, | |||
| 397 | xtc3_context->has_large_ints[i*3+m],"large inter"); | |||
| 398 | } | |||
| 399 | } | |||
| 400 | ||||
| 401 | static void flush_large(struct xtc3_context *xtc3_context, | |||
| 402 | int n) /* How many to flush. */ | |||
| 403 | { | |||
| 404 | int i; | |||
| 405 | i=0; | |||
| 406 | while (i<n) | |||
| 407 | { | |||
| 408 | int j,k; | |||
| 409 | /* If the first large is of a different kind than the currently used we must | |||
| 410 | emit an "instruction" to change the large type. */ | |||
| 411 | large_instruction_change(xtc3_context,i); | |||
| 412 | /* How many large of the same kind in a row? */ | |||
| 413 | for (j=0; | |||
| 414 | (i+j<n) && | |||
| 415 | (xtc3_context->has_large_type[i+j]==xtc3_context->has_large_type[i]); | |||
| 416 | j++); | |||
| 417 | if (j<3) | |||
| 418 | { | |||
| 419 | for (k=0; k<j; k++) | |||
| 420 | { | |||
| 421 | insert_value_in_array(&xtc3_context->instructions, | |||
| 422 | &xtc3_context->ninstr, | |||
| 423 | &xtc3_context->ninstr_alloc, | |||
| 424 | INSTR_ONLY_LARGE2U,"instr"); | |||
| 425 | write_three_large(xtc3_context,i+k); | |||
| 426 | } | |||
| 427 | } | |||
| 428 | else | |||
| 429 | { | |||
| 430 | insert_value_in_array(&xtc3_context->instructions, | |||
| 431 | &xtc3_context->ninstr, | |||
| 432 | &xtc3_context->ninstr_alloc, | |||
| 433 | INSTR_LARGE_RLE5U,"instr"); | |||
| 434 | insert_value_in_array(&xtc3_context->rle, | |||
| 435 | &xtc3_context->nrle, | |||
| 436 | &xtc3_context->nrle_alloc, | |||
| 437 | (unsigned int)j,"rle (large)"); | |||
| 438 | for (k=0; k<j; k++) | |||
| 439 | write_three_large(xtc3_context,i+k); | |||
| 440 | } | |||
| 441 | i+=j; | |||
| 442 | } | |||
| 443 | if ((xtc3_context->has_large-n)!=0) | |||
| 444 | { | |||
| 445 | int j; | |||
| 446 | for (i=0; i<xtc3_context->has_large-n; i++) | |||
| 447 | { | |||
| 448 | xtc3_context->has_large_type[i]=xtc3_context->has_large_type[i+n]; | |||
| 449 | for (j=0; j<3; j++) | |||
| 450 | xtc3_context->has_large_ints[i*3+j]=xtc3_context->has_large_ints[(i+n)*3+j]; | |||
| 451 | } | |||
| 452 | } | |||
| 453 | xtc3_context->has_large-=n; /* Number of remaining large atoms in buffer */ | |||
| 454 | } | |||
| 455 | ||||
| 456 | static double compute_intlen(unsigned int *ints) | |||
| 457 | { | |||
| 458 | /* The largest value. */ | |||
| 459 | unsigned int m=ints[0]; | |||
| 460 | if (ints[1]>m) | |||
| 461 | m=ints[1]; | |||
| 462 | if (ints[2]>m) | |||
| 463 | m=ints[2]; | |||
| 464 | return (double)m; | |||
| 465 | } | |||
| 466 | ||||
| 467 | static void buffer_large(struct xtc3_context *xtc3_context, int *input, int inpdata, | |||
| 468 | int natoms, int intradelta_ok) | |||
| 469 | { | |||
| 470 | unsigned int direct[3], intradelta[3]={0,}, interdelta[3]={0,}; | |||
| 471 | double minlen; | |||
| 472 | int best_type; | |||
| 473 | int frame=inpdata/(natoms*3); | |||
| 474 | int atomframe=inpdata%(natoms*3); | |||
| 475 | /* If it is full we must write them all. */ | |||
| 476 | if (xtc3_context->has_large==MAX_LARGE_RLE1024) | |||
| 477 | flush_large(xtc3_context,xtc3_context->has_large); /* Flush all. */ | |||
| 478 | /* Find out which is the best choice for the large integer. Direct coding, or some | |||
| 479 | kind of delta coding? */ | |||
| 480 | /* First create direct coding. */ | |||
| 481 | direct[0]=(unsigned int)(input[inpdata]-xtc3_context->minint[0]); | |||
| 482 | direct[1]=(unsigned int)(input[inpdata+1]-xtc3_context->minint[1]); | |||
| 483 | direct[2]=(unsigned int)(input[inpdata+2]-xtc3_context->minint[2]); | |||
| 484 | minlen=compute_intlen(direct); | |||
| 485 | best_type=0; /* Direct. */ | |||
| 486 | #if 1 | |||
| 487 | /* Then try intra coding if we can. */ | |||
| 488 | if ((intradelta_ok) && (atomframe>=3)) | |||
| 489 | { | |||
| 490 | double thislen; | |||
| 491 | intradelta[0]=positive_int(input[inpdata]-input[inpdata-3]); | |||
| 492 | intradelta[1]=positive_int(input[inpdata+1]-input[inpdata-2]); | |||
| 493 | intradelta[2]=positive_int(input[inpdata+2]-input[inpdata-1]); | |||
| 494 | thislen=compute_intlen(intradelta); | |||
| 495 | if (thislen*TRESHOLD_INTRA_INTER_DIRECT1.5<minlen) | |||
| 496 | { | |||
| 497 | minlen=thislen; | |||
| 498 | best_type=1; /* Intra delta */ | |||
| 499 | } | |||
| 500 | } | |||
| 501 | #endif | |||
| 502 | #if 1 | |||
| 503 | /* Then try inter coding if we can. */ | |||
| 504 | if (frame>0) | |||
| 505 | { | |||
| 506 | double thislen; | |||
| 507 | interdelta[0]=positive_int(input[inpdata]-input[inpdata-natoms*3]); | |||
| 508 | interdelta[1]=positive_int(input[inpdata+1]-input[inpdata-natoms*3+1]); | |||
| 509 | interdelta[2]=positive_int(input[inpdata+2]-input[inpdata-natoms*3+2]); | |||
| 510 | thislen=compute_intlen(interdelta); | |||
| 511 | if (thislen*TRESHOLD_INTRA_INTER_DIRECT1.5<minlen) | |||
| 512 | { | |||
| 513 | best_type=2; /* Inter delta */ | |||
| 514 | } | |||
| 515 | } | |||
| 516 | #endif | |||
| 517 | xtc3_context->has_large_type[xtc3_context->has_large]=best_type; | |||
| 518 | if (best_type==0) | |||
| 519 | { | |||
| 520 | xtc3_context->has_large_ints[xtc3_context->has_large*3]=direct[0]; | |||
| 521 | xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=direct[1]; | |||
| 522 | xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=direct[2]; | |||
| 523 | } | |||
| 524 | else if (best_type==1) | |||
| 525 | { | |||
| 526 | xtc3_context->has_large_ints[xtc3_context->has_large*3]=intradelta[0]; | |||
| 527 | xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=intradelta[1]; | |||
| 528 | xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=intradelta[2]; | |||
| 529 | } | |||
| 530 | else if (best_type==2) | |||
| 531 | { | |||
| 532 | xtc3_context->has_large_ints[xtc3_context->has_large*3]=interdelta[0]; | |||
| 533 | xtc3_context->has_large_ints[xtc3_context->has_large*3+1]=interdelta[1]; | |||
| 534 | xtc3_context->has_large_ints[xtc3_context->has_large*3+2]=interdelta[2]; | |||
| 535 | } | |||
| 536 | xtc3_context->has_large++; | |||
| 537 | } | |||
| 538 | ||||
| 539 | static void output_int(unsigned char *output,int *outdata, unsigned int n) | |||
| 540 | { | |||
| 541 | output[(*outdata)++]=((unsigned int)n)&0xFFU; | |||
| 542 | output[(*outdata)++]=(((unsigned int)n)>>8)&0xFFU; | |||
| 543 | output[(*outdata)++]=(((unsigned int)n)>>16)&0xFFU; | |||
| 544 | output[(*outdata)++]=(((unsigned int)n)>>24)&0xFFU; | |||
| 545 | } | |||
| 546 | ||||
| 547 | #if 0 | |||
| 548 | static void printarray(unsigned int *a, int n, char *name) | |||
| 549 | { | |||
| 550 | int i; | |||
| 551 | for (i=0; i<n; i++) | |||
| 552 | fprintf(stderrstderr,"%u %s\n",a[i],name); | |||
| 553 | } | |||
| 554 | #endif | |||
| 555 | ||||
| 556 | /* The base_compress routine first compresses all x coordinates, then | |||
| 557 | y and finally z. The bases used for each can be different. The | |||
| 558 | MAXBASEVALS value determines how many coordinates are compressed | |||
| 559 | into a single number. Only resulting whole bytes are dealt with for | |||
| 560 | simplicity. MAXMAXBASEVALS is the insanely large value to accept | |||
| 561 | files written with that value. BASEINTERVAL determines how often a | |||
| 562 | new base is actually computed and stored in the output | |||
| 563 | file. MAXBASEVALS*BASEINTERVAL values are stored using the same | |||
| 564 | base in BASEINTERVAL different integers. Note that the primarily | |||
| 565 | the decompression using a large MAXBASEVALS becomes very slow. */ | |||
| 566 | #define MAXMAXBASEVALS16384U 16384U | |||
| 567 | #define MAXBASEVALS24U 24U | |||
| 568 | #define BASEINTERVAL8 8 | |||
| 569 | ||||
| 570 | /* How many bytes are needed to store n values in base base */ | |||
| 571 | static int base_bytes(unsigned int base, int n) | |||
| 572 | { | |||
| 573 | int i,j; | |||
| 574 | unsigned int largeint[MAXMAXBASEVALS16384U+1]; | |||
| 575 | unsigned int largeint_tmp[MAXMAXBASEVALS16384U+1]; | |||
| 576 | int numbytes=0; | |||
| 577 | for (i=0; i<n+1; i++) | |||
| 578 | largeint[i]=0U; | |||
| 579 | for (i=0; i<n; i++) | |||
| 580 | { | |||
| 581 | if (i!=0) | |||
| 582 | { | |||
| 583 | Ptngc_largeint_mul(base,largeint,largeint_tmp,n+1); | |||
| 584 | for (j=0; j<n+1; j++) | |||
| 585 | largeint[j]=largeint_tmp[j]; | |||
| 586 | } | |||
| 587 | Ptngc_largeint_add(base-1U,largeint,n+1); | |||
| 588 | } | |||
| 589 | for (i=0; i<n; i++) | |||
| 590 | if (largeint[i]) | |||
| 591 | for (j=0; j<4; j++) | |||
| 592 | if ((largeint[i]>>(j*8))&0xFFU) | |||
| 593 | numbytes=i*4+j+1; | |||
| 594 | return numbytes; | |||
| 595 | } | |||
| 596 | ||||
| 597 | static void base_compress(unsigned int *data, int len, unsigned char *output, int *outlen) | |||
| 598 | { | |||
| 599 | unsigned int largeint[MAXBASEVALS24U+1]; | |||
| 600 | unsigned int largeint_tmp[MAXBASEVALS24U+1]; | |||
| 601 | int ixyz, i; | |||
| 602 | unsigned int j; | |||
| 603 | int nwrittenout=0; | |||
| 604 | unsigned int numbytes=0; | |||
| 605 | /* Store the MAXBASEVALS value in the output. */ | |||
| 606 | output[nwrittenout++]=(unsigned char)(MAXBASEVALS24U&0xFFU); | |||
| 607 | output[nwrittenout++]=(unsigned char)((MAXBASEVALS24U>>8)&0xFFU); | |||
| 608 | /* Store the BASEINTERVAL value in the output. */ | |||
| 609 | output[nwrittenout++]=(unsigned char)(BASEINTERVAL8&0xFFU); | |||
| 610 | for (ixyz=0; ixyz<3; ixyz++) | |||
| 611 | { | |||
| 612 | unsigned int base=0U; | |||
| 613 | int nvals=0; | |||
| 614 | int basegiven=0; | |||
| 615 | for (j=0; j<MAXBASEVALS24U+1; j++) | |||
| 616 | largeint[j]=0U; | |||
| 617 | for (i=ixyz; i<len; i+=3) | |||
| 618 | { | |||
| 619 | if (nvals==0) | |||
| 620 | { | |||
| 621 | int basecheckvals=0; | |||
| 622 | int k; | |||
| 623 | if (basegiven==0) | |||
| 624 | { | |||
| 625 | base=0U; | |||
| 626 | /* Find the largest value for this particular coordinate. */ | |||
| 627 | for (k=i; k<len; k+=3) | |||
| 628 | { | |||
| 629 | if (data[k]>base) | |||
| 630 | base=data[k]; | |||
| 631 | basecheckvals++; | |||
| 632 | if (basecheckvals==MAXBASEVALS24U*BASEINTERVAL8) | |||
| 633 | break; | |||
| 634 | } | |||
| 635 | /* The base is one larger than the largest values. */ | |||
| 636 | base++; | |||
| 637 | if (base<2) | |||
| 638 | base=2; | |||
| 639 | /* Store the base in the output. */ | |||
| 640 | output[nwrittenout++]=(unsigned char)(base&0xFFU); | |||
| 641 | output[nwrittenout++]=(unsigned char)((base>>8)&0xFFU); | |||
| 642 | output[nwrittenout++]=(unsigned char)((base>>16)&0xFFU); | |||
| 643 | output[nwrittenout++]=(unsigned char)((base>>24)&0xFFU); | |||
| 644 | basegiven=BASEINTERVAL8; | |||
| 645 | /* How many bytes is needed to store MAXBASEVALS values using this base? */ | |||
| 646 | numbytes=base_bytes(base,MAXBASEVALS24U); | |||
| 647 | } | |||
| 648 | basegiven--; | |||
| 649 | #ifdef SHOWIT | |||
| 650 | fprintf(stderrstderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,MAXBASEVALS24U); | |||
| 651 | #endif | |||
| 652 | } | |||
| 653 | if (nvals!=0) | |||
| 654 | { | |||
| 655 | Ptngc_largeint_mul(base,largeint,largeint_tmp,MAXBASEVALS24U+1); | |||
| 656 | for (j=0; j<MAXBASEVALS24U+1; j++) | |||
| 657 | largeint[j]=largeint_tmp[j]; | |||
| 658 | } | |||
| 659 | Ptngc_largeint_add(data[i],largeint,MAXBASEVALS24U+1); | |||
| 660 | #ifdef SHOWIT | |||
| 661 | fprintf(stderrstderr,"outputting value %u\n",data[i]); | |||
| 662 | #endif | |||
| 663 | nvals++; | |||
| 664 | if (nvals==MAXBASEVALS24U) | |||
| 665 | { | |||
| 666 | #ifdef SHOWIT | |||
| 667 | fprintf(stderrstderr,"Writing largeint: "); | |||
| 668 | #endif | |||
| 669 | for (j=0; j<numbytes; j++) | |||
| 670 | { | |||
| 671 | int ilarge=j/4; | |||
| 672 | int ibyte=j%4; | |||
| 673 | output[nwrittenout++]=(unsigned char)((largeint[ilarge]>>(ibyte*8))&(0xFFU)); | |||
| 674 | #ifdef SHOWIT | |||
| 675 | fprintf(stderrstderr,"%02x",(unsigned int)output[nwrittenout-1]); | |||
| 676 | #endif | |||
| 677 | } | |||
| 678 | #ifdef SHOWIT | |||
| 679 | fprintf(stderrstderr,"\n"); | |||
| 680 | #endif | |||
| 681 | nvals=0; | |||
| 682 | for (j=0; j<MAXBASEVALS24U+1; j++) | |||
| 683 | largeint[j]=0U; | |||
| 684 | } | |||
| 685 | } | |||
| 686 | if (nvals) | |||
| 687 | { | |||
| 688 | numbytes=base_bytes(base,nvals); | |||
| 689 | #ifdef SHOWIT | |||
| 690 | fprintf(stderrstderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,nvals); | |||
| 691 | #endif | |||
| 692 | for (j=0; j<numbytes; j++) | |||
| 693 | { | |||
| 694 | int ilarge=j/4; | |||
| 695 | int ibyte=j%4; | |||
| 696 | output[nwrittenout++]=(unsigned char)((largeint[ilarge]>>(ibyte*8))&(0xFFU)); | |||
| 697 | } | |||
| 698 | } | |||
| 699 | } | |||
| 700 | *outlen=nwrittenout; | |||
| 701 | } | |||
| 702 | ||||
| 703 | static void base_decompress(unsigned char *input, int len, unsigned int *output) | |||
| 704 | { | |||
| 705 | unsigned int largeint[MAXMAXBASEVALS16384U+1]; | |||
| 706 | unsigned int largeint_tmp[MAXMAXBASEVALS16384U+1]; | |||
| 707 | int ixyz, i, j; | |||
| 708 | int maxbasevals=(int)((unsigned int)(input[0])|(((unsigned int)(input[1]))<<8)); | |||
| 709 | int baseinterval=(int)input[2]; | |||
| 710 | if (maxbasevals>(int)MAXMAXBASEVALS16384U) | |||
| 711 | { | |||
| 712 | fprintf(stderrstderr,"Read a larger maxbasevals value from the file than I can handle. Fix" | |||
| 713 | " by increasing MAXMAXBASEVALS to at least %d. Although, this is" | |||
| 714 | " probably a bug in TRAJNG, since MAXMAXBASEVALS should already be insanely large enough.\n",maxbasevals); | |||
| 715 | exit(EXIT_FAILURE1); | |||
| 716 | } | |||
| 717 | input+=3; | |||
| 718 | for (ixyz=0; ixyz<3; ixyz++) | |||
| 719 | { | |||
| 720 | int numbytes=0; | |||
| 721 | int nvals_left=len/3; | |||
| 722 | int outvals=ixyz; | |||
| 723 | int basegiven=0; | |||
| 724 | unsigned int base=0U; | |||
| 725 | #ifdef SHOWIT | |||
| 726 | fprintf(stderrstderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,maxbasevals); | |||
| 727 | #endif | |||
| 728 | while (nvals_left) | |||
| 729 | { | |||
| 730 | int n; | |||
| 731 | if (basegiven==0) | |||
| 732 | { | |||
| 733 | base=(unsigned int)(input[0])| | |||
| 734 | (((unsigned int)(input[1]))<<8)| | |||
| 735 | (((unsigned int)(input[2]))<<16)| | |||
| 736 | (((unsigned int)(input[3]))<<24); | |||
| 737 | input+=4; | |||
| 738 | basegiven=baseinterval; | |||
| 739 | /* How many bytes is needed to store maxbasevals values using this base? */ | |||
| 740 | numbytes=base_bytes(base,maxbasevals); | |||
| 741 | } | |||
| 742 | basegiven--; | |||
| 743 | if (nvals_left<maxbasevals) | |||
| 744 | { | |||
| 745 | numbytes=base_bytes(base,nvals_left); | |||
| 746 | #ifdef SHOWIT | |||
| 747 | fprintf(stderrstderr,"Base for %d is %u. I need %d bytes for %d values.\n",ixyz,base,numbytes,nvals_left); | |||
| 748 | #endif | |||
| 749 | } | |||
| 750 | for (j=0; j<maxbasevals+1; j++) | |||
| 751 | largeint[j]=0U; | |||
| 752 | #ifdef SHOWIT | |||
| 753 | fprintf(stderrstderr,"Reading largeint: "); | |||
| 754 | #endif | |||
| 755 | if (numbytes/4 < maxbasevals+1) | |||
| 756 | { | |||
| 757 | for (j=0; j<numbytes; j++) | |||
| 758 | { | |||
| 759 | int ilarge=j/4; | |||
| 760 | int ibyte=j%4; | |||
| 761 | largeint[ilarge]|=((unsigned int)input[j])<<(ibyte*8); | |||
| ||||
| 762 | #ifdef SHOWIT | |||
| 763 | fprintf(stderrstderr,"%02x",(unsigned int)input[j]); | |||
| 764 | #endif | |||
| 765 | } | |||
| 766 | } | |||
| 767 | #ifdef SHOWIT | |||
| 768 | fprintf(stderrstderr,"\n"); | |||
| 769 | #endif | |||
| 770 | input+=numbytes; | |||
| 771 | /* Do the long division required to get the output values. */ | |||
| 772 | n=maxbasevals; | |||
| 773 | if (n>nvals_left) | |||
| 774 | n=nvals_left; | |||
| 775 | for (i=n-1; i>=0; i--) | |||
| 776 | { | |||
| 777 | output[outvals+i*3]=Ptngc_largeint_div(base,largeint,largeint_tmp,maxbasevals+1); | |||
| 778 | for (j=0; j<maxbasevals+1; j++) | |||
| 779 | largeint[j]=largeint_tmp[j]; | |||
| 780 | } | |||
| 781 | #ifdef SHOWIT | |||
| 782 | for (i=0; i<n; i++) | |||
| 783 | fprintf(stderrstderr,"outputting value %u\n",output[outvals+i*3]); | |||
| 784 | #endif | |||
| 785 | outvals+=n*3; | |||
| 786 | nvals_left-=n; | |||
| 787 | } | |||
| 788 | } | |||
| 789 | } | |||
| 790 | ||||
| 791 | /* If a large proportion of the integers are large (More than 10\% are >14 bits) we return 0, otherwise 1 */ | |||
| 792 | static int heuristic_bwlzh(unsigned int *ints, int nints) | |||
| 793 | { | |||
| 794 | int i,num; | |||
| 795 | num=0; | |||
| 796 | for (i=0; i<nints; i++) | |||
| 797 | if (ints[i]>=16384) | |||
| 798 | num++; | |||
| 799 | if (num>nints/10) | |||
| 800 | return 0; | |||
| 801 | else | |||
| 802 | return 1; | |||
| 803 | } | |||
| 804 | ||||
| 805 | /* Speed selects how careful to try to find the most efficient compression. The BWLZH algo is expensive! | |||
| 806 | Speed <=2 always avoids BWLZH everywhere it is possible. | |||
| 807 | Speed 3 and 4 and 5 use heuristics (check proportion of large value). This should mostly be safe. | |||
| 808 | Speed 5 enables the LZ77 component of BWLZH. | |||
| 809 | Speed 6 always tests if BWLZH is better and if it is uses it. This can be very slow. | |||
| 810 | */ | |||
| 811 | unsigned char *Ptngc_pack_array_xtc3(int *input, int *length, int natoms, int speed) | |||
| 812 | { | |||
| 813 | unsigned char *output=NULL((void*)0); | |||
| 814 | int i,ienc,j; | |||
| 815 | int outdata=0; | |||
| 816 | /* Pack triplets. */ | |||
| 817 | int ntriplets=*length/3; | |||
| 818 | int intmax; | |||
| 819 | int max_small; | |||
| 820 | int small_index; | |||
| 821 | int max_large_index; | |||
| 822 | int large_index[3]; | |||
| 823 | int prevcoord[3]; | |||
| 824 | int runlength=0; /* Initial runlength. "Stupidly" set to zero for | |||
| 825 | simplicity and explicity */ | |||
| 826 | int swapatoms=0; /* Initial guess is that we should not swap the | |||
| 827 | first two atoms in each large+small | |||
| 828 | transition */ | |||
| 829 | int didswap; /* Whether swapping was actually done. */ | |||
| 830 | int inpdata=0; | |||
| 831 | int encode_ints[3+MAX_SMALL_RLE12*3]; /* Up to 3 large + 24 small ints can be encoded at once */ | |||
| 832 | int nencode; | |||
| 833 | int ntriplets_left=ntriplets; | |||
| 834 | int refused=0; | |||
| 835 | unsigned char *bwlzh_buf=NULL((void*)0); | |||
| 836 | int bwlzh_buf_len; | |||
| 837 | unsigned char *base_buf=NULL((void*)0); | |||
| 838 | int base_buf_len; | |||
| 839 | ||||
| 840 | struct xtc3_context xtc3_context; | |||
| 841 | init_xtc3_context(&xtc3_context); | |||
| 842 | ||||
| 843 | xtc3_context.maxint[0]=xtc3_context.minint[0]=input[0]; | |||
| 844 | xtc3_context.maxint[1]=xtc3_context.minint[1]=input[1]; | |||
| 845 | xtc3_context.maxint[2]=xtc3_context.minint[2]=input[2]; | |||
| 846 | ||||
| 847 | /* Values of speed should be sane. */ | |||
| 848 | if (speed<1) | |||
| 849 | speed=1; | |||
| 850 | if (speed>6) | |||
| 851 | speed=6; | |||
| 852 | ||||
| 853 | #ifdef SHOWIT | |||
| 854 | nbits_sum=0; | |||
| 855 | nvalues_sum=0; | |||
| 856 | #endif | |||
| 857 | /* Allocate enough memory for output */ | |||
| 858 | if (*length < 48) | |||
| 859 | output=warnmalloc(8*48*sizeof *output)Ptngc_warnmalloc_x(8*48*sizeof *output,"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,859); | |||
| 860 | else | |||
| 861 | output=warnmalloc(8* *length*sizeof *output)Ptngc_warnmalloc_x(8* *length*sizeof *output,"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,861); | |||
| 862 | ||||
| 863 | ||||
| 864 | for (i=1; i<ntriplets; i++) | |||
| 865 | for (j=0; j<3; j++) | |||
| 866 | { | |||
| 867 | if (input[i*3+j]>xtc3_context.maxint[j]) | |||
| 868 | xtc3_context.maxint[j]=input[i*3+j]; | |||
| 869 | if (input[i*3+j]<xtc3_context.minint[j]) | |||
| 870 | xtc3_context.minint[j]=input[i*3+j]; | |||
| 871 | } | |||
| 872 | ||||
| 873 | large_index[0]=Ptngc_find_magic_index(xtc3_context.maxint[0]-xtc3_context.minint[0]+1); | |||
| 874 | large_index[1]=Ptngc_find_magic_index(xtc3_context.maxint[1]-xtc3_context.minint[1]+1); | |||
| 875 | large_index[2]=Ptngc_find_magic_index(xtc3_context.maxint[2]-xtc3_context.minint[2]+1); | |||
| 876 | max_large_index=large_index[0]; | |||
| 877 | if (large_index[1]>max_large_index) | |||
| 878 | max_large_index=large_index[1]; | |||
| 879 | if (large_index[2]>max_large_index) | |||
| 880 | max_large_index=large_index[2]; | |||
| 881 | ||||
| 882 | #ifdef SHOWIT | |||
| 883 | for (j=0; j<3; j++) | |||
| 884 | fprintf(stderrstderr,"minint[%d]=%d. maxint[%d]=%d large_index[%d]=%d value=%d\n",j,xtc3_context.minint[j],j,xtc3_context.maxint[j], | |||
| 885 | j,large_index[j],Ptngc_magic(large_index[j])); | |||
| 886 | #endif | |||
| 887 | ||||
| 888 | /* Guess initial small index */ | |||
| 889 | small_index=max_large_index/2; | |||
| 890 | ||||
| 891 | /* Find the largest value that is not large. Not large is half index of | |||
| 892 | large. */ | |||
| 893 | max_small=Ptngc_magic(small_index); | |||
| 894 | intmax=0; | |||
| 895 | for (i=0; i<*length; i++) | |||
| 896 | { | |||
| 897 | int item=input[i]; | |||
| 898 | int s=positive_int(item); | |||
| 899 | if (s>intmax) | |||
| 900 | if (s<max_small) | |||
| 901 | intmax=s; | |||
| 902 | } | |||
| 903 | /* This value is not critical, since if I guess wrong, the code will | |||
| 904 | just insert instructions to increase this value. */ | |||
| 905 | small_index=Ptngc_find_magic_index(intmax); | |||
| 906 | #ifdef SHOWIT | |||
| 907 | fprintf(stderrstderr,"initial small_index=%d value=%d\n",small_index,Ptngc_magic(small_index)); | |||
| 908 | #endif | |||
| 909 | ||||
| 910 | output_int(output,&outdata,positive_int(xtc3_context.minint[0])); | |||
| 911 | output_int(output,&outdata,positive_int(xtc3_context.minint[1])); | |||
| 912 | output_int(output,&outdata,positive_int(xtc3_context.minint[2])); | |||
| 913 | ||||
| 914 | #if 0 | |||
| 915 | #ifdef SHOWIT | |||
| 916 | for (i=0; i<ntriplets_left; i++) | |||
| 917 | fprintf(stderrstderr,"VALUE:%d %6d %6d %6d\n", | |||
| 918 | i, | |||
| 919 | input[inpdata+i*3], | |||
| 920 | input[inpdata+i*3+1], | |||
| 921 | input[inpdata+i*3+2]); | |||
| 922 | #endif | |||
| 923 | #endif | |||
| 924 | ||||
| 925 | /* Initial prevcoord is the minimum integers. */ | |||
| 926 | prevcoord[0]=xtc3_context.minint[0]; | |||
| 927 | prevcoord[1]=xtc3_context.minint[1]; | |||
| 928 | prevcoord[2]=xtc3_context.minint[2]; | |||
| 929 | ||||
| 930 | while (ntriplets_left) | |||
| 931 | { | |||
| 932 | if (ntriplets_left<0) | |||
| 933 | { | |||
| 934 | fprintf(stderrstderr,"TRAJNG: BUG! ntriplets_left<0!\n"); | |||
| 935 | exit(EXIT_FAILURE1); | |||
| 936 | } | |||
| 937 | /* If only less than three atoms left we just write them all as large integers. Here no swapping is done! */ | |||
| 938 | if (ntriplets_left<3) | |||
| 939 | { | |||
| 940 | for (ienc=0; ienc<ntriplets_left; ienc++) | |||
| 941 | { | |||
| 942 | buffer_large(&xtc3_context,input,inpdata,natoms,1); | |||
| 943 | inpdata+=3; | |||
| 944 | ntriplets_left--; | |||
| 945 | } | |||
| 946 | flush_large(&xtc3_context,xtc3_context.has_large); /* Flush all */ | |||
| 947 | } | |||
| 948 | else | |||
| 949 | { | |||
| 950 | int min_runlength=0; | |||
| 951 | int largest_required_base; | |||
| 952 | int largest_runlength_base; | |||
| 953 | int largest_required_index; | |||
| 954 | int largest_runlength_index; | |||
| 955 | int new_runlength; | |||
| 956 | int new_small_index; | |||
| 957 | int iter_runlength; | |||
| 958 | int iter_small_index; | |||
| 959 | int rle_index_dep; | |||
| 960 | didswap=0; | |||
| 961 | /* Insert the next batch of integers to be encoded into the buffer */ | |||
| 962 | #ifdef SHOWIT | |||
| 963 | fprintf(stderrstderr,"Initial batch\n"); | |||
| 964 | #endif | |||
| 965 | insert_batch(input+inpdata,ntriplets_left,prevcoord,encode_ints,0,&nencode); | |||
| 966 | ||||
| 967 | /* First we must decide if the next value is large (does not reasonably fit in current small encoding) | |||
| 968 | Also, if we have not written any values yet, we must begin by writing a large atom. */ | |||
| 969 | if ((inpdata==0) || (is_quite_large(encode_ints,small_index,max_large_index)) || (refused)) | |||
| 970 | { | |||
| 971 | /* If any of the next two atoms are large we should probably write them as large and not swap them */ | |||
| 972 | int no_swap=0; | |||
| 973 | if ((is_quite_large(encode_ints+3,small_index,max_large_index)) || (is_quite_large(encode_ints+6,small_index,max_large_index))) | |||
| 974 | no_swap=1; | |||
| 975 | #if 1 | |||
| 976 | if (!no_swap) | |||
| 977 | { | |||
| 978 | /* If doing inter-frame coding results in smaller values we should not do any swapping either. */ | |||
| 979 | int frame=inpdata/(natoms*3); | |||
| 980 | if (frame>0) | |||
| 981 | { | |||
| 982 | unsigned int delta[3], delta2[3]; | |||
| 983 | delta[0]=positive_int(input[inpdata+3]-input[inpdata-natoms*3+3]); | |||
| 984 | delta[1]=positive_int(input[inpdata+4]-input[inpdata-natoms*3+4]); | |||
| 985 | delta[2]=positive_int(input[inpdata+5]-input[inpdata-natoms*3+5]); | |||
| 986 | delta2[0]=positive_int(encode_ints[3]); | |||
| 987 | delta2[1]=positive_int(encode_ints[4]); | |||
| 988 | delta2[2]=positive_int(encode_ints[5]); | |||
| 989 | #ifdef SHOWIT | |||
| 990 | fprintf(stderrstderr,"A1: inter delta: %u %u %u. intra delta=%u %u %u\n", | |||
| 991 | delta[0],delta[1],delta[2], | |||
| 992 | delta2[0],delta2[1],delta2[2]); | |||
| 993 | #endif | |||
| 994 | if (compute_intlen(delta)*TRESHOLD_INTER_INTRA5.0<compute_intlen(delta2)) | |||
| 995 | { | |||
| 996 | delta[0]=positive_int(input[inpdata+6]-input[inpdata-natoms*3+6]); | |||
| 997 | delta[1]=positive_int(input[inpdata+7]-input[inpdata-natoms*3+7]); | |||
| 998 | delta[2]=positive_int(input[inpdata+8]-input[inpdata-natoms*3+8]); | |||
| 999 | delta2[0]=positive_int(encode_ints[6]); | |||
| 1000 | delta2[1]=positive_int(encode_ints[7]); | |||
| 1001 | delta2[2]=positive_int(encode_ints[8]); | |||
| 1002 | #ifdef SHOWIT | |||
| 1003 | fprintf(stderrstderr,"A2: inter delta: %u %u %u. intra delta=%u %u %u\n", | |||
| 1004 | delta[0],delta[1],delta[2], | |||
| 1005 | delta2[0],delta2[1],delta2[2]); | |||
| 1006 | #endif | |||
| 1007 | if (compute_intlen(delta)*TRESHOLD_INTER_INTRA5.0<compute_intlen(delta2)) | |||
| 1008 | { | |||
| 1009 | no_swap=1; | |||
| 1010 | #ifdef SHOWIT | |||
| 1011 | fprintf(stderrstderr,"SETTING NO SWAP!\n"); | |||
| 1012 | #endif | |||
| 1013 | } | |||
| 1014 | } | |||
| 1015 | } | |||
| 1016 | } | |||
| 1017 | #endif | |||
| 1018 | if (!no_swap) | |||
| 1019 | { | |||
| 1020 | /* Next we must decide if we should swap the first | |||
| 1021 | two values. */ | |||
| 1022 | #if 1 | |||
| 1023 | swapdecide(&xtc3_context,input+inpdata,&swapatoms,large_index,xtc3_context.minint); | |||
| 1024 | #else | |||
| 1025 | swapatoms=0; | |||
| 1026 | #endif | |||
| 1027 | /* If we should do the integer swapping manipulation we should do it now */ | |||
| 1028 | if (swapatoms) | |||
| 1029 | { | |||
| 1030 | didswap=1; | |||
| 1031 | for (i=0; i<3; i++) | |||
| 1032 | { | |||
| 1033 | int in[3], out[3]; | |||
| 1034 | in[0]=input[inpdata+i]; | |||
| 1035 | in[1]=input[inpdata+3+i]-input[inpdata+i]; | |||
| 1036 | in[2]=input[inpdata+6+i]-input[inpdata+3+i]; | |||
| 1037 | swap_ints(in,out); | |||
| 1038 | encode_ints[i]=out[0]; | |||
| 1039 | encode_ints[3+i]=out[1]; | |||
| 1040 | encode_ints[6+i]=out[2]; | |||
| 1041 | } | |||
| 1042 | /* We have swapped atoms, so the minimum run-length is 2 */ | |||
| 1043 | #ifdef SHOWIT | |||
| 1044 | fprintf(stderrstderr,"Swap atoms results in:\n"); | |||
| 1045 | for (i=0; i<3; i++) | |||
| 1046 | fprintf(stderrstderr,"%d: %6d %6d %6d\t\t%6d %6d %6d\n",i*3, | |||
| 1047 | encode_ints[i*3], | |||
| 1048 | encode_ints[i*3+1], | |||
| 1049 | encode_ints[i*3+2], | |||
| 1050 | positive_int(encode_ints[i*3]), | |||
| 1051 | positive_int(encode_ints[i*3+1]), | |||
| 1052 | positive_int(encode_ints[i*3+2])); | |||
| 1053 | ||||
| 1054 | #endif | |||
| 1055 | min_runlength=2; | |||
| 1056 | } | |||
| 1057 | } | |||
| 1058 | /* Cache large value for later possible combination with | |||
| 1059 | a sequence of small integers. */ | |||
| 1060 | if ((swapatoms) && (didswap)) | |||
| 1061 | { | |||
| 1062 | buffer_large(&xtc3_context,input,inpdata+3,natoms,0); /* This is a swapped integer, so inpdata is one atom later and | |||
| 1063 | intra coding is not ok. */ | |||
| 1064 | for (ienc=0; ienc<3; ienc++) | |||
| 1065 | prevcoord[ienc]=input[inpdata+3+ienc]; | |||
| 1066 | } | |||
| 1067 | else | |||
| 1068 | { | |||
| 1069 | buffer_large(&xtc3_context,input,inpdata,natoms,1); | |||
| 1070 | for (ienc=0; ienc<3; ienc++) | |||
| 1071 | prevcoord[ienc]=input[inpdata+ienc]; | |||
| 1072 | } | |||
| 1073 | ||||
| 1074 | ||||
| 1075 | #ifdef SHOWIT | |||
| 1076 | fprintf(stderrstderr,"Prevcoord after packing of large: %d %d %d\n", | |||
| 1077 | prevcoord[0],prevcoord[1],prevcoord[2]); | |||
| 1078 | #endif | |||
| 1079 | ||||
| 1080 | /* We have written a large integer so we have one less atoms to worry about */ | |||
| 1081 | inpdata+=3; | |||
| 1082 | ntriplets_left--; | |||
| 1083 | ||||
| 1084 | refused=0; | |||
| 1085 | ||||
| 1086 | /* Insert the next batch of integers to be encoded into the buffer */ | |||
| 1087 | #ifdef SHOWIT | |||
| 1088 | fprintf(stderrstderr,"Update batch due to large int.\n"); | |||
| 1089 | #endif | |||
| 1090 | if ((swapatoms) && (didswap)) | |||
| 1091 | { | |||
| 1092 | /* Keep swapped values. */ | |||
| 1093 | for (i=0; i<2; i++) | |||
| 1094 | for (ienc=0; ienc<3; ienc++) | |||
| 1095 | encode_ints[i*3+ienc]=encode_ints[(i+1)*3+ienc]; | |||
| 1096 | } | |||
| 1097 | insert_batch(input+inpdata,ntriplets_left,prevcoord,encode_ints,min_runlength,&nencode); | |||
| 1098 | } | |||
| 1099 | /* Here we should only have differences for the atom coordinates. */ | |||
| 1100 | /* Convert the ints to positive ints */ | |||
| 1101 | for (ienc=0; ienc<nencode; ienc++) | |||
| 1102 | { | |||
| 1103 | int pint=positive_int(encode_ints[ienc]); | |||
| 1104 | encode_ints[ienc]=pint; | |||
| 1105 | } | |||
| 1106 | /* Now we must decide what base and runlength to do. If we have swapped atoms it will be at least 2. | |||
| 1107 | If even the next atom is large, we will not do anything. */ | |||
| 1108 | largest_required_base=0; | |||
| 1109 | /* Determine required base */ | |||
| 1110 | for (ienc=0; ienc<min_runlength*3; ienc++) | |||
| 1111 | if (encode_ints[ienc]>largest_required_base) | |||
| 1112 | largest_required_base=encode_ints[ienc]; | |||
| 1113 | /* Also compute what the largest base is for the current runlength setting! */ | |||
| 1114 | largest_runlength_base=0; | |||
| 1115 | for (ienc=0; (ienc<runlength*3) && (ienc<nencode); ienc++) | |||
| 1116 | if (encode_ints[ienc]>largest_runlength_base) | |||
| 1117 | largest_runlength_base=encode_ints[ienc]; | |||
| 1118 | ||||
| 1119 | largest_required_index=Ptngc_find_magic_index(largest_required_base); | |||
| 1120 | largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base); | |||
| 1121 | ||||
| 1122 | if (largest_required_index<largest_runlength_index) | |||
| 1123 | { | |||
| 1124 | new_runlength=min_runlength; | |||
| 1125 | new_small_index=largest_required_index; | |||
| 1126 | } | |||
| 1127 | else | |||
| 1128 | { | |||
| 1129 | new_runlength=runlength; | |||
| 1130 | new_small_index=largest_runlength_index; | |||
| 1131 | } | |||
| 1132 | ||||
| 1133 | /* Only allow increase of runlength wrt min_runlength */ | |||
| 1134 | if (new_runlength<min_runlength) | |||
| 1135 | new_runlength=min_runlength; | |||
| 1136 | ||||
| 1137 | /* If the current runlength is longer than the number of | |||
| 1138 | triplets left stop it from being so. */ | |||
| 1139 | if (new_runlength>ntriplets_left) | |||
| 1140 | new_runlength=ntriplets_left; | |||
| 1141 | ||||
| 1142 | /* We must at least try to get some small integers going. */ | |||
| 1143 | if (new_runlength==0) | |||
| 1144 | { | |||
| 1145 | new_runlength=1; | |||
| 1146 | new_small_index=small_index; | |||
| 1147 | } | |||
| 1148 | ||||
| 1149 | iter_runlength=new_runlength; | |||
| 1150 | iter_small_index=new_small_index; | |||
| 1151 | ||||
| 1152 | /* Iterate to find optimal encoding and runlength */ | |||
| 1153 | #ifdef SHOWIT | |||
| 1154 | fprintf(stderrstderr,"Entering iterative loop.\n"); | |||
| 1155 | fflush(stderrstderr); | |||
| 1156 | #endif | |||
| 1157 | ||||
| 1158 | do { | |||
| 1159 | new_runlength=iter_runlength; | |||
| 1160 | new_small_index=iter_small_index; | |||
| 1161 | ||||
| 1162 | #ifdef SHOWIT | |||
| 1163 | fprintf(stderrstderr,"Test new_small_index=%d Base=%d\n",new_small_index,Ptngc_magic(new_small_index)); | |||
| 1164 | #endif | |||
| 1165 | /* What is the largest runlength | |||
| 1166 | we can do with the currently | |||
| 1167 | selected encoding? Also the max supported runlength is MAX_SMALL_RLE triplets! */ | |||
| 1168 | for (ienc=0; ienc<nencode && ienc<MAX_SMALL_RLE12*3; ienc++) | |||
| 1169 | { | |||
| 1170 | int test_index=Ptngc_find_magic_index(encode_ints[ienc]); | |||
| 1171 | if (test_index>new_small_index) | |||
| 1172 | break; | |||
| 1173 | } | |||
| 1174 | if (ienc/3>new_runlength) | |||
| 1175 | { | |||
| 1176 | iter_runlength=ienc/3; | |||
| 1177 | #ifdef SHOWIT | |||
| 1178 | fprintf(stderrstderr,"I found a new possible runlength: %d\n",iter_runlength); | |||
| 1179 | #endif | |||
| 1180 | } | |||
| 1181 | ||||
| 1182 | /* How large encoding do we have to use? */ | |||
| 1183 | largest_runlength_base=0; | |||
| 1184 | for (ienc=0; ienc<iter_runlength*3; ienc++) | |||
| 1185 | if (encode_ints[ienc]>largest_runlength_base) | |||
| 1186 | largest_runlength_base=encode_ints[ienc]; | |||
| 1187 | largest_runlength_index=Ptngc_find_magic_index(largest_runlength_base); | |||
| 1188 | if (largest_runlength_index!=new_small_index) | |||
| 1189 | { | |||
| 1190 | iter_small_index=largest_runlength_index; | |||
| 1191 | #ifdef SHOWIT | |||
| 1192 | fprintf(stderrstderr,"I found a new possible small index: %d Base=%d\n",iter_small_index,Ptngc_magic(iter_small_index)); | |||
| 1193 | #endif | |||
| 1194 | } | |||
| 1195 | } while ((new_runlength!=iter_runlength) || | |||
| 1196 | (new_small_index!=iter_small_index)); | |||
| 1197 | ||||
| 1198 | #ifdef SHOWIT | |||
| 1199 | fprintf(stderrstderr,"Exit iterative loop.\n"); | |||
| 1200 | fflush(stderrstderr); | |||
| 1201 | #endif | |||
| 1202 | ||||
| 1203 | /* Verify that we got something good. We may have caught a | |||
| 1204 | substantially larger atom. If so we should just bail | |||
| 1205 | out and let the loop get on another lap. We may have a | |||
| 1206 | minimum runlength though and then we have to fulfill | |||
| 1207 | the request to write out these atoms! */ | |||
| 1208 | rle_index_dep=0; | |||
| 1209 | if (new_runlength<3) | |||
| 1210 | rle_index_dep=IS_LARGE6; | |||
| 1211 | else if (new_runlength<6) | |||
| 1212 | rle_index_dep=QUITE_LARGE3; | |||
| 1213 | if ((min_runlength) | |||
| 1214 | || ((new_small_index<small_index+IS_LARGE6) && (new_small_index+rle_index_dep<max_large_index)) | |||
| 1215 | #if 1 | |||
| 1216 | || (new_small_index+IS_LARGE6<max_large_index) | |||
| 1217 | #endif | |||
| 1218 | ) | |||
| 1219 | { | |||
| 1220 | /* If doing inter-frame coding of large integers results | |||
| 1221 | in smaller values than the small value we should not | |||
| 1222 | produce a sequence of small values here. */ | |||
| 1223 | int frame=inpdata/(natoms*3); | |||
| 1224 | int numsmaller=0; | |||
| 1225 | #if 1 | |||
| 1226 | if ((!swapatoms) && (frame>0)) | |||
| 1227 | { | |||
| 1228 | for (i=0; i<new_runlength; i++) | |||
| 1229 | { | |||
| 1230 | unsigned int delta[3]; | |||
| 1231 | unsigned int delta2[3]; | |||
| 1232 | delta[0]=positive_int(input[inpdata+i*3]-input[inpdata-natoms*3+i*3]); | |||
| 1233 | delta[1]=positive_int(input[inpdata+i*3+1]-input[inpdata-natoms*3+i*3+1]); | |||
| 1234 | delta[2]=positive_int(input[inpdata+i*3+2]-input[inpdata-natoms*3+i*3+2]); | |||
| 1235 | delta2[0]=positive_int(encode_ints[i*3]); | |||
| 1236 | delta2[1]=positive_int(encode_ints[i*3+1]); | |||
| 1237 | delta2[2]=positive_int(encode_ints[i*3+2]); | |||
| 1238 | if (compute_intlen(delta)*TRESHOLD_INTER_INTRA5.0<compute_intlen(delta2)) | |||
| 1239 | numsmaller++; | |||
| 1240 | } | |||
| 1241 | } | |||
| 1242 | #endif | |||
| 1243 | /* Most of the values should become smaller, otherwise | |||
| 1244 | we should encode them with intra coding. */ | |||
| 1245 | if ((!swapatoms) && (numsmaller>=2*new_runlength/3)) | |||
| 1246 | { | |||
| 1247 | /* Put all the values in large arrays, instead of the small array */ | |||
| 1248 | if (new_runlength) | |||
| 1249 | { | |||
| 1250 | for (i=0; i<new_runlength; i++) | |||
| 1251 | buffer_large(&xtc3_context,input,inpdata+i*3,natoms,1); | |||
| 1252 | for (i=0; i<3; i++) | |||
| 1253 | prevcoord[i]=input[inpdata+(new_runlength-1)*3+i]; | |||
| 1254 | inpdata+=3*new_runlength; | |||
| 1255 | ntriplets_left-=new_runlength; | |||
| 1256 | } | |||
| 1257 | } | |||
| 1258 | else | |||
| 1259 | { | |||
| 1260 | if ((new_runlength!=runlength) || (new_small_index!=small_index)) | |||
| 1261 | { | |||
| 1262 | int change=new_small_index-small_index; | |||
| 1263 | ||||
| 1264 | if (new_small_index<=0) | |||
| 1265 | change=0; | |||
| 1266 | ||||
| 1267 | if (change<0) | |||
| 1268 | { | |||
| 1269 | int ixx; | |||
| 1270 | for (ixx=0; ixx<new_runlength; ixx++) | |||
| 1271 | { | |||
| 1272 | int rejected; | |||
| 1273 | do { | |||
| 1274 | int ixyz; | |||
| 1275 | double isum=0.; /* ints can be almost 32 bit so multiplication will overflow. So do doubles. */ | |||
| 1276 | for (ixyz=0; ixyz<3; ixyz++) | |||
| 1277 | { | |||
| 1278 | /* encode_ints is already positive (and multiplied by 2 versus the original, just as magic ints) */ | |||
| 1279 | double id=encode_ints[ixx*3+ixyz]; | |||
| 1280 | isum+=id*id; | |||
| 1281 | } | |||
| 1282 | rejected=0; | |||
| 1283 | #ifdef SHOWIT | |||
| 1284 | fprintf(stderrstderr,"Tested decrease %d of index: %g>=%g?\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change)); | |||
| 1285 | #endif | |||
| 1286 | if (isum>(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change)) | |||
| 1287 | { | |||
| 1288 | #ifdef SHOWIT | |||
| 1289 | fprintf(stderrstderr,"Rejected decrease %d of index due to length of vector: %g>=%g\n",change,isum,(double)Ptngc_magic(small_index+change)*(double)Ptngc_magic(small_index+change)); | |||
| 1290 | #endif | |||
| 1291 | rejected=1; | |||
| 1292 | change++; | |||
| 1293 | } | |||
| 1294 | } while ((change<0) && (rejected)); | |||
| 1295 | if (change==0) | |||
| 1296 | break; | |||
| 1297 | } | |||
| 1298 | } | |||
| 1299 | ||||
| 1300 | /* Always accept the new small indices here. */ | |||
| 1301 | small_index=new_small_index; | |||
| 1302 | /* If we have a new runlength emit it */ | |||
| 1303 | if (runlength!=new_runlength) | |||
| 1304 | { | |||
| 1305 | runlength=new_runlength; | |||
| 1306 | insert_value_in_array(&xtc3_context.instructions, | |||
| 1307 | &xtc3_context.ninstr, | |||
| 1308 | &xtc3_context.ninstr_alloc, | |||
| 1309 | INSTR_SMALL_RUNLENGTH1U,"instr"); | |||
| 1310 | insert_value_in_array(&xtc3_context.rle, | |||
| 1311 | &xtc3_context.nrle, | |||
| 1312 | &xtc3_context.nrle_alloc, | |||
| 1313 | (unsigned int)runlength,"rle (small)"); | |||
| 1314 | } | |||
| 1315 | #ifdef SHOWIT | |||
| 1316 | fprintf(stderrstderr,"Current small index: %d Base=%d\n",small_index,Ptngc_magic(small_index)); | |||
| 1317 | #endif | |||
| 1318 | } | |||
| 1319 | /* If we have a large previous integer we can combine it with a sequence of small ints. */ | |||
| 1320 | if (xtc3_context.has_large) | |||
| 1321 | { | |||
| 1322 | /* If swapatoms is set to 1 but we did actually not | |||
| 1323 | do any swapping, we must first write out the | |||
| 1324 | large atom and then the small. If swapatoms is 1 | |||
| 1325 | and we did swapping we can use the efficient | |||
| 1326 | encoding. */ | |||
| 1327 | if ((swapatoms) && (!didswap)) | |||
| 1328 | { | |||
| 1329 | #ifdef SHOWIT | |||
| 1330 | fprintf(stderrstderr,"Swapatoms was set to 1 but we did not do swapping!\n"); | |||
| 1331 | fprintf(stderrstderr,"Only one large integer.\n"); | |||
| 1332 | #endif | |||
| 1333 | /* Flush all large atoms. */ | |||
| 1334 | flush_large(&xtc3_context,xtc3_context.has_large); | |||
| 1335 | #ifdef SHOWIT | |||
| 1336 | fprintf(stderrstderr,"Sequence of only small integers.\n"); | |||
| 1337 | #endif | |||
| 1338 | insert_value_in_array(&xtc3_context.instructions, | |||
| 1339 | &xtc3_context.ninstr, | |||
| 1340 | &xtc3_context.ninstr_alloc, | |||
| 1341 | INSTR_ONLY_SMALL3U,"instr"); | |||
| 1342 | } | |||
| 1343 | else | |||
| 1344 | { | |||
| 1345 | ||||
| 1346 | #ifdef SHOWIT | |||
| 1347 | fprintf(stderrstderr,"Sequence of one large and small integers (good compression).\n"); | |||
| 1348 | #endif | |||
| 1349 | /* Flush all large atoms but one! */ | |||
| 1350 | if (xtc3_context.has_large>1) | |||
| 1351 | flush_large(&xtc3_context,xtc3_context.has_large-1); | |||
| 1352 | ||||
| 1353 | /* Here we must check if we should emit a large | |||
| 1354 | type change instruction. */ | |||
| 1355 | large_instruction_change(&xtc3_context,0); | |||
| 1356 | ||||
| 1357 | insert_value_in_array(&xtc3_context.instructions, | |||
| 1358 | &xtc3_context.ninstr, | |||
| 1359 | &xtc3_context.ninstr_alloc, | |||
| 1360 | INSTR_DEFAULT0U,"instr"); | |||
| 1361 | ||||
| 1362 | write_three_large(&xtc3_context,0); | |||
| 1363 | xtc3_context.has_large=0; | |||
| 1364 | } | |||
| 1365 | } | |||
| 1366 | else | |||
| 1367 | { | |||
| 1368 | #ifdef SHOWIT | |||
| 1369 | fprintf(stderrstderr,"Sequence of only small integers.\n"); | |||
| 1370 | #endif | |||
| 1371 | insert_value_in_array(&xtc3_context.instructions, | |||
| 1372 | &xtc3_context.ninstr, | |||
| 1373 | &xtc3_context.ninstr_alloc, | |||
| 1374 | INSTR_ONLY_SMALL3U,"instr"); | |||
| 1375 | } | |||
| 1376 | /* Insert the small integers into the small integer array. */ | |||
| 1377 | for (ienc=0; ienc<runlength*3; ienc++) | |||
| 1378 | insert_value_in_array(&xtc3_context.smallintra, | |||
| 1379 | &xtc3_context.nsmallintra, | |||
| 1380 | &xtc3_context.nsmallintra_alloc, | |||
| 1381 | (unsigned int)encode_ints[ienc],"smallintra"); | |||
| 1382 | ||||
| 1383 | #ifdef SHOWIT | |||
| 1384 | for (ienc=0; ienc<runlength; ienc++) | |||
| 1385 | fprintf(stderrstderr,"Small: %d %d %d\n", | |||
| 1386 | encode_ints[ienc*3], | |||
| 1387 | encode_ints[ienc*3+1], | |||
| 1388 | encode_ints[ienc*3+2]); | |||
| 1389 | #endif | |||
| 1390 | /* Update prevcoord. */ | |||
| 1391 | for (ienc=0; ienc<runlength; ienc++) | |||
| 1392 | { | |||
| 1393 | #ifdef SHOWIT | |||
| 1394 | fprintf(stderrstderr,"Prevcoord in packing: %d %d %d\n", | |||
| 1395 | prevcoord[0],prevcoord[1],prevcoord[2]); | |||
| 1396 | #endif | |||
| 1397 | prevcoord[0]+=unpositive_int(encode_ints[ienc*3]); | |||
| 1398 | prevcoord[1]+=unpositive_int(encode_ints[ienc*3+1]); | |||
| 1399 | prevcoord[2]+=unpositive_int(encode_ints[ienc*3+2]); | |||
| 1400 | } | |||
| 1401 | #ifdef SHOWIT | |||
| 1402 | fprintf(stderrstderr,"Prevcoord in packing: %d %d %d\n", | |||
| 1403 | prevcoord[0],prevcoord[1],prevcoord[2]); | |||
| 1404 | #endif | |||
| 1405 | ||||
| 1406 | inpdata+=3*runlength; | |||
| 1407 | ntriplets_left-=runlength; | |||
| 1408 | #if 1 | |||
| 1409 | } | |||
| 1410 | #endif | |||
| 1411 | } | |||
| 1412 | else | |||
| 1413 | { | |||
| 1414 | #ifdef SHOWIT | |||
| 1415 | fprintf(stderrstderr,"Refused value: %d old is %d max is %d\n",new_small_index,small_index,max_large_index); | |||
| 1416 | fflush(stderrstderr); | |||
| 1417 | #endif | |||
| 1418 | refused=1; | |||
| 1419 | } | |||
| 1420 | } | |||
| 1421 | #ifdef SHOWIT | |||
| 1422 | fprintf(stderrstderr,"Number of triplets left is %d\n",ntriplets_left); | |||
| 1423 | #endif | |||
| 1424 | } | |||
| 1425 | ||||
| 1426 | /* If we have large previous integers we must flush them now. */ | |||
| 1427 | if (xtc3_context.has_large) | |||
| 1428 | flush_large(&xtc3_context,xtc3_context.has_large); | |||
| 1429 | ||||
| 1430 | /* Now it is time to compress all the data in the buffers with the bwlzh or base algo. */ | |||
| 1431 | ||||
| 1432 | #if 0 | |||
| 1433 | /* Inspect the data. */ | |||
| 1434 | printarray(xtc3_context.instructions,xtc3_context.ninstr,"A instr"); | |||
| 1435 | printarray(xtc3_context.rle,xtc3_context.nrle,"A rle"); | |||
| 1436 | printarray(xtc3_context.large_direct,xtc3_context.nlargedir,"A largedir"); | |||
| 1437 | printarray(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,"A largeintra"); | |||
| 1438 | printarray(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,"A largeinter"); | |||
| 1439 | printarray(xtc3_context.smallintra,xtc3_context.nsmallintra,"A smallintra"); | |||
| 1440 | exit(0); | |||
| 1441 | #endif | |||
| 1442 | ||||
| 1443 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1444 | fprintf(stderrstderr,"instructions: %d\n",xtc3_context.ninstr); | |||
| 1445 | #endif | |||
| 1446 | ||||
| 1447 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1448 | #define bwlzh_compress bwlzh_compress_verbose | |||
| 1449 | #define bwlzh_compress_no_lz77 bwlzh_compress_no_lz77_verbose | |||
| 1450 | #endif | |||
| 1451 | ||||
| 1452 | output_int(output,&outdata,(unsigned int)xtc3_context.ninstr); | |||
| 1453 | if (xtc3_context.ninstr) | |||
| 1454 | { | |||
| 1455 | bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.ninstr))Ptngc_warnmalloc_x(bwlzh_get_buflen(xtc3_context.ninstr),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1455); | |||
| 1456 | if (speed>=5) | |||
| 1457 | bwlzh_compress(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len); | |||
| 1458 | else | |||
| 1459 | bwlzh_compress_no_lz77(xtc3_context.instructions,xtc3_context.ninstr,bwlzh_buf,&bwlzh_buf_len); | |||
| 1460 | output_int(output,&outdata,(unsigned int)bwlzh_buf_len); | |||
| 1461 | memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); | |||
| 1462 | outdata+=bwlzh_buf_len; | |||
| 1463 | free(bwlzh_buf); | |||
| 1464 | } | |||
| 1465 | ||||
| 1466 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1467 | fprintf(stderrstderr,"rle: %d\n",xtc3_context.nrle); | |||
| 1468 | #endif | |||
| 1469 | ||||
| 1470 | output_int(output,&outdata,(unsigned int)xtc3_context.nrle); | |||
| 1471 | if (xtc3_context.nrle) | |||
| 1472 | { | |||
| 1473 | bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nrle))Ptngc_warnmalloc_x(bwlzh_get_buflen(xtc3_context.nrle),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1473); | |||
| 1474 | if (speed>=5) | |||
| 1475 | bwlzh_compress(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len); | |||
| 1476 | else | |||
| 1477 | bwlzh_compress_no_lz77(xtc3_context.rle,xtc3_context.nrle,bwlzh_buf,&bwlzh_buf_len); | |||
| 1478 | output_int(output,&outdata,(unsigned int)bwlzh_buf_len); | |||
| 1479 | memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); | |||
| 1480 | outdata+=bwlzh_buf_len; | |||
| 1481 | free(bwlzh_buf); | |||
| 1482 | } | |||
| 1483 | ||||
| 1484 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1485 | fprintf(stderrstderr,"large direct: %d\n",xtc3_context.nlargedir); | |||
| 1486 | #endif | |||
| 1487 | ||||
| 1488 | output_int(output,&outdata,(unsigned int)xtc3_context.nlargedir); | |||
| 1489 | if (xtc3_context.nlargedir) | |||
| 1490 | { | |||
| 1491 | if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_direct,xtc3_context.nlargedir)))) | |||
| 1492 | { | |||
| 1493 | bwlzh_buf=NULL((void*)0); | |||
| 1494 | bwlzh_buf_len=INT_MAX2147483647; | |||
| 1495 | } | |||
| 1496 | else | |||
| 1497 | { | |||
| 1498 | bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargedir))Ptngc_warnmalloc_x(bwlzh_get_buflen(xtc3_context.nlargedir),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1498); | |||
| 1499 | if (speed>=5) | |||
| 1500 | bwlzh_compress(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len); | |||
| 1501 | else | |||
| 1502 | bwlzh_compress_no_lz77(xtc3_context.large_direct,xtc3_context.nlargedir,bwlzh_buf,&bwlzh_buf_len); | |||
| 1503 | } | |||
| 1504 | /* If this can be written smaller using base compression we should do that. */ | |||
| 1505 | base_buf=warnmalloc((xtc3_context.nlargedir+3)*sizeof(int))Ptngc_warnmalloc_x((xtc3_context.nlargedir+3)*sizeof(int),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1505); | |||
| 1506 | base_compress(xtc3_context.large_direct,xtc3_context.nlargedir,base_buf,&base_buf_len); | |||
| 1507 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1508 | fprintf(stderrstderr,"Large direct: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); | |||
| 1509 | #endif | |||
| 1510 | if (base_buf_len<bwlzh_buf_len) | |||
| 1511 | { | |||
| 1512 | output[outdata++]=0U; | |||
| 1513 | output_int(output,&outdata,(unsigned int)base_buf_len); | |||
| 1514 | memcpy(output+outdata,base_buf,base_buf_len); | |||
| 1515 | outdata+=base_buf_len; | |||
| 1516 | } | |||
| 1517 | else | |||
| 1518 | { | |||
| 1519 | output[outdata++]=1U; | |||
| 1520 | output_int(output,&outdata,(unsigned int)bwlzh_buf_len); | |||
| 1521 | memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); | |||
| 1522 | outdata+=bwlzh_buf_len; | |||
| 1523 | } | |||
| 1524 | free(bwlzh_buf); | |||
| 1525 | free(base_buf); | |||
| 1526 | } | |||
| 1527 | ||||
| 1528 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1529 | fprintf(stderrstderr,"large intra: %d\n",xtc3_context.nlargeintra); | |||
| 1530 | #endif | |||
| 1531 | ||||
| 1532 | output_int(output,&outdata,(unsigned int)xtc3_context.nlargeintra); | |||
| 1533 | if (xtc3_context.nlargeintra) | |||
| 1534 | { | |||
| 1535 | if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_intra_delta,xtc3_context.nlargeintra)))) | |||
| 1536 | { | |||
| 1537 | bwlzh_buf=NULL((void*)0); | |||
| 1538 | bwlzh_buf_len=INT_MAX2147483647; | |||
| 1539 | } | |||
| 1540 | else | |||
| 1541 | { | |||
| 1542 | bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargeintra))Ptngc_warnmalloc_x(bwlzh_get_buflen(xtc3_context.nlargeintra) ,"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1542); | |||
| 1543 | if (speed>=5) | |||
| 1544 | bwlzh_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len); | |||
| 1545 | else | |||
| 1546 | bwlzh_compress_no_lz77(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,bwlzh_buf,&bwlzh_buf_len); | |||
| 1547 | } | |||
| 1548 | /* If this can be written smaller using base compression we should do that. */ | |||
| 1549 | base_buf=warnmalloc((xtc3_context.nlargeintra+3)*sizeof(int))Ptngc_warnmalloc_x((xtc3_context.nlargeintra+3)*sizeof(int),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1549); | |||
| 1550 | base_compress(xtc3_context.large_intra_delta,xtc3_context.nlargeintra,base_buf,&base_buf_len); | |||
| 1551 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1552 | fprintf(stderrstderr,"Large intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); | |||
| 1553 | #endif | |||
| 1554 | if (base_buf_len<bwlzh_buf_len) | |||
| 1555 | { | |||
| 1556 | output[outdata++]=0U; | |||
| 1557 | output_int(output,&outdata,(unsigned int)base_buf_len); | |||
| 1558 | memcpy(output+outdata,base_buf,base_buf_len); | |||
| 1559 | outdata+=base_buf_len; | |||
| 1560 | } | |||
| 1561 | else | |||
| 1562 | { | |||
| 1563 | output[outdata++]=1U; | |||
| 1564 | output_int(output,&outdata,(unsigned int)bwlzh_buf_len); | |||
| 1565 | memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); | |||
| 1566 | outdata+=bwlzh_buf_len; | |||
| 1567 | } | |||
| 1568 | free(bwlzh_buf); | |||
| 1569 | free(base_buf); | |||
| 1570 | } | |||
| 1571 | ||||
| 1572 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1573 | fprintf(stderrstderr,"large inter: %d\n",xtc3_context.nlargeinter); | |||
| 1574 | #endif | |||
| 1575 | ||||
| 1576 | output_int(output,&outdata,(unsigned int)xtc3_context.nlargeinter); | |||
| 1577 | if (xtc3_context.nlargeinter) | |||
| 1578 | { | |||
| 1579 | if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.large_inter_delta,xtc3_context.nlargeinter)))) | |||
| 1580 | { | |||
| 1581 | bwlzh_buf=NULL((void*)0); | |||
| 1582 | bwlzh_buf_len=INT_MAX2147483647; | |||
| 1583 | } | |||
| 1584 | else | |||
| 1585 | { | |||
| 1586 | bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nlargeinter))Ptngc_warnmalloc_x(bwlzh_get_buflen(xtc3_context.nlargeinter) ,"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1586); | |||
| 1587 | if (speed>=5) | |||
| 1588 | bwlzh_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len); | |||
| 1589 | else | |||
| 1590 | bwlzh_compress_no_lz77(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,bwlzh_buf,&bwlzh_buf_len); | |||
| 1591 | } | |||
| 1592 | /* If this can be written smaller using base compression we should do that. */ | |||
| 1593 | base_buf=warnmalloc((xtc3_context.nlargeinter+3)*sizeof(int))Ptngc_warnmalloc_x((xtc3_context.nlargeinter+3)*sizeof(int),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1593); | |||
| 1594 | base_compress(xtc3_context.large_inter_delta,xtc3_context.nlargeinter,base_buf,&base_buf_len); | |||
| 1595 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1596 | fprintf(stderrstderr,"Large inter: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); | |||
| 1597 | #endif | |||
| 1598 | if (base_buf_len<bwlzh_buf_len) | |||
| 1599 | { | |||
| 1600 | output[outdata++]=0U; | |||
| 1601 | output_int(output,&outdata,(unsigned int)base_buf_len); | |||
| 1602 | memcpy(output+outdata,base_buf,base_buf_len); | |||
| 1603 | outdata+=base_buf_len; | |||
| 1604 | } | |||
| 1605 | else | |||
| 1606 | { | |||
| 1607 | output[outdata++]=1U; | |||
| 1608 | output_int(output,&outdata,(unsigned int)bwlzh_buf_len); | |||
| 1609 | memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); | |||
| 1610 | outdata+=bwlzh_buf_len; | |||
| 1611 | } | |||
| 1612 | free(bwlzh_buf); | |||
| 1613 | free(base_buf); | |||
| 1614 | } | |||
| 1615 | ||||
| 1616 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1617 | fprintf(stderrstderr,"small intra: %d\n",xtc3_context.nsmallintra); | |||
| 1618 | #endif | |||
| 1619 | ||||
| 1620 | output_int(output,&outdata,(unsigned int)xtc3_context.nsmallintra); | |||
| 1621 | if (xtc3_context.nsmallintra) | |||
| 1622 | { | |||
| 1623 | if ((speed<=2) || ((speed<=5) && (!heuristic_bwlzh(xtc3_context.smallintra,xtc3_context.nsmallintra)))) | |||
| 1624 | { | |||
| 1625 | bwlzh_buf=NULL((void*)0); | |||
| 1626 | bwlzh_buf_len=INT_MAX2147483647; | |||
| 1627 | } | |||
| 1628 | else | |||
| 1629 | { | |||
| 1630 | bwlzh_buf=warnmalloc(bwlzh_get_buflen(xtc3_context.nsmallintra))Ptngc_warnmalloc_x(bwlzh_get_buflen(xtc3_context.nsmallintra) ,"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1630); | |||
| 1631 | if (speed>=5) | |||
| 1632 | bwlzh_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len); | |||
| 1633 | else | |||
| 1634 | bwlzh_compress_no_lz77(xtc3_context.smallintra,xtc3_context.nsmallintra,bwlzh_buf,&bwlzh_buf_len); | |||
| 1635 | } | |||
| 1636 | /* If this can be written smaller using base compression we should do that. */ | |||
| 1637 | base_buf=warnmalloc((xtc3_context.nsmallintra+3)*sizeof(int))Ptngc_warnmalloc_x((xtc3_context.nsmallintra+3)*sizeof(int),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1637); | |||
| 1638 | base_compress(xtc3_context.smallintra,xtc3_context.nsmallintra,base_buf,&base_buf_len); | |||
| 1639 | #if defined(SHOWIT) || defined(SHOWIT_LIGHT) | |||
| 1640 | fprintf(stderrstderr,"Small intra: Base len=%d. BWLZH len=%d\n",base_buf_len,bwlzh_buf_len); | |||
| 1641 | #endif | |||
| 1642 | if (base_buf_len<bwlzh_buf_len) | |||
| 1643 | { | |||
| 1644 | output[outdata++]=0U; | |||
| 1645 | output_int(output,&outdata,(unsigned int)base_buf_len); | |||
| 1646 | memcpy(output+outdata,base_buf,base_buf_len); | |||
| 1647 | outdata+=base_buf_len; | |||
| 1648 | } | |||
| 1649 | else | |||
| 1650 | { | |||
| 1651 | output[outdata++]=1U; | |||
| 1652 | output_int(output,&outdata,(unsigned int)bwlzh_buf_len); | |||
| 1653 | memcpy(output+outdata,bwlzh_buf,bwlzh_buf_len); | |||
| 1654 | outdata+=bwlzh_buf_len; | |||
| 1655 | } | |||
| 1656 | free(bwlzh_buf); | |||
| 1657 | free(base_buf); | |||
| 1658 | } | |||
| 1659 | *length=outdata; | |||
| 1660 | ||||
| 1661 | free_xtc3_context(&xtc3_context); | |||
| 1662 | return output; | |||
| 1663 | } | |||
| 1664 | ||||
| 1665 | static void decompress_bwlzh_block(unsigned char **ptr, | |||
| 1666 | int nvals, | |||
| 1667 | unsigned int **vals) | |||
| 1668 | { | |||
| 1669 | int bwlzh_buf_len=(int)(((unsigned int)(*ptr)[0]) | | |||
| 1670 | (((unsigned int)(*ptr)[1])<<8) | | |||
| 1671 | (((unsigned int)(*ptr)[2])<<16) | | |||
| 1672 | (((unsigned int)(*ptr)[3])<<24)); | |||
| 1673 | (*ptr)+=4; | |||
| 1674 | *vals=warnmalloc(nvals*sizeof (**vals))Ptngc_warnmalloc_x(nvals*sizeof (**vals),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1674); | |||
| 1675 | bwlzh_decompress(*ptr,nvals,*vals); | |||
| 1676 | (*ptr)+=bwlzh_buf_len; | |||
| 1677 | } | |||
| 1678 | ||||
| 1679 | static void decompress_base_block(unsigned char **ptr, | |||
| 1680 | int nvals, | |||
| 1681 | unsigned int **vals) | |||
| 1682 | { | |||
| 1683 | int base_buf_len=(int)(((unsigned int)(*ptr)[0]) | | |||
| 1684 | (((unsigned int)(*ptr)[1])<<8) | | |||
| 1685 | (((unsigned int)(*ptr)[2])<<16) | | |||
| 1686 | (((unsigned int)(*ptr)[3])<<24)); | |||
| 1687 | (*ptr)+=4; | |||
| 1688 | *vals=warnmalloc(nvals*sizeof (**vals))Ptngc_warnmalloc_x(nvals*sizeof (**vals),"/home/alexxy/Develop/gromacs/src/external/tng_io/src/compression/xtc3.c" ,1688); | |||
| 1689 | base_decompress(*ptr,nvals,*vals); | |||
| 1690 | (*ptr)+=base_buf_len; | |||
| 1691 | } | |||
| 1692 | ||||
| 1693 | static void unpack_one_large(struct xtc3_context *xtc3_context, | |||
| 1694 | int *ilargedir, int *ilargeintra, | |||
| 1695 | int *ilargeinter, int *prevcoord, | |||
| 1696 | int *minint, int *output, | |||
| 1697 | int outdata, int didswap, | |||
| 1698 | int natoms, int current_large_type) | |||
| 1699 | { | |||
| 1700 | int large_ints[3]={0,0,0}; | |||
| 1701 | if (current_large_type==0 && xtc3_context->large_direct) | |||
| 1702 | { | |||
| 1703 | large_ints[0]=(int)xtc3_context->large_direct[(*ilargedir)]+minint[0]; | |||
| 1704 | large_ints[1]=(int)xtc3_context->large_direct[(*ilargedir)+1]+minint[1]; | |||
| 1705 | large_ints[2]=(int)xtc3_context->large_direct[(*ilargedir)+2]+minint[2]; | |||
| 1706 | (*ilargedir)+=3; | |||
| 1707 | } | |||
| 1708 | else if (current_large_type==1 && xtc3_context->large_intra_delta) | |||
| 1709 | { | |||
| 1710 | large_ints[0]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)])+prevcoord[0]; | |||
| 1711 | large_ints[1]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+1])+prevcoord[1]; | |||
| 1712 | large_ints[2]=unpositive_int(xtc3_context->large_intra_delta[(*ilargeintra)+2])+prevcoord[2]; | |||
| 1713 | (*ilargeintra)+=3; | |||
| 1714 | } | |||
| 1715 | else if (xtc3_context->large_inter_delta) | |||
| 1716 | { | |||
| 1717 | large_ints[0]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)]) | |||
| 1718 | +output[outdata-natoms*3+didswap*3]; | |||
| 1719 | large_ints[1]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+1]) | |||
| 1720 | +output[outdata-natoms*3+1+didswap*3]; | |||
| 1721 | large_ints[2]=unpositive_int(xtc3_context->large_inter_delta[(*ilargeinter)+2]) | |||
| 1722 | +output[outdata-natoms*3+2+didswap*3]; | |||
| 1723 | (*ilargeinter)+=3; | |||
| 1724 | } | |||
| 1725 | prevcoord[0]=large_ints[0]; | |||
| 1726 | prevcoord[1]=large_ints[1]; | |||
| 1727 | prevcoord[2]=large_ints[2]; | |||
| 1728 | output[outdata]=large_ints[0]; | |||
| 1729 | output[outdata+1]=large_ints[1]; | |||
| 1730 | output[outdata+2]=large_ints[2]; | |||
| 1731 | #ifdef SHOWIT | |||
| 1732 | fprintf(stderrstderr,"Unpack one large: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]); | |||
| 1733 | #endif | |||
| 1734 | } | |||
| 1735 | ||||
| 1736 | ||||
| 1737 | int Ptngc_unpack_array_xtc3(unsigned char *packed,int *output, int length, int natoms) | |||
| 1738 | { | |||
| 1739 | int i; | |||
| 1740 | int minint[3]; | |||
| 1741 | unsigned char *ptr=packed; | |||
| 1742 | int prevcoord[3]; | |||
| 1743 | int outdata=0; | |||
| 1744 | int ntriplets_left=length/3; | |||
| 1745 | int swapatoms=0; | |||
| 1746 | int runlength=0; | |||
| 1747 | int current_large_type=0; | |||
| 1748 | int iinstr=0; | |||
| 1749 | int irle=0; | |||
| 1750 | int ilargedir=0; | |||
| 1751 | int ilargeintra=0; | |||
| 1752 | int ilargeinter=0; | |||
| 1753 | int ismallintra=0; | |||
| 1754 | ||||
| 1755 | struct xtc3_context xtc3_context; | |||
| 1756 | init_xtc3_context(&xtc3_context); | |||
| 1757 | ||||
| 1758 | for (i=0; i<3; i++) | |||
| ||||
| 1759 | { | |||
| 1760 | minint[i]=unpositive_int((int)(((unsigned int)ptr[0]) | | |||
| 1761 | (((unsigned int)ptr[1])<<8) | | |||
| 1762 | (((unsigned int)ptr[2])<<16) | | |||
| 1763 | (((unsigned int)ptr[3])<<24))); | |||
| 1764 | ptr+=4; | |||
| 1765 | } | |||
| 1766 | ||||
| 1767 | xtc3_context.ninstr=(int)(((unsigned int)ptr[0]) | | |||
| 1768 | (((unsigned int)ptr[1])<<8) | | |||
| 1769 | (((unsigned int)ptr[2])<<16) | | |||
| 1770 | (((unsigned int)ptr[3])<<24)); | |||
| 1771 | ptr+=4; | |||
| 1772 | if (xtc3_context.ninstr) | |||
| 1773 | decompress_bwlzh_block(&ptr,xtc3_context.ninstr,&xtc3_context.instructions); | |||
| 1774 | ||||
| 1775 | xtc3_context.nrle=(int)(((unsigned int)ptr[0]) | | |||
| 1776 | (((unsigned int)ptr[1])<<8) | | |||
| 1777 | (((unsigned int)ptr[2])<<16) | | |||
| 1778 | (((unsigned int)ptr[3])<<24)); | |||
| 1779 | ptr+=4; | |||
| 1780 | if (xtc3_context.nrle) | |||
| 1781 | decompress_bwlzh_block(&ptr,xtc3_context.nrle,&xtc3_context.rle); | |||
| 1782 | ||||
| 1783 | xtc3_context.nlargedir=(int)(((unsigned int)ptr[0]) | | |||
| 1784 | (((unsigned int)ptr[1])<<8) | | |||
| 1785 | (((unsigned int)ptr[2])<<16) | | |||
| 1786 | (((unsigned int)ptr[3])<<24)); | |||
| 1787 | ptr+=4; | |||
| 1788 | if (xtc3_context.nlargedir) | |||
| 1789 | { | |||
| 1790 | if (*ptr++==1) | |||
| 1791 | decompress_bwlzh_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct); | |||
| 1792 | else | |||
| 1793 | decompress_base_block(&ptr,xtc3_context.nlargedir,&xtc3_context.large_direct); | |||
| 1794 | } | |||
| 1795 | ||||
| 1796 | xtc3_context.nlargeintra=(int)(((unsigned int)ptr[0]) | | |||
| 1797 | (((unsigned int)ptr[1])<<8) | | |||
| 1798 | (((unsigned int)ptr[2])<<16) | | |||
| 1799 | (((unsigned int)ptr[3])<<24)); | |||
| 1800 | ptr+=4; | |||
| 1801 | if (xtc3_context.nlargeintra) | |||
| 1802 | { | |||
| 1803 | if (*ptr++==1) | |||
| 1804 | decompress_bwlzh_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta); | |||
| 1805 | else | |||
| 1806 | decompress_base_block(&ptr,xtc3_context.nlargeintra,&xtc3_context.large_intra_delta); | |||
| 1807 | } | |||
| 1808 | ||||
| 1809 | xtc3_context.nlargeinter=(int)(((unsigned int)ptr[0]) | | |||
| 1810 | (((unsigned int)ptr[1])<<8) | | |||
| 1811 | (((unsigned int)ptr[2])<<16) | | |||
| 1812 | (((unsigned int)ptr[3])<<24)); | |||
| 1813 | ptr+=4; | |||
| 1814 | if (xtc3_context.nlargeinter) | |||
| 1815 | { | |||
| 1816 | if (*ptr++==1) | |||
| 1817 | decompress_bwlzh_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta); | |||
| 1818 | else | |||
| 1819 | decompress_base_block(&ptr,xtc3_context.nlargeinter,&xtc3_context.large_inter_delta); | |||
| 1820 | } | |||
| 1821 | ||||
| 1822 | xtc3_context.nsmallintra=(int)(((unsigned int)ptr[0]) | | |||
| 1823 | (((unsigned int)ptr[1])<<8) | | |||
| 1824 | (((unsigned int)ptr[2])<<16) | | |||
| 1825 | (((unsigned int)ptr[3])<<24)); | |||
| 1826 | ptr+=4; | |||
| 1827 | if (xtc3_context.nsmallintra) | |||
| 1828 | { | |||
| 1829 | if (*ptr++==1) | |||
| 1830 | decompress_bwlzh_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra); | |||
| 1831 | else | |||
| 1832 | decompress_base_block(&ptr,xtc3_context.nsmallintra,&xtc3_context.smallintra); | |||
| 1833 | } | |||
| 1834 | ||||
| 1835 | /* Initial prevcoord is the minimum integers. */ | |||
| 1836 | prevcoord[0]=minint[0]; | |||
| 1837 | prevcoord[1]=minint[1]; | |||
| 1838 | prevcoord[2]=minint[2]; | |||
| 1839 | ||||
| 1840 | while (ntriplets_left>0 && iinstr<xtc3_context.ninstr) | |||
| 1841 | { | |||
| 1842 | int instr=xtc3_context.instructions[iinstr++]; | |||
| 1843 | #ifdef SHOWIT | |||
| 1844 | fprintf(stderrstderr,"instr=%d @ %d\n",instr,iinstr-1); | |||
| 1845 | #endif | |||
| 1846 | #ifdef SHOWIT | |||
| 1847 | fprintf(stderrstderr,"ntriplets left=%d\n",ntriplets_left); | |||
| 1848 | #endif | |||
| 1849 | if ((instr==INSTR_DEFAULT0U) /* large+small */ | |||
| 1850 | || (instr==INSTR_ONLY_LARGE2U) /* only large */ | |||
| 1851 | || (instr==INSTR_ONLY_SMALL3U)) /* only small */ | |||
| 1852 | { | |||
| 1853 | if (instr!=INSTR_ONLY_SMALL3U) | |||
| 1854 | { | |||
| 1855 | int didswap=0; | |||
| 1856 | if ((instr==INSTR_DEFAULT0U) && (swapatoms)) | |||
| 1857 | didswap=1; | |||
| 1858 | unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter, | |||
| 1859 | prevcoord, minint, output, outdata, didswap, | |||
| 1860 | natoms, current_large_type); | |||
| 1861 | ntriplets_left--; | |||
| 1862 | outdata+=3; | |||
| 1863 | } | |||
| 1864 | if (instr!=INSTR_ONLY_LARGE2U) | |||
| 1865 | { | |||
| 1866 | for (i=0; i<runlength; i++) | |||
| 1867 | { | |||
| 1868 | prevcoord[0]+=unpositive_int(xtc3_context.smallintra[ismallintra]); | |||
| 1869 | prevcoord[1]+=unpositive_int(xtc3_context.smallintra[ismallintra+1]); | |||
| 1870 | prevcoord[2]+=unpositive_int(xtc3_context.smallintra[ismallintra+2]); | |||
| 1871 | ismallintra+=3; | |||
| 1872 | output[outdata+i*3]=prevcoord[0]; | |||
| 1873 | output[outdata+i*3+1]=prevcoord[1]; | |||
| 1874 | output[outdata+i*3+2]=prevcoord[2]; | |||
| 1875 | #ifdef SHOWIT | |||
| 1876 | fprintf(stderrstderr,"Unpack small: %d %d %d\n",prevcoord[0],prevcoord[1],prevcoord[2]); | |||
| 1877 | #endif | |||
| 1878 | } | |||
| 1879 | if ((instr==INSTR_DEFAULT0U) && (swapatoms)) | |||
| 1880 | { | |||
| 1881 | for (i=0; i<3; i++) | |||
| 1882 | { | |||
| 1883 | int tmp=output[outdata-3+i]; | |||
| 1884 | output[outdata-3+i]=output[outdata+i]; | |||
| 1885 | output[outdata+i]=tmp; | |||
| 1886 | } | |||
| 1887 | #ifdef SHOWIT | |||
| 1888 | fprintf(stderrstderr,"Unswap results in\n"); | |||
| 1889 | for (i=0; i<3; i++) | |||
| 1890 | fprintf(stderrstderr,"%d %d %d\n",output[outdata-3+i*3],output[outdata-2+i*3],output[outdata-1+i*3]); | |||
| 1891 | #endif | |||
| 1892 | } | |||
| 1893 | ntriplets_left-=runlength; | |||
| 1894 | outdata+=runlength*3; | |||
| 1895 | } | |||
| 1896 | } | |||
| 1897 | else if (instr==INSTR_LARGE_RLE5U && irle<xtc3_context.nrle) | |||
| 1898 | { | |||
| 1899 | int large_rle=xtc3_context.rle[irle++]; | |||
| 1900 | #ifdef SHOWIT | |||
| 1901 | fprintf(stderrstderr,"large_rle=%d @ %d\n",large_rle,irle-1); | |||
| 1902 | #endif | |||
| 1903 | for (i=0; i<large_rle; i++) | |||
| 1904 | { | |||
| 1905 | unpack_one_large(&xtc3_context,&ilargedir, &ilargeintra, &ilargeinter, | |||
| 1906 | prevcoord, minint, output, outdata, 0, | |||
| 1907 | natoms, current_large_type); | |||
| 1908 | ntriplets_left--; | |||
| 1909 | outdata+=3; | |||
| 1910 | } | |||
| 1911 | } | |||
| 1912 | else if (instr==INSTR_SMALL_RUNLENGTH1U && irle<xtc3_context.nrle) | |||
| 1913 | { | |||
| 1914 | runlength=xtc3_context.rle[irle++]; | |||
| 1915 | #ifdef SHOWIT | |||
| 1916 | fprintf(stderrstderr,"small_rle=%d @ %d\n",runlength,irle-1); | |||
| 1917 | #endif | |||
| 1918 | } | |||
| 1919 | else if (instr==INSTR_FLIP4U) | |||
| 1920 | { | |||
| 1921 | swapatoms=1-swapatoms; | |||
| 1922 | #ifdef SHOWIT | |||
| 1923 | fprintf(stderrstderr,"new flip=%d\n",swapatoms); | |||
| 1924 | #endif | |||
| 1925 | } | |||
| 1926 | else if (instr==INSTR_LARGE_DIRECT6U) | |||
| 1927 | { | |||
| 1928 | current_large_type=0; | |||
| 1929 | #ifdef SHOWIT | |||
| 1930 | fprintf(stderrstderr,"large direct\n"); | |||
| 1931 | #endif | |||
| 1932 | } | |||
| 1933 | else if (instr==INSTR_LARGE_INTRA_DELTA7U) | |||
| 1934 | { | |||
| 1935 | current_large_type=1; | |||
| 1936 | #ifdef SHOWIT | |||
| 1937 | fprintf(stderrstderr,"large intra delta\n"); | |||
| 1938 | #endif | |||
| 1939 | } | |||
| 1940 | else if (instr==INSTR_LARGE_INTER_DELTA8U) | |||
| 1941 | { | |||
| 1942 | current_large_type=2; | |||
| 1943 | #ifdef SHOWIT | |||
| 1944 | fprintf(stderrstderr,"large inter delta\n"); | |||
| 1945 | #endif | |||
| 1946 | } | |||
| 1947 | } | |||
| 1948 | if (ntriplets_left<0) | |||
| 1949 | { | |||
| 1950 | fprintf(stderrstderr,"TRAJNG XTC3: A bug has been found. At end ntriplets_left<0\n"); | |||
| 1951 | exit(EXIT_FAILURE1); | |||
| 1952 | } | |||
| 1953 | free_xtc3_context(&xtc3_context); | |||
| 1954 | return 0; | |||
| 1955 | } |