source: thirdparty/SZ/sz/src/sz_uint16.c @ 9ee2ce3

Revision 9ee2ce3, 40.1 KB checked in by Hal Finkel <hfinkel@…>, 6 years ago (diff)

importing new SZ files
  • Property mode set to 100644
RevLine 
[2c47b73]1/**
2 *  @file sz_uint16.c
3 *  @author Sheng Di
4 *  @date Aug, 2017
5 *  @brief sz_uint16, Compression and Decompression functions
6 *  (C) 2017 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
7 *      See COPYRIGHT in top-level directory.
8 */
9
10
11#include <stdio.h>
12#include <stdlib.h>
13#include <string.h>
14#include <unistd.h>
15#include <math.h>
16#include "sz.h"
17#include "CompressElement.h"
18#include "DynamicByteArray.h"
19#include "DynamicIntArray.h"
20#include "zlib.h"
21#include "rw.h"
22#include "TightDataPointStorageI.h"
23#include "sz_uint16.h"
[9ee2ce3]24#include "utility.h"
[2c47b73]25
26unsigned int optimize_intervals_uint16_1D(uint16_t *oriData, size_t dataLength, double realPrecision)
27{       
28        size_t i = 0, radiusIndex;
29        int64_t pred_value = 0, pred_err;
30        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
31        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
32        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
33        for(i=2;i<dataLength;i++)
34        {
35                if(i%confparams_cpr->sampleDistance==0)
36                {
37                        //pred_value = 2*oriData[i-1] - oriData[i-2];
38                        pred_value = oriData[i-1];
39                        pred_err = llabs(pred_value - oriData[i]);
40                        radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
41                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
42                                radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
43                        intervals[radiusIndex]++;
44                }
45        }
46        //compute the appropriate number
47        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
48        size_t sum = 0;
49        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
50        {
51                sum += intervals[i];
52                if(sum>targetCount)
53                        break;
54        }
55        if(i>=confparams_cpr->maxRangeRadius)
56                i = confparams_cpr->maxRangeRadius-1;
57               
58        unsigned int accIntervals = 2*(i+1);
59        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
60       
61        if(powerOf2<32)
62                powerOf2 = 32;
63       
64        free(intervals);
65        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
66        return powerOf2;
67}
68
69unsigned int optimize_intervals_uint16_2D(uint16_t *oriData, size_t r1, size_t r2, double realPrecision)
70{       
71        size_t i,j, index;
72        size_t radiusIndex;
73        int64_t pred_value = 0, pred_err;
74        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
75        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
76        size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
77        for(i=1;i<r1;i++)
78        {
79                for(j=1;j<r2;j++)
80                {
81                        if((i+j)%confparams_cpr->sampleDistance==0)
82                        {
83                                index = i*r2+j;
84                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
85                                pred_err = llabs(pred_value - oriData[index]);
86                                radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
87                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
88                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
89                                intervals[radiusIndex]++;
90                        }                       
91                }
92        }
93        //compute the appropriate number
94        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
95        size_t sum = 0;
96        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
97        {
98                sum += intervals[i];
99                if(sum>targetCount)
100                        break;
101        }
102        if(i>=confparams_cpr->maxRangeRadius)
103                i = confparams_cpr->maxRangeRadius-1;
104        unsigned int accIntervals = 2*(i+1);
105        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
106
107        if(powerOf2<32)
108                powerOf2 = 32;
109
110        free(intervals);
111        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
112        return powerOf2;
113}
114
115unsigned int optimize_intervals_uint16_3D(uint16_t *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
116{       
117        size_t i,j,k, index;
118        size_t radiusIndex;
119        size_t r23=r2*r3;
120        int64_t pred_value = 0, pred_err;
121        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
122        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
123        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
124        for(i=1;i<r1;i++)
125        {
126                for(j=1;j<r2;j++)
127                {
128                        for(k=1;k<r3;k++)
129                        {                       
130                                if((i+j+k)%confparams_cpr->sampleDistance==0)
131                                {
132                                        index = i*r23+j*r3+k;
133                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23] 
134                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
135                                        pred_err = llabs(pred_value - oriData[index]);
136                                        radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
137                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
138                                        {
139                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
140                                                //printf("radiusIndex=%d\n", radiusIndex);
141                                        }
142                                        intervals[radiusIndex]++;
143                                }
144                        }
145                }
146        }
147        //compute the appropriate number
148        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
149        size_t sum = 0;
150        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
151        {
152                sum += intervals[i];
153                if(sum>targetCount)
154                        break;
155        }
156        if(i>=confparams_cpr->maxRangeRadius)
157                i = confparams_cpr->maxRangeRadius-1;
158        unsigned int accIntervals = 2*(i+1);
159        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
160
161        if(powerOf2<32)
162                powerOf2 = 32;
163       
164        free(intervals);
165        //printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
166        return powerOf2;
167}
168
169
170unsigned int optimize_intervals_uint16_4D(uint16_t *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision)
171{
172        size_t i,j,k,l, index;
173        size_t radiusIndex;
174        size_t r234=r2*r3*r4;
175        size_t r34=r3*r4;
176        int64_t pred_value = 0, pred_err;
177        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
178        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
179        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)*(r4-1)/confparams_cpr->sampleDistance;
180        for(i=1;i<r1;i++)
181        {
182                for(j=1;j<r2;j++)
183                {
184                        for(k=1;k<r3;k++)
185                        {
186                                for (l=1;l<r4;l++)
187                                {
188                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
189                                        {
190                                                index = i*r234+j*r34+k*r4+l;
191                                                pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r34]
192                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
193                                                pred_err = llabs(pred_value - oriData[index]);
194                                                radiusIndex = (uint64_t)((pred_err/realPrecision+1)/2);
195                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
196                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
197                                                intervals[radiusIndex]++;
198                                        }
199                                }
200                        }
201                }
202        }
203        //compute the appropriate number
204        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
205        size_t sum = 0;
206        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
207        {
208                sum += intervals[i];
209                if(sum>targetCount)
210                        break;
211        }
212        if(i>=confparams_cpr->maxRangeRadius)
213                i = confparams_cpr->maxRangeRadius-1;
214
215        unsigned int accIntervals = 2*(i+1);
216        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
217
218        if(powerOf2<32)
219                powerOf2 = 32;
220
221        free(intervals);
222        return powerOf2;
223}
224
225TightDataPointStorageI* SZ_compress_uint16_1D_MDQ(uint16_t *oriData, size_t dataLength, double realPrecision, int64_t valueRangeSize, int64_t minValue)
226{
227        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
228        int byteSize = computeByteSizePerIntValue(valueRangeSize);
229        unsigned int quantization_intervals;
230        if(exe_params->optQuantMode==1)
231                quantization_intervals = optimize_intervals_uint16_1D(oriData, dataLength, realPrecision);
232        else
233                quantization_intervals = exe_params->intvCapacity;
234        updateQuantizationInfo(quantization_intervals); 
235        size_t i;
236
237        int* type = (int*) malloc(dataLength*sizeof(int));
238               
239        uint16_t* spaceFillingValue = oriData; //
240       
241        DynamicByteArray *exactDataByteArray;
242        new_DBA(&exactDataByteArray, DynArrayInitLen);
243               
244        int64_t last3CmprsData[3] = {0,0,0};
245                               
246        //add the first data   
247        type[0] = 0;
248        compressUInt16Value(spaceFillingValue[0], minValue, byteSize, bytes);
249        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
250        listAdd_int(last3CmprsData, spaceFillingValue[0]);
251               
252        type[1] = 0;
253        compressUInt16Value(spaceFillingValue[1], minValue, byteSize, bytes);
254        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
255        listAdd_int(last3CmprsData, spaceFillingValue[1]);
256        //printf("%.30G\n",last3CmprsData[0]); 
257       
258        int state;
259        double checkRadius = (exe_params->intvCapacity-1)*realPrecision;
260        int64_t curData;
261        int64_t pred, predAbsErr;
262        double interval = 2*realPrecision;
263       
264        for(i=2;i<dataLength;i++)
265        {
266                curData = spaceFillingValue[i];
267                //pred = 2*last3CmprsData[0] - last3CmprsData[1];
268                pred = last3CmprsData[0];
269                predAbsErr = llabs(curData - pred);     
[9ee2ce3]270                if(predAbsErr<checkRadius)
[2c47b73]271                {
272                        state = (predAbsErr/realPrecision+1)/2;
273                        if(curData>=pred)
274                        {
275                                type[i] = exe_params->intvRadius+state;
276                                pred = pred + state*interval;
277                        }
278                        else //curData<pred
279                        {
280                                type[i] = exe_params->intvRadius-state;
281                                pred = pred - state*interval;
282                        }
283                        if(pred>SZ_UINT16_MAX) pred = SZ_UINT16_MAX;
284                        if(pred<SZ_UINT16_MIN) pred = SZ_UINT16_MIN;                   
285                        listAdd_int(last3CmprsData, pred);                                     
286                        continue;
287                }
288               
289                //unpredictable data processing         
290                type[i] = 0;
291                compressUInt16Value(curData, minValue, byteSize, bytes);
292                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
293                listAdd_int(last3CmprsData, curData);
294        }//end of for
295               
296        size_t exactDataNum = exactDataByteArray->size / byteSize;
297       
298        TightDataPointStorageI* tdps;   
299                       
300        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
301                        type, exactDataByteArray->array, exactDataByteArray->size, 
302                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
303
304//sdi:Debug
305/*      int sum =0;
306        for(i=0;i<dataLength;i++)
307                if(type[i]==0) sum++;
308        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);*/
309       
310        //free memory
311        free(type);     
312        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
313       
314        return tdps;
315}
316
317void SZ_compress_args_uint16_StoreOriData(uint16_t* oriData, size_t dataLength, TightDataPointStorageI* tdps, 
318unsigned char** newByteData, size_t *outSize)
319{
320        int intSize=sizeof(uint16_t);   
321        size_t k = 0, i;
322        tdps->isLossless = 1;
323        size_t totalByteLength = 3 + MetaDataByteLength + exe_params->SZ_SIZE_TYPE + 1 + intSize*dataLength;
324        *newByteData = (unsigned char*)malloc(totalByteLength);
325       
326        unsigned char dsLengthBytes[8];
327        for (i = 0; i < 3; i++)//3
328                (*newByteData)[k++] = versionNumber[i]; 
329
330        if(exe_params->SZ_SIZE_TYPE==4)//1
331                (*newByteData)[k++] = 16; //00010000
332        else
333                (*newByteData)[k++] = 80;       //01010000: 01000000 indicates the SZ_SIZE_TYPE=8
334       
335        convertSZParamsToBytes(confparams_cpr, &((*newByteData)[k]));
336        k = k + MetaDataByteLength;     
337       
338        sizeToBytes(dsLengthBytes,dataLength); //SZ_SIZE_TYPE: 4 or 8   
339        for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)
340                (*newByteData)[k++] = dsLengthBytes[i];
341               
342        if(sysEndianType==BIG_ENDIAN_SYSTEM)
343                memcpy((*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, oriData, dataLength*intSize);
344        else
345        {
346                unsigned char* p = (*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
347                for(i=0;i<dataLength;i++,p+=intSize)
348                        int16ToBytes_bigEndian(p, oriData[i]);
349        }       
350        *outSize = totalByteLength;
351}
352
353void SZ_compress_args_uint16_NoCkRngeNoGzip_1D(unsigned char** newByteData, uint16_t *oriData, 
354size_t dataLength, double realPrecision, size_t *outSize, int64_t valueRangeSize, uint16_t minValue)
355{
356        TightDataPointStorageI* tdps = SZ_compress_uint16_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, minValue);
357        //TODO: return bytes....
358        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
359        if(*outSize > dataLength*sizeof(uint16_t))
360                SZ_compress_args_uint16_StoreOriData(oriData, dataLength+2, tdps, newByteData, outSize);
361        free_TightDataPointStorageI(tdps);
362}
363
364TightDataPointStorageI* SZ_compress_uint16_2D_MDQ(uint16_t *oriData, size_t r1, size_t r2, double realPrecision, int64_t valueRangeSize, int64_t minValue)
365{
366        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
367        int byteSize = computeByteSizePerIntValue(valueRangeSize);
368       
369        unsigned int quantization_intervals;
370        if(exe_params->optQuantMode==1)
371        {
372                quantization_intervals = optimize_intervals_uint16_2D(oriData, r1, r2, realPrecision);
373                updateQuantizationInfo(quantization_intervals);
374        }       
375        else
376                quantization_intervals = exe_params->intvCapacity;
377        size_t i,j; 
378        int64_t pred1D, pred2D, curValue, tmp;
379        int diff = 0.0;
380        double itvNum = 0;
381        uint16_t *P0, *P1;
382               
383        size_t dataLength = r1*r2;     
384       
385        P0 = (uint16_t*)malloc(r2*sizeof(uint16_t));
386        memset(P0, 0, r2*sizeof(uint16_t));
387        P1 = (uint16_t*)malloc(r2*sizeof(uint16_t));
388        memset(P1, 0, r2*sizeof(uint16_t));
389               
390        int* type = (int*) malloc(dataLength*sizeof(int));
391        //type[dataLength]=0;
392               
393        uint16_t* spaceFillingValue = oriData; //
394       
395        DynamicByteArray *exactDataByteArray;
396        new_DBA(&exactDataByteArray, DynArrayInitLen); 
397
398        type[0] = 0;
399        curValue = P1[0] = spaceFillingValue[0];
400        compressUInt16Value(curValue, minValue, byteSize, bytes);
401        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
402
403        /* Process Row-0 data 1*/
404        pred1D = P1[0];
405        diff = spaceFillingValue[1] - pred1D;
406
407        itvNum =  llabs(diff)/realPrecision + 1;
408
409        if (itvNum < exe_params->intvCapacity)
410        {
411                if (diff < 0) itvNum = -itvNum;
412                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
413                tmp = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
414                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
415                        P1[1] = tmp;
416                else if(tmp < SZ_UINT16_MIN)
417                        P1[1] = SZ_UINT16_MIN;
418                else
419                        P1[1] = SZ_UINT16_MAX;
420        }
421        else
422        {
423                type[1] = 0;
424                curValue = P1[1] = spaceFillingValue[1];
425                compressUInt16Value(curValue, minValue, byteSize, bytes);
426                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
427        }
428
429    /* Process Row-0 data 2 --> data r2-1 */
430        for (j = 2; j < r2; j++)
431        {
432                pred1D = 2*P1[j-1] - P1[j-2];
433                diff = spaceFillingValue[j] - pred1D;
434
435                itvNum = llabs(diff)/realPrecision + 1;
436
437                if (itvNum < exe_params->intvCapacity)
438                {
439                        if (diff < 0) itvNum = -itvNum;
440                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
441                        tmp = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
442                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
443                                P1[j] = tmp;
444                        else if(tmp < SZ_UINT16_MIN)
445                                P1[j] = SZ_UINT16_MIN;
446                        else
447                                P1[j] = SZ_UINT16_MAX;                 
448                }
449                else
450                {
451                        type[j] = 0;
452                        curValue = P1[j] = spaceFillingValue[j];
453                        compressUInt16Value(curValue, minValue, byteSize, bytes);
454                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
455                }
456        }
457
458        /* Process Row-1 --> Row-r1-1 */
459        size_t index;
460        for (i = 1; i < r1; i++)
461        {       
462                /* Process row-i data 0 */
463                index = i*r2;
464                pred1D = P1[0];
465                diff = spaceFillingValue[index] - pred1D;
466
467                itvNum = llabs(diff)/realPrecision + 1;
468
469                if (itvNum < exe_params->intvCapacity)
470                {
471                        if (diff < 0) itvNum = -itvNum;
472                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
473                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
474                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
475                                P0[0] = tmp;
476                        else if(tmp < SZ_UINT16_MIN)
477                                P0[0] = SZ_UINT16_MIN;
478                        else
479                                P0[0] = SZ_UINT16_MAX;                 
480                }
481                else
482                {
483                        type[index] = 0;
484                        curValue = P0[0] = spaceFillingValue[index];
485                        compressUInt16Value(curValue, minValue, byteSize, bytes);
486                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
487                }
488                                                                       
489                /* Process row-i data 1 --> r2-1*/
490                for (j = 1; j < r2; j++)
491                {
492                        index = i*r2+j;
493                        pred2D = P0[j-1] + P1[j] - P1[j-1];
494
495                        diff = spaceFillingValue[index] - pred2D;
496
497                        itvNum = llabs(diff)/realPrecision + 1;
498
499                        if (itvNum < exe_params->intvCapacity)
500                        {
501                                if (diff < 0) itvNum = -itvNum;
502                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
503                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
504                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
505                                        P0[j] = tmp;
506                                else if(tmp < SZ_UINT16_MIN)
507                                        P0[j] = SZ_UINT16_MIN;
508                                else
509                                        P0[j] = SZ_UINT16_MAX;                                         
510                        }
511                        else
512                        {
513                                type[index] = 0;
514                                curValue = P0[j] = spaceFillingValue[index];
515                                compressUInt16Value(curValue, minValue, byteSize, bytes);
516                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
517                        }
518                }
519
520                uint16_t *Pt;
521                Pt = P1;
522                P1 = P0;
523                P0 = Pt;
524        }
525       
526        if(r2!=1)
527                free(P0);
528        free(P1);                       
529       
530        size_t exactDataNum = exactDataByteArray->size;
531       
532        TightDataPointStorageI* tdps;   
533                       
534        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
535                        type, exactDataByteArray->array, exactDataByteArray->size, 
536                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
537                       
538        //free memory
539        free(type);     
540        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
541       
542        return tdps;   
543}
544
545/**
546 *
547 * Note: @r1 is high dimension
548 *               @r2 is low dimension
549 * */
550void SZ_compress_args_uint16_NoCkRngeNoGzip_2D(unsigned char** newByteData, uint16_t *oriData, size_t r1, size_t r2, double realPrecision, size_t *outSize, 
551int64_t valueRangeSize, uint16_t minValue)
552{
553        TightDataPointStorageI* tdps = SZ_compress_uint16_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, minValue);
554
555        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
556
557        size_t dataLength = r1*r2;
558        if(*outSize>dataLength*sizeof(uint16_t))
559                SZ_compress_args_uint16_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
560       
561        free_TightDataPointStorageI(tdps);     
562}
563
564TightDataPointStorageI* SZ_compress_uint16_3D_MDQ(uint16_t *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, int64_t valueRangeSize, int64_t minValue)
565{
566        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
567        int byteSize = computeByteSizePerIntValue(valueRangeSize);
568       
569        unsigned int quantization_intervals;
570        if(exe_params->optQuantMode==1)
571        {
572                quantization_intervals = optimize_intervals_uint16_3D(oriData, r1, r2, r3, realPrecision);
573                updateQuantizationInfo(quantization_intervals);
574        }       
575        else
576                quantization_intervals = exe_params->intvCapacity;
577        size_t i,j,k; 
578        int64_t pred1D, pred2D, pred3D, curValue, tmp;
579        int diff = 0.0;
580        double itvNum = 0;
581        uint16_t *P0, *P1;
582               
583        size_t dataLength = r1*r2*r3;           
584
585        size_t r23 = r2*r3;
586        P0 = (uint16_t*)malloc(r23*sizeof(uint16_t));
587        P1 = (uint16_t*)malloc(r23*sizeof(uint16_t));
588
589        int* type = (int*) malloc(dataLength*sizeof(int));
590
591        uint16_t* spaceFillingValue = oriData; //
592       
593        DynamicByteArray *exactDataByteArray;
594        new_DBA(&exactDataByteArray, DynArrayInitLen); 
595
596        type[0] = 0;
597        P1[0] = spaceFillingValue[0];
598        compressUInt16Value(spaceFillingValue[0], minValue, byteSize, bytes);
599        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
600
601        /* Process Row-0 data 1*/
602        pred1D = P1[0];
603        diff = spaceFillingValue[1] - pred1D;
604
605        itvNum = llabs(diff)/realPrecision + 1;
606
607        if (itvNum < exe_params->intvCapacity)
608        {
609                if (diff < 0) itvNum = -itvNum;
610                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
611                tmp = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
612                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
613                        P1[1] = tmp;
614                else if(tmp < SZ_UINT16_MIN)
615                        P1[1] = SZ_UINT16_MIN;
616                else
617                        P1[1] = SZ_UINT16_MAX;         
618        }
619        else
620        {
621                type[1] = 0;
622                curValue = P1[1] = spaceFillingValue[1];
623                compressUInt16Value(curValue, minValue, byteSize, bytes);
624                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
625        }
626
627    /* Process Row-0 data 2 --> data r3-1 */
628        for (j = 2; j < r3; j++)
629        {
630                pred1D = 2*P1[j-1] - P1[j-2];
631                diff = spaceFillingValue[j] - pred1D;
632
633                itvNum = llabs(diff)/realPrecision + 1;
634
635                if (itvNum < exe_params->intvCapacity)
636                {
637                        if (diff < 0) itvNum = -itvNum;
638                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
639                        tmp = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
640                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
641                                P1[j] = tmp;
642                        else if(tmp < SZ_UINT16_MIN)
643                                P1[j] = SZ_UINT16_MIN;
644                        else
645                                P1[j] = SZ_UINT16_MAX;                 
646                }
647                else
648                {
649                        type[j] = 0;
650                        curValue = P1[j] = spaceFillingValue[j];
651                        compressUInt16Value(curValue, minValue, byteSize, bytes);
652                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
653                }
654        }
655
656        /* Process Row-1 --> Row-r2-1 */
657        size_t index;
658        for (i = 1; i < r2; i++)
659        {
660                /* Process row-i data 0 */
661                index = i*r3;   
662                pred1D = P1[index-r3];
663                diff = spaceFillingValue[index] - pred1D;
664
665                itvNum = llabs(diff)/realPrecision + 1;
666
667                if (itvNum < exe_params->intvCapacity)
668                {
669                        if (diff < 0) itvNum = -itvNum;
670                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
671                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
672                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
673                                P1[index] = tmp;
674                        else if(tmp < SZ_UINT16_MIN)
675                                P1[index] = SZ_UINT16_MIN;
676                        else
677                                P1[index] = SZ_UINT16_MAX;                     
678                }
679                else
680                {
681                        type[index] = 0;
682                        curValue = P1[index] = spaceFillingValue[index];
683                        compressUInt16Value(curValue, minValue, byteSize, bytes);
684                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
685                }
686
687                /* Process row-i data 1 --> data r3-1*/
688                for (j = 1; j < r3; j++)
689                {
690                        index = i*r3+j;
691                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];
692
693                        diff = spaceFillingValue[index] - pred2D;
694
695                        itvNum = llabs(diff)/realPrecision + 1;
696
697                        if (itvNum < exe_params->intvCapacity)
698                        {
699                                if (diff < 0) itvNum = -itvNum;
700                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
701                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
702                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
703                                        P1[index] = tmp;
704                                else if(tmp < SZ_UINT16_MIN)
705                                        P1[index] = SZ_UINT16_MIN;
706                                else
707                                        P1[index] = SZ_UINT16_MAX;                             
708                        }
709                        else
710                        {
711                                type[index] = 0;
712                                curValue = P1[index] = spaceFillingValue[index];
713                                compressUInt16Value(curValue, minValue, byteSize, bytes);
714                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
715                        }
716                }
717        }
718
719
720        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////
721
722        for (k = 1; k < r1; k++)
723        {
724                /* Process Row-0 data 0*/
725                index = k*r23;
726                pred1D = P1[0];
727                diff = spaceFillingValue[index] - pred1D;
728
729                itvNum = llabs(diff)/realPrecision + 1;
730
731                if (itvNum < exe_params->intvCapacity)
732                {
733                        if (diff < 0) itvNum = -itvNum;
734                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
735                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
736                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
737                                P0[0] = tmp;
738                        else if(tmp < SZ_UINT16_MIN)
739                                P0[0] = SZ_UINT16_MIN;
740                        else
741                                P0[0] = SZ_UINT16_MAX;
742                }
743                else
744                {
745                        type[index] = 0;
746                        curValue = P0[0] = spaceFillingValue[index];
747                        compressUInt16Value(curValue, minValue, byteSize, bytes);
748                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
749                }
750
751
752            /* Process Row-0 data 1 --> data r3-1 */
753                for (j = 1; j < r3; j++)
754                {
755                        //index = k*r2*r3+j;
756                        index ++;
757                        pred2D = P0[j-1] + P1[j] - P1[j-1];
758                        diff = spaceFillingValue[index] - pred2D;
759
760                        itvNum = llabs(diff)/realPrecision + 1;
761
762                        if (itvNum < exe_params->intvCapacity)
763                        {
764                                if (diff < 0) itvNum = -itvNum;
765                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
766                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
767                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
768                                        P0[j] = tmp;
769                                else if(tmp < SZ_UINT16_MIN)
770                                        P0[j] = SZ_UINT16_MIN;
771                                else
772                                        P0[j] = SZ_UINT16_MAX;                         
773                        }
774                        else
775                        {
776                                type[index] = 0;
777                                curValue = P0[j] = spaceFillingValue[index];
778                                compressUInt16Value(curValue, minValue, byteSize, bytes);
779                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
780                        }
781                }
782
783            /* Process Row-1 --> Row-r2-1 */
784                size_t index2D;
785                for (i = 1; i < r2; i++)
786                {
787                        /* Process Row-i data 0 */
788                        index = k*r23 + i*r3;
789                        index2D = i*r3;         
790                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
791                        diff = spaceFillingValue[index] - pred2D;
792
793                        itvNum = llabs(diff)/realPrecision + 1;
794
795                        if (itvNum < exe_params->intvCapacity)
796                        {
797                                if (diff < 0) itvNum = -itvNum;
798                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
799                                tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
800                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
801                                        P0[index2D] = tmp;
802                                else if(tmp < SZ_UINT16_MIN)
803                                        P0[index2D] = SZ_UINT16_MIN;
804                                else
805                                        P0[index2D] = SZ_UINT16_MAX;
806                        }
807                        else
808                        {
809                                type[index] = 0;
810                                curValue = P0[index2D] = spaceFillingValue[index];
811                                compressUInt16Value(curValue, minValue, byteSize, bytes);
812                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
813                        }
814
815                        /* Process Row-i data 1 --> data r3-1 */
816                        for (j = 1; j < r3; j++)
817                        {
818//                              if(k==63&&i==43&&j==27)
819//                                      printf("i=%d\n", i);
820                                //index = k*r2*r3 + i*r3 + j;                   
821                                index ++;
822                                index2D = i*r3 + j;
823                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
824                                diff = spaceFillingValue[index] - pred3D;
825
826                                itvNum = llabs(diff)/realPrecision + 1;
827
828                                if (itvNum < exe_params->intvCapacity)
829                                {
830                                        if (diff < 0) itvNum = -itvNum;
831                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
832                                        tmp = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
833                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
834                                                P0[index2D] = tmp;
835                                        else if(tmp < SZ_UINT16_MIN)
836                                                P0[index2D] = SZ_UINT16_MIN;
837                                        else
838                                                P0[index2D] = SZ_UINT16_MAX;
839                                }
840                                else
841                                {
842                                        type[index] = 0;
843                                        curValue = P0[index2D] = spaceFillingValue[index];
844                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
845                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
846                                }
847                        }
848                }
849
850                uint16_t *Pt;
851                Pt = P1;
852                P1 = P0;
853                P0 = Pt;
854        }
855        if(r23!=1)
856                free(P0);
857        free(P1);
858
859        size_t exactDataNum = exactDataByteArray->size;
860       
861        TightDataPointStorageI* tdps;   
862                       
863        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
864                        type, exactDataByteArray->array, exactDataByteArray->size, 
865                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
866                       
867        //free memory
868        free(type);     
869        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
870       
871        return tdps;   
872}
873
874
875void SZ_compress_args_uint16_NoCkRngeNoGzip_3D(unsigned char** newByteData, uint16_t *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, size_t *outSize, 
876int64_t valueRangeSize, int64_t minValue)
877{       
878        TightDataPointStorageI* tdps = SZ_compress_uint16_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, minValue);
879
880        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
881
882        size_t dataLength = r1*r2*r3;
883        if(*outSize>dataLength*sizeof(uint16_t))
884                SZ_compress_args_uint16_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
885       
886        free_TightDataPointStorageI(tdps);     
887}
888
889
890TightDataPointStorageI* SZ_compress_uint16_4D_MDQ(uint16_t *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, int64_t valueRangeSize, int64_t minValue)
891{
892        unsigned char bytes[8] = {0,0,0,0,0,0,0,0};
893        int byteSize = computeByteSizePerIntValue(valueRangeSize);
894       
895        unsigned int quantization_intervals;
896        if(exe_params->optQuantMode==1)
897        {
898                quantization_intervals = optimize_intervals_uint16_4D(oriData, r1, r2, r3, r4, realPrecision);
899                updateQuantizationInfo(quantization_intervals);
900        }       
901        else
902                quantization_intervals = exe_params->intvCapacity;
903        size_t i,j,k; 
904        int64_t pred1D, pred2D, pred3D, curValue, tmp;
905        int diff = 0.0;
906        double itvNum = 0;
907        uint16_t *P0, *P1;
908               
909        size_t dataLength = r1*r2*r3*r4;               
910
911        size_t r234 = r2*r3*r4;
912        size_t r34 = r3*r4;
913
914        P0 = (uint16_t*)malloc(r34*sizeof(uint16_t));
915        P1 = (uint16_t*)malloc(r34*sizeof(uint16_t));
916       
917        int* type = (int*) malloc(dataLength*sizeof(int));
918
919        uint16_t* spaceFillingValue = oriData; //
920       
921        DynamicByteArray *exactDataByteArray;
922        new_DBA(&exactDataByteArray, DynArrayInitLen); 
923
924        size_t l;
925        for (l = 0; l < r1; l++)
926        {
927
928                ///////////////////////////     Process layer-0 ///////////////////////////
929                /* Process Row-0 data 0*/
930                size_t index = l*r234;
931                size_t index2D = 0;
932
933                type[index] = 0;
934                curValue = P1[index2D] = spaceFillingValue[index];
935                compressUInt16Value(curValue, minValue, byteSize, bytes);
936                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
937
938                /* Process Row-0 data 1*/
939                index = l*r234+1;
940                index2D = 1;
941
942                pred1D = P1[index2D-1];
943                diff = curValue - pred1D;
944
945                itvNum = llabs(diff)/realPrecision + 1;
946
947                if (itvNum < exe_params->intvCapacity)
948                {
949                        if (diff < 0) itvNum = -itvNum;
950                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
951                        tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
952                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
953                                P1[index2D] = tmp;
954                        else if(tmp < SZ_UINT16_MIN)
955                                P1[index2D] = SZ_UINT16_MIN;
956                        else
957                                P1[index2D] = SZ_UINT16_MAX;                   
958                }
959                else
960                {
961                        type[index] = 0;
962
963                        curValue = P1[index2D] = spaceFillingValue[0];
964                        compressUInt16Value(curValue, minValue, byteSize, bytes);
965                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
966                }
967
968                /* Process Row-0 data 2 --> data r4-1 */
969                for (j = 2; j < r4; j++)
970                {
971                        index = l*r234+j;
972                        index2D = j;
973
974                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
975                        diff = spaceFillingValue[index] - pred1D;
976
977                        itvNum = llabs(diff)/realPrecision + 1;
978
979                        if (itvNum < exe_params->intvCapacity)
980                        {
981                                if (diff < 0) itvNum = -itvNum;
982                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
983                                tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
984                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
985                                        P1[index2D] = tmp;
986                                else if(tmp < SZ_UINT16_MIN)
987                                        P1[index2D] = SZ_UINT16_MIN;
988                                else
989                                        P1[index2D] = SZ_UINT16_MAX;                                   
990                        }
991                        else
992                        {
993                                type[index] = 0;
994
995                                curValue = P1[index2D] = spaceFillingValue[0];
996                                compressUInt16Value(curValue, minValue, byteSize, bytes);
997                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
998                        }
999                }
1000
1001                /* Process Row-1 --> Row-r3-1 */
1002                for (i = 1; i < r3; i++)
1003                {
1004                        /* Process row-i data 0 */
1005                        index = l*r234+i*r4;
1006                        index2D = i*r4;
1007
1008                        pred1D = P1[index2D-r4];
1009                        diff = spaceFillingValue[index] - pred1D;
1010
1011                        itvNum = llabs(diff)/realPrecision + 1;
1012
1013                        if (itvNum < exe_params->intvCapacity)
1014                        {
1015                                if (diff < 0) itvNum = -itvNum;
1016                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1017                                tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1018                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1019                                        P1[index2D] = tmp;
1020                                else if(tmp < SZ_UINT16_MIN)
1021                                        P1[index2D] = SZ_UINT16_MIN;
1022                                else
1023                                        P1[index2D] = SZ_UINT16_MAX;                                   
1024                        }
1025                        else
1026                        {
1027                                type[index] = 0;
1028
1029                                curValue = P1[index2D] = spaceFillingValue[0];
1030                                compressUInt16Value(curValue, minValue, byteSize, bytes);
1031                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1032                        }
1033
1034                        /* Process row-i data 1 --> data r4-1*/
1035                        for (j = 1; j < r4; j++)
1036                        {
1037                                index = l*r234+i*r4+j;
1038                                index2D = i*r4+j;
1039
1040                                pred2D = P1[index2D-1] + P1[index2D-r4] - P1[index2D-r4-1];
1041
1042                                diff = spaceFillingValue[index] - pred2D;
1043
1044                                itvNum = llabs(diff)/realPrecision + 1;
1045
1046                                if (itvNum < exe_params->intvCapacity)
1047                                {
1048                                        if (diff < 0) itvNum = -itvNum;
1049                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1050                                        tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1051                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1052                                                P1[index2D] = tmp;
1053                                        else if(tmp < SZ_UINT16_MIN)
1054                                                P1[index2D] = SZ_UINT16_MIN;
1055                                        else
1056                                                P1[index2D] = SZ_UINT16_MAX;                                           
1057                                }
1058                                else
1059                                {
1060                                        type[index] = 0;
1061
1062                                        curValue = P1[index2D] = spaceFillingValue[0];
1063                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
1064                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1065                                }
1066                        }
1067                }
1068
1069
1070                ///////////////////////////     Process layer-1 --> layer-r2-1 ///////////////////////////
1071
1072                for (k = 1; k < r2; k++)
1073                {
1074                        /* Process Row-0 data 0*/
1075                        index = l*r234+k*r34;
1076                        index2D = 0;
1077
1078                        pred1D = P1[index2D];
1079                        diff = spaceFillingValue[index] - pred1D;
1080
1081                        itvNum = llabs(diff)/realPrecision + 1;
1082
1083                        if (itvNum < exe_params->intvCapacity)
1084                        {
1085                                if (diff < 0) itvNum = -itvNum;
1086                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1087                                tmp = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1088                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1089                                        P0[index2D] = tmp;
1090                                else if(tmp < SZ_UINT16_MIN)
1091                                        P0[index2D] = SZ_UINT16_MIN;
1092                                else
1093                                        P0[index2D] = SZ_UINT16_MAX;                                   
1094                        }
1095                        else
1096                        {
1097                                type[index] = 0;
1098
1099                                curValue = P0[index2D] = spaceFillingValue[0];
1100                                compressUInt16Value(curValue, minValue, byteSize, bytes);
1101                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1102                        }
1103
1104                        /* Process Row-0 data 1 --> data r4-1 */
1105                        for (j = 1; j < r4; j++)
1106                        {
1107                                index = l*r234+k*r34+j;
1108                                index2D = j;
1109
1110                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
1111                                diff = spaceFillingValue[index] - pred2D;
1112
1113                                itvNum = llabs(diff)/realPrecision + 1;
1114
1115                                if (itvNum < exe_params->intvCapacity)
1116                                {
1117                                        if (diff < 0) itvNum = -itvNum;
1118                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1119                                        tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1120                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1121                                                P0[index2D] = tmp;
1122                                        else if(tmp < SZ_UINT16_MIN)
1123                                                P0[index2D] = SZ_UINT16_MIN;
1124                                        else
1125                                                P0[index2D] = SZ_UINT16_MAX;                                           
1126                                }
1127                                else
1128                                {
1129                                        type[index] = 0;
1130
1131                                        curValue = P0[index2D] = spaceFillingValue[0];
1132                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
1133                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1134                                }
1135                        }
1136
1137                        /* Process Row-1 --> Row-r3-1 */
1138                        for (i = 1; i < r3; i++)
1139                        {
1140                                /* Process Row-i data 0 */
1141                                index = l*r234+k*r34+i*r4;
1142                                index2D = i*r4;
1143
1144                                pred2D = P0[index2D-r4] + P1[index2D] - P1[index2D-r4];
1145                                diff = spaceFillingValue[index] - pred2D;
1146
1147                                itvNum = llabs(diff)/realPrecision + 1;
1148
1149                                if (itvNum < exe_params->intvCapacity)
1150                                {
1151                                        if (diff < 0) itvNum = -itvNum;
1152                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1153                                        tmp = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1154                                        if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1155                                                P0[index2D] = tmp;
1156                                        else if(tmp < SZ_UINT16_MIN)
1157                                                P0[index2D] = SZ_UINT16_MIN;
1158                                        else
1159                                                P0[index2D] = SZ_UINT16_MAX;                                           
1160                                }
1161                                else
1162                                {
1163                                        type[index] = 0;
1164
1165                                        curValue = P0[index2D] = spaceFillingValue[0];
1166                                        compressUInt16Value(curValue, minValue, byteSize, bytes);
1167                                        memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1168                                }
1169
1170                                /* Process Row-i data 1 --> data r4-1 */
1171                                for (j = 1; j < r4; j++)
1172                                {
1173                                        index = l*r234+k*r34+i*r4+j;
1174                                        index2D = i*r4+j;
1175
1176                                        pred3D = P0[index2D-1] + P0[index2D-r4]+ P1[index2D] - P0[index2D-r4-1] - P1[index2D-r4] - P1[index2D-1] + P1[index2D-r4-1];
1177                                        diff = spaceFillingValue[index] - pred3D;
1178
1179
1180                                        itvNum = llabs(diff)/realPrecision + 1;
1181
1182                                        if (itvNum < exe_params->intvCapacity)
1183                                        {
1184                                                if (diff < 0) itvNum = -itvNum;
1185                                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
1186                                                tmp = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
1187                                                if(tmp >= SZ_UINT16_MIN&&tmp<SZ_UINT16_MAX)
1188                                                        P0[index2D] = tmp;
1189                                                else if(tmp < SZ_UINT16_MIN)
1190                                                        P0[index2D] = SZ_UINT16_MIN;
1191                                                else
1192                                                        P0[index2D] = SZ_UINT16_MAX;                                                   
1193                                        }
1194                                        else
1195                                        {
1196                                                type[index] = 0;
1197
1198                                                curValue = P0[index2D] = spaceFillingValue[0];
1199                                                compressUInt16Value(curValue, minValue, byteSize, bytes);
1200                                                memcpyDBA_Data(exactDataByteArray, bytes, byteSize);
1201                                        }
1202                                }
1203                        }
1204
1205                        uint16_t *Pt;
1206                        Pt = P1;
1207                        P1 = P0;
1208                        P0 = Pt;
1209                }
1210        }
1211
1212        free(P0);
1213        free(P1);
1214
1215        size_t exactDataNum = exactDataByteArray->size;
1216       
1217        TightDataPointStorageI* tdps;   
1218                       
1219        new_TightDataPointStorageI(&tdps, dataLength, exactDataNum, byteSize, 
1220                        type, exactDataByteArray->array, exactDataByteArray->size, 
1221                        realPrecision, minValue, quantization_intervals, SZ_UINT16);
1222                       
1223        //free memory
1224        free(type);     
1225        free(exactDataByteArray); //exactDataByteArray->array has been released in free_TightDataPointStorageF(tdps);
1226       
1227        return tdps;   
1228}
1229
1230void SZ_compress_args_uint16_NoCkRngeNoGzip_4D(unsigned char** newByteData, uint16_t *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, 
1231size_t *outSize, int64_t valueRangeSize, int64_t minValue)
1232{
1233        TightDataPointStorageI* tdps = SZ_compress_uint16_4D_MDQ(oriData, r1, r2, r3, r4, realPrecision, valueRangeSize, minValue);
1234
1235        convertTDPStoFlatBytes_int(tdps, newByteData, outSize);
1236
1237        size_t dataLength = r1*r2*r3*r4;
1238        if(*outSize>dataLength*sizeof(uint16_t))
1239                SZ_compress_args_uint16_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
1240
1241        free_TightDataPointStorageI(tdps);
1242}
1243
1244void SZ_compress_args_uint16_withinRange(unsigned char** newByteData, uint16_t *oriData, size_t dataLength, size_t *outSize)
1245{
1246        TightDataPointStorageI* tdps = (TightDataPointStorageI*) malloc(sizeof(TightDataPointStorageI));
1247        tdps->typeArray = NULL; 
1248       
1249        tdps->allSameData = 1;
1250        tdps->dataSeriesLength = dataLength;
1251        tdps->exactDataBytes = (unsigned char*)malloc(sizeof(unsigned char)*2);
1252        tdps->isLossless = 0;
1253        //tdps->exactByteSize = 4;
1254        tdps->exactDataNum = 1;
1255        tdps->exactDataBytes_size = 2;
1256       
1257        uint16_t value = oriData[0];
1258        int16ToBytes_bigEndian(tdps->exactDataBytes, value);
1259       
1260        size_t tmpOutSize;
1261        convertTDPStoFlatBytes_int(tdps, newByteData, &tmpOutSize);
1262
1263        *outSize = tmpOutSize;//3+1+sizeof(uint16_t)+SZ_SIZE_TYPE; //8==3+1+4(uint16_size)
1264        free_TightDataPointStorageI(tdps);     
1265}
1266
1267int SZ_compress_args_uint16_wRngeNoGzip(unsigned char** newByteData, uint16_t *oriData, 
1268size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1269int errBoundMode, double absErr_Bound, double relBoundRatio)
1270{
1271        int status = SZ_SCES;
1272        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1273        int64_t valueRangeSize = 0;
1274       
1275        uint16_t minValue = computeRangeSize_int(oriData, SZ_UINT16, dataLength, &valueRangeSize);
1276        double realPrecision = getRealPrecision_int(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1277               
1278        if(valueRangeSize <= realPrecision)
1279        {
1280                SZ_compress_args_uint16_withinRange(newByteData, oriData, dataLength, outSize);
1281        }
1282        else
1283        {
1284//              SZ_compress_args_uint16_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize);
1285                if(r5==0&&r4==0&&r3==0&&r2==0)
1286                {
1287                        SZ_compress_args_uint16_NoCkRngeNoGzip_1D(newByteData, oriData, r1, realPrecision, outSize, valueRangeSize, minValue);
1288                }
1289                else if(r5==0&&r4==0&&r3==0)
1290                {
1291                        SZ_compress_args_uint16_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize, valueRangeSize, minValue);
1292                }
1293                else if(r5==0&&r4==0)
1294                {
1295                        SZ_compress_args_uint16_NoCkRngeNoGzip_3D(newByteData, oriData, r3, r2, r1, realPrecision, outSize, valueRangeSize, minValue);
1296                }
1297                else if(r5==0)
1298                {
1299                        SZ_compress_args_uint16_NoCkRngeNoGzip_3D(newByteData, oriData, r4*r3, r2, r1, realPrecision, outSize, valueRangeSize, minValue);
1300                }
1301        }
1302        return status;
1303}
1304
1305int SZ_compress_args_uint16(unsigned char** newByteData, uint16_t *oriData, 
1306size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
1307int errBoundMode, double absErr_Bound, double relBoundRatio)
1308{
1309        confparams_cpr->errorBoundMode = errBoundMode;
1310       
1311        if(errBoundMode>=PW_REL)
1312        {
1313                printf("Error: Current SZ version doesn't support integer data compression with point-wise relative error bound being based on pwrType=AVG\n");
1314                exit(0);
1315                return SZ_NSCS;
1316        }
1317        int status = SZ_SCES;
1318        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
1319        int64_t valueRangeSize = 0;
1320
1321        uint16_t minValue = (uint16_t)computeRangeSize_int(oriData, SZ_UINT16, dataLength, &valueRangeSize);
1322        double realPrecision = 0; 
1323       
1324        if(confparams_cpr->errorBoundMode==PSNR)
1325        {
1326                confparams_cpr->errorBoundMode = ABS;
1327                realPrecision = confparams_cpr->absErrBound = computeABSErrBoundFromPSNR(confparams_cpr->psnr, (double)confparams_cpr->predThreshold, (double)valueRangeSize);
1328                //printf("realPrecision=%lf\n", realPrecision);
1329        }
1330        else
1331                realPrecision = getRealPrecision_int(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
1332
1333        if(valueRangeSize <= realPrecision)
1334        {
1335                SZ_compress_args_uint16_withinRange(newByteData, oriData, dataLength, outSize);
1336        }
1337        else
1338        {
1339                size_t tmpOutSize = 0;
1340                unsigned char* tmpByteData;
1341                if (r2==0)
1342                {
1343                        SZ_compress_args_uint16_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1344                }
1345                else
1346                if (r3==0)
1347                {
1348                        SZ_compress_args_uint16_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1349                }
1350                else
1351                if (r4==0)
1352                {
1353                        SZ_compress_args_uint16_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1354                }
1355                else
1356                if (r5==0)
1357                {
1358                        SZ_compress_args_uint16_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, minValue);
1359                }
1360                else
1361                {
1362                        printf("Error: doesn't support 5 dimensions for now.\n");
1363                        status = SZ_DERR; //dimension error
1364                }
1365                //Call Gzip to do the further compression.
1366                if(confparams_cpr->szMode==SZ_BEST_SPEED)
1367                {
1368                        *outSize = tmpOutSize;
1369                        *newByteData = tmpByteData;
1370                }
1371                else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
1372                {
[9ee2ce3]1373                        *outSize = sz_lossless_compress(confparams_cpr->losslessCompressor, confparams_cpr->gzipMode, tmpByteData, tmpOutSize, newByteData);
[2c47b73]1374                        free(tmpByteData);
1375                }
1376                else
1377                {
1378                        printf("Error: Wrong setting of confparams_cpr->szMode in the uint16_t compression.\n");
1379                        status = SZ_MERR; //mode error                 
1380                }
1381        }
1382       
1383        return status;
1384}
Note: See TracBrowser for help on using the repository browser.