source: thirdparty/SZ/sz/src/sz_double.c @ 2c47b73

/**
 *  @file sz_double.c
 *  @author Sheng Di and Dingwen Tao
 *  @date Aug, 2016
 *  @brief SZ_Init, Compression and Decompression functions
 *  (C) 2016 by Mathematics and Computer Science (MCS), Argonne National Laboratory.
 *      See COPYRIGHT in top-level directory.
 */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <math.h>
#include "sz.h"
#include "CompressElement.h"
#include "DynamicByteArray.h"
#include "DynamicIntArray.h"
#include "TightDataPointStorageD.h"
#include "sz_double.h"
#include "sz_double_pwr.h"
#include "szd_double.h"
#include "szd_double_pwr.h"
#include "zlib.h"
#include "rw.h"
#include "sz_double_ts.h"

unsigned char* SZ_skip_compress_double(double* data, size_t dataLength, size_t* outSize)
{
        *outSize = dataLength*sizeof(double);
        unsigned char* out = (unsigned char*)malloc(dataLength*sizeof(double));
        memcpy(out, data, dataLength*sizeof(double));
        return out;
}

void computeReqLength_double(double realPrecision, short radExpo, int* reqLength, double* medianValue)
{
        short reqExpo = getPrecisionReqLength_double(realPrecision);
        *reqLength = 12+radExpo - reqExpo; //radExpo-reqExpo == reqMantiLength
        if(*reqLength<12)
                *reqLength = 12;
        if(*reqLength>64)
        {
                *reqLength = 64;
                *medianValue = 0;
        }
}

unsigned int optimize_intervals_double_1D(double *oriData, size_t dataLength, double realPrecision)
{       
        size_t i = 0, radiusIndex;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
        for(i=2;i<dataLength;i++)
        {
                if(i%confparams_cpr->sampleDistance==0)
                {
                        //pred_value = 2*oriData[i-1] - oriData[i-2];
                        pred_value = oriData[i-1];
                        pred_err = fabs(pred_value - oriData[i]);
                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
                        intervals[radiusIndex]++;
                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }

        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
        return powerOf2;
}

unsigned int optimize_intervals_double_2D(double *oriData, size_t r1, size_t r2, double realPrecision)
{       
        size_t i,j, index;
        size_t radiusIndex;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = (r1-1)*(r2-1)/confparams_cpr->sampleDistance;
        for(i=1;i<r1;i++)
        {
                for(j=1;j<r2;j++)
                {
                        if((i+j)%confparams_cpr->sampleDistance==0)
                        {
                                index = i*r2+j;
                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
                                pred_err = fabs(pred_value - oriData[index]);
                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
                                intervals[radiusIndex]++;
                        }                       
                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;   
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

unsigned int optimize_intervals_double_3D(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision)
{       
        size_t i,j,k, index;
        size_t radiusIndex;
        size_t r23=r2*r3;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)/confparams_cpr->sampleDistance;
        for(i=1;i<r1;i++)
        {
                for(j=1;j<r2;j++)
                {
                        for(k=1;k<r3;k++)
                        {
                                if((i+j+k)%confparams_cpr->sampleDistance==0)
                                {
                                        index = i*r23+j*r3+k;
                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23] 
                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
                                        pred_err = fabs(pred_value - oriData[index]);
                                        radiusIndex = (pred_err/realPrecision+1)/2;
                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
                                        intervals[radiusIndex]++;
                                }                               
                        }
                        
                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
                
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        //printf("confparams_cpr->maxRangeRadius = %d, accIntervals=%d, powerOf2=%d\n", confparams_cpr->maxRangeRadius, accIntervals, powerOf2);
        return powerOf2;
}

unsigned int optimize_intervals_double_4D(double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision)
{
        size_t i,j,k,l, index;
        size_t radiusIndex;
        size_t r234=r2*r3*r4;
        size_t r34=r3*r4;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = (r1-1)*(r2-1)*(r3-1)*(r4-1)/confparams_cpr->sampleDistance;
        for(i=1;i<r1;i++)
        {
                for(j=1;j<r2;j++)
                {
                        for(k=1;k<r3;k++)
                        {
                                for (l=1;l<r4;l++)
                                {
                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
                                        {
                                                index = i*r234+j*r34+k*r4+l;
                                                pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r34]
                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
                                                pred_err = fabs(pred_value - oriData[index]);
                                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
                                                intervals[radiusIndex]++;
                                        }
                                }
                        }
                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;

        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

TightDataPointStorageD* SZ_compress_double_1D_MDQ(double *oriData, 
size_t dataLength, double realPrecision, double valueRangeSize, double medianValue_d)
{
#ifdef HAVE_TIMECMPR
        double* decData = NULL; 
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData = (double*)(multisteps->hist_data);
#endif  
        
        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
                quantization_intervals = optimize_intervals_double_1D_opt(oriData, dataLength, realPrecision);
        else
                quantization_intervals = exe_params->intvCapacity;
        updateQuantizationInfo(quantization_intervals); 

        size_t i;
        int reqLength;
        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);

        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);      

        int* type = (int*) malloc(dataLength*sizeof(int));
                
        double* spaceFillingValue = oriData; //
        
        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);
        
        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);
        
        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);
        
        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;
        double last3CmprsData[3] = {0};

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));                       
                                
        //add the first data    
        type[0] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        listAdd_double(last3CmprsData, vce->data);
#ifdef HAVE_TIMECMPR    
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData[0] = vce->data;
#endif          
                
        //add the second data
        type[1] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        listAdd_double(last3CmprsData, vce->data);
#ifdef HAVE_TIMECMPR    
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData[1] = vce->data;
#endif
        int state;
        double checkRadius;
        double curData;
        double pred;
        double predAbsErr;
        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
        double interval = 2*realPrecision;

        for(i=2;i<dataLength;i++)
        {                               
                //printf("%.30G\n",last3CmprsData[0]);
                curData = spaceFillingValue[i];
                //pred = 2*last3CmprsData[0] - last3CmprsData[1];
                pred = last3CmprsData[0];
                predAbsErr = fabs(curData - pred);      
                if(predAbsErr<=checkRadius)
                {
                        state = (predAbsErr/realPrecision+1)/2;
                        if(curData>=pred)
                        {
                                type[i] = exe_params->intvRadius+state;
                                pred = pred + state*interval;
                        }
                        else //curData<pred
                        {
                                type[i] = exe_params->intvRadius-state;
                                pred = pred - state*interval;
                        }
                        listAdd_double(last3CmprsData, pred);
#ifdef HAVE_TIMECMPR                                    
                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                decData[i] = pred;                      
#endif  
                        continue;
                }
                
                //unpredictable data processing
                type[i] = 0;            
                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                                        
                listAdd_double(last3CmprsData, vce->data);
#ifdef HAVE_TIMECMPR
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[i] = vce->data;
#endif  
                
        }//end of for
                
        int exactDataNum = exactLeadNumArray->size;
        
        TightDataPointStorageD* tdps;
                        
        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum, 
                        type, exactMidByteArray->array, exactMidByteArray->size,  
                        exactLeadNumArray->array,  
                        resiBitArray->array, resiBitArray->size, 
                        resiBitsLength, 
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);
        
//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n", 
//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
        
        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);      
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);     
        
        return tdps;    
}

void SZ_compress_args_double_StoreOriData(double* oriData, size_t dataLength, TightDataPointStorageD* tdps, 
unsigned char** newByteData, size_t *outSize)
{
        int doubleSize = sizeof(double);
        size_t k = 0, i;
        tdps->isLossless = 1;
        size_t totalByteLength = 3 + MetaDataByteLength + exe_params->SZ_SIZE_TYPE + 1 + doubleSize*dataLength;
        *newByteData = (unsigned char*)malloc(totalByteLength);
        
        unsigned char dsLengthBytes[8];
        for (i = 0; i < 3; i++)//3
                (*newByteData)[k++] = versionNumber[i];
        
        if(exe_params->SZ_SIZE_TYPE==4)//1
                (*newByteData)[k++] = 16; //00010000
        else
                (*newByteData)[k++] = 80;       //01010000: 01000000 indicates the SZ_SIZE_TYPE=8

        convertSZParamsToBytes(confparams_cpr, &((*newByteData)[k]));
        k = k + MetaDataByteLength;

        sizeToBytes(dsLengthBytes,dataLength);
        for (i = 0; i < exe_params->SZ_SIZE_TYPE; i++)//ST: 4 or 8
                (*newByteData)[k++] = dsLengthBytes[i];

        if(sysEndianType==BIG_ENDIAN_SYSTEM)
                memcpy((*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE, oriData, dataLength*doubleSize);
        else
        {
                unsigned char* p = (*newByteData)+4+MetaDataByteLength+exe_params->SZ_SIZE_TYPE;
                for(i=0;i<dataLength;i++,p+=doubleSize)
                        doubleToBytes(p, oriData[i]);
        }
        *outSize = totalByteLength;
}


char SZ_compress_args_double_NoCkRngeNoGzip_1D(unsigned char** newByteData, double *oriData, 
size_t dataLength, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
{
        char compressionType = 0;       
        TightDataPointStorageD* tdps = NULL;    
#ifdef HAVE_TIMECMPR
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
        {
                int timestep = sz_tsc->currentStep;
                if(timestep % confparams_cpr->snapshotCmprStep != 0)
                {
                        tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
                        compressionType = 1; //time-series based compression 
                }
                else
                {       
                        tdps = SZ_compress_double_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_d);
                        compressionType = 0; //snapshot-based compression
                        multisteps->lastSnapshotStep = timestep;
                }               
        }
        else
#endif
                tdps = SZ_compress_double_1D_MDQ(oriData, dataLength, realPrecision, valueRangeSize, medianValue_d);                    
        
        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
        
        if(*outSize>dataLength*sizeof(double))
                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);
        
        free_TightDataPointStorageD(tdps);      
        return compressionType;
}

TightDataPointStorageD* SZ_compress_double_2D_MDQ(double *oriData, size_t r1, size_t r2, double realPrecision, double valueRangeSize, double medianValue_d)
{
#ifdef HAVE_TIMECMPR    
        double* decData = NULL;
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData = (double*)(multisteps->hist_data);
#endif  
        
        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
        {
                quantization_intervals = optimize_intervals_double_2D_opt(oriData, r1, r2, realPrecision);
                updateQuantizationInfo(quantization_intervals);
        }
        else
                quantization_intervals = exe_params->intvCapacity;      
        size_t i,j; 
        int reqLength;
        double pred1D, pred2D;
        double diff = 0.0;
        double itvNum = 0;
        double *P0, *P1;
                
        size_t dataLength = r1*r2;      
        
        P0 = (double*)malloc(r2*sizeof(double));
        memset(P0, 0, r2*sizeof(double));
        P1 = (double*)malloc(r2*sizeof(double));
        memset(P1, 0, r2*sizeof(double));
                
        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);
        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);      

        int* type = (int*) malloc(dataLength*sizeof(int));
        //type[dataLength]=0;
                
        double* spaceFillingValue = oriData; //
        
        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);
        
        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);
        
        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);
        
        type[0] = 0;
        
        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);
        
        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));
                        
        /* Process Row-0 data 0*/
        type[0] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        P1[0] = vce->data;
#ifdef HAVE_TIMECMPR    
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData[0] = vce->data;
#endif  

        /* Process Row-0 data 1*/
        pred1D = P1[0];
        diff = spaceFillingValue[1] - pred1D;

        itvNum =  fabs(diff)/realPrecision + 1;

        if (itvNum < exe_params->intvCapacity)
        {
                if (diff < 0) itvNum = -itvNum;
                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
        }
        else
        {
                type[1] = 0;
                compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                P1[1] = vce->data;
        }
#ifdef HAVE_TIMECMPR    
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData[1] = P1[1];
#endif

    /* Process Row-0 data 2 --> data r2-1 */
        for (j = 2; j < r2; j++)
        {
                pred1D = 2*P1[j-1] - P1[j-2];
                diff = spaceFillingValue[j] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[j] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[j] = vce->data;
                }
#ifdef HAVE_TIMECMPR    
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[j] = P1[j];
#endif          
        }

        /* Process Row-1 --> Row-r1-1 */
        size_t index;
        for (i = 1; i < r1; i++)
        {       
                /* Process row-i data 0 */
                index = i*r2;
                pred1D = P1[0];
                diff = spaceFillingValue[index] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[index] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P0[0] = vce->data;
                }
#ifdef HAVE_TIMECMPR    
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[index] = P0[0];
#endif
                                                                        
                /* Process row-i data 1 --> r2-1*/
                for (j = 1; j < r2; j++)
                {
                        index = i*r2+j;
                        pred2D = P0[j-1] + P1[j] - P1[j-1];

                        diff = spaceFillingValue[index] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[j] = vce->data;
                        }
#ifdef HAVE_TIMECMPR    
                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                decData[index] = P0[j];
#endif                  
                }

                double *Pt;
                Pt = P1;
                P1 = P0;
                P0 = Pt;
        }
                
        if(r2!=1)       
                free(P0);
        free(P1);
        size_t exactDataNum = exactLeadNumArray->size;
        
        TightDataPointStorageD* tdps;
                        
        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum, 
                        type, exactMidByteArray->array, exactMidByteArray->size,  
                        exactLeadNumArray->array,  
                        resiBitArray->array, resiBitArray->size, 
                        resiBitsLength, 
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

/*      int sum =0;
        for(i=0;i<dataLength;i++)
                if(type[i]==0) sum++;
        printf("opt_quantizations=%d, exactDataNum=%d, sum=%d\n",quantization_intervals, exactDataNum, sum);

        for(i=0;i<dataLength;i++)
                printf("%d ", type[i]);
        printf("\n");*/

//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n", 
//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);
        
//      for(i = 3800;i<3844;i++)
//              printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);
        
        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);     
        free(vce);
        free(lce);      
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);
        
        return tdps;
}

/**
 * 
 * Note: @r1 is high dimension
 *               @r2 is low dimension 
 * */
char SZ_compress_args_double_NoCkRngeNoGzip_2D(unsigned char** newByteData, double *oriData, size_t r1, size_t r2, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
{
        size_t dataLength = r1*r2;
        char compressionType = 0;       
        TightDataPointStorageD* tdps = NULL;    
#ifdef HAVE_TIMECMPR
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
        {
                int timestep = sz_tsc->currentStep;
                if(timestep % confparams_cpr->snapshotCmprStep != 0)
                {
                        tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
                        compressionType = 1; //time-series based compression 
                }
                else
                {       
                        tdps = SZ_compress_double_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_d);
                        compressionType = 0; //snapshot-based compression
                        multisteps->lastSnapshotStep = timestep;
                }               
        }
        else
#endif
                tdps = SZ_compress_double_2D_MDQ(oriData, r1, r2, realPrecision, valueRangeSize, medianValue_d);        
        
        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);
        
        if(*outSize>dataLength*sizeof(double))
                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);  
        
        free_TightDataPointStorageD(tdps);
        return compressionType;
}

TightDataPointStorageD* SZ_compress_double_3D_MDQ(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, double valueRangeSize, double medianValue_d)
{
#ifdef HAVE_TIMECMPR
        double* decData = NULL;
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData = (double*)(multisteps->hist_data);
#endif          

        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
        {
                quantization_intervals = optimize_intervals_double_3D_opt(oriData, r1, r2, r3, realPrecision);
                updateQuantizationInfo(quantization_intervals);
        }       
        else
                quantization_intervals = exe_params->intvCapacity;
        size_t i,j,k; 
        int reqLength;
        double pred1D, pred2D, pred3D;
        double diff = 0.0;
        double itvNum = 0;
        double *P0, *P1;

        size_t dataLength = r1*r2*r3;

        size_t r23 = r2*r3;

        P0 = (double*)malloc(r23*sizeof(double));
        P1 = (double*)malloc(r23*sizeof(double));

        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);
        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);      

        int* type = (int*) malloc(dataLength*sizeof(int));
        //type[dataLength]=0;

        double* spaceFillingValue = oriData; //

        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);

        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);

        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        type[0] = 0;

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);

        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));


        ///////////////////////////     Process layer-0 ///////////////////////////
        /* Process Row-0 data 0*/
        type[0] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        P1[0] = vce->data;
#ifdef HAVE_TIMECMPR    
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[0] = P1[0];
#endif

        /* Process Row-0 data 1*/
        pred1D = P1[0];
        diff = spaceFillingValue[1] - pred1D;

        itvNum = fabs(diff)/realPrecision + 1;

        if (itvNum < exe_params->intvCapacity)
        {
                if (diff < 0) itvNum = -itvNum;
                type[1] = (int) (itvNum/2) + exe_params->intvRadius;
                P1[1] = pred1D + 2 * (type[1] - exe_params->intvRadius) * realPrecision;
        }
        else
        {
                type[1] = 0;
                compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                P1[1] = vce->data;
        }
#ifdef HAVE_TIMECMPR    
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                decData[1] = P1[1];
#endif

    /* Process Row-0 data 2 --> data r3-1 */
        for (j = 2; j < r3; j++)
        {
                pred1D = 2*P1[j-1] - P1[j-2];
                diff = spaceFillingValue[j] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[j] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[j] = pred1D + 2 * (type[j] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[j] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[j], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[j] = vce->data;
                }
#ifdef HAVE_TIMECMPR    
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[j] = P1[j];
#endif          
        }

        /* Process Row-1 --> Row-r2-1 */
        size_t index;
        for (i = 1; i < r2; i++)
        {
                /* Process row-i data 0 */
                index = i*r3;
                pred1D = P1[index-r3];
                diff = spaceFillingValue[index] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[index] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[index] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[index] = vce->data;
                }
#ifdef HAVE_TIMECMPR    
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[index] = P1[index];
#endif          

                /* Process row-i data 1 --> data r3-1*/
                for (j = 1; j < r3; j++)
                {
                        index = i*r3+j;
                        pred2D = P1[index-1] + P1[index-r3] - P1[index-r3-1];

                        diff = spaceFillingValue[index] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P1[index] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P1[index] = vce->data;
                        }
#ifdef HAVE_TIMECMPR    
                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                decData[index] = P1[index];
#endif                  
                }
        }


        ///////////////////////////     Process layer-1 --> layer-r1-1 ///////////////////////////

        for (k = 1; k < r1; k++)
        {
                /* Process Row-0 data 0*/
                index = k*r23;
                pred1D = P1[0];
                diff = spaceFillingValue[index] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                        P0[0] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[index] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P0[0] = vce->data;
                }
#ifdef HAVE_TIMECMPR    
                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                        decData[index] = P0[0];
#endif

            /* Process Row-0 data 1 --> data r3-1 */
                for (j = 1; j < r3; j++)
                {
                        //index = k*r2*r3+j;
                        index ++;
                        pred2D = P0[j-1] + P1[j] - P1[j-1];
                        diff = spaceFillingValue[index] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[j] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[j] = vce->data;
                        }
#ifdef HAVE_TIMECMPR    
                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                decData[index] = P0[j];
#endif                  
                }

            /* Process Row-1 --> Row-r2-1 */
                size_t index2D;
                for (i = 1; i < r2; i++)
                {
                        /* Process Row-i data 0 */
                        index = k*r23 + i*r3;
                        index2D = i*r3;
                        pred2D = P0[index2D-r3] + P1[index2D] - P1[index2D-r3];
                        diff = spaceFillingValue[index] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[index2D] = vce->data;
                        }
#ifdef HAVE_TIMECMPR    
                        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                decData[index] = P0[index2D];
#endif                  

                        /* Process Row-i data 1 --> data r3-1 */
                        for (j = 1; j < r3; j++)
                        {
                                //index = k*r2*r3 + i*r3 + j;
                                index ++;
                                index2D = i*r3 + j;
                                pred3D = P0[index2D-1] + P0[index2D-r3]+ P1[index2D] - P0[index2D-r3-1] - P1[index2D-r3] - P1[index2D-1] + P1[index2D-r3-1];
                                diff = spaceFillingValue[index] - pred3D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[index] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P0[index2D] = vce->data;
                                }
#ifdef HAVE_TIMECMPR    
                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                        decData[index] = P0[index2D];
#endif                          
                        }
                }

                double *Pt;
                Pt = P1;
                P1 = P0;
                P0 = Pt;
        }
        if(r23!=1)
                free(P0);
        free(P1);
        size_t exactDataNum = exactLeadNumArray->size;

        TightDataPointStorageD* tdps;

        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
                        type, exactMidByteArray->array, exactMidByteArray->size,
                        exactLeadNumArray->array,
                        resiBitArray->array, resiBitArray->size,
                        resiBitsLength, 
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

//      printf("exactDataNum=%d, expSegmentsInBytes_size=%d, exactMidByteArray->size=%d\n",
//                      exactDataNum, expSegmentsInBytes_size, exactMidByteArray->size);

//      for(i = 3800;i<3844;i++)
//              printf("exactLeadNumArray->array[%d]=%d\n",i,exactLeadNumArray->array[i]);

        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);      
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);     
        
        return tdps;    
}


char SZ_compress_args_double_NoCkRngeNoGzip_3D(unsigned char** newByteData, double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
{
        size_t dataLength = r1*r2*r3;
        char compressionType = 0;       
        TightDataPointStorageD* tdps = NULL;    
#ifdef HAVE_TIMECMPR
        if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
        {
                int timestep = sz_tsc->currentStep;
                if(timestep % confparams_cpr->snapshotCmprStep != 0)
                {
                        tdps = SZ_compress_double_1D_MDQ_ts(oriData, dataLength, multisteps, realPrecision, valueRangeSize, medianValue_d);
                        compressionType = 1; //time-series based compression 
                }
                else
                {       
                        tdps = SZ_compress_double_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_d);
                        compressionType = 0; //snapshot-based compression
                        multisteps->lastSnapshotStep = timestep;
                }               
        }
        else
#endif
                tdps = SZ_compress_double_3D_MDQ(oriData, r1, r2, r3, realPrecision, valueRangeSize, medianValue_d);            

        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);

        if(*outSize>dataLength*sizeof(double))
                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);

        free_TightDataPointStorageD(tdps);
        return compressionType;
}

TightDataPointStorageD* SZ_compress_double_4D_MDQ(double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, double valueRangeSize, double medianValue_d)
{
        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
        {
                quantization_intervals = optimize_intervals_double_4D(oriData, r1, r2, r3, r4, realPrecision);
                updateQuantizationInfo(quantization_intervals);
        }
        else
                quantization_intervals = exe_params->intvCapacity;

        size_t i,j,k; 
        int reqLength;
        double pred1D, pred2D, pred3D;
        double diff = 0.0;
        double itvNum = 0;
        double *P0, *P1;

        size_t dataLength = r1*r2*r3*r4;

        size_t r234 = r2*r3*r4;
        size_t r34 = r3*r4;

        P0 = (double*)malloc(r34*sizeof(double));
        P1 = (double*)malloc(r34*sizeof(double));

        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);
        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);

        int* type = (int*) malloc(dataLength*sizeof(int));

        double* spaceFillingValue = oriData; //

        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);

        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);

        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);

        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));


        size_t l;
        for (l = 0; l < r1; l++)
        {

                ///////////////////////////     Process layer-0 ///////////////////////////
                /* Process Row-0 data 0*/
                size_t index = l*r234;
                size_t index2D = 0;

                type[index] = 0;
                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                P1[index2D] = vce->data;

                /* Process Row-0 data 1*/
                index = l*r234+1;
                index2D = 1;

                pred1D = P1[index2D-1];
                diff = spaceFillingValue[index] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[index] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[index2D] = vce->data;
                }

                /* Process Row-0 data 2 --> data r4-1 */
                for (j = 2; j < r4; j++)
                {
                        index = l*r234+j;
                        index2D = j;

                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
                        diff = spaceFillingValue[index] - pred1D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P1[index2D] = vce->data;
                        }
                }

                /* Process Row-1 --> Row-r3-1 */
                for (i = 1; i < r3; i++)
                {
                        /* Process row-i data 0 */
                        index = l*r234+i*r4;
                        index2D = i*r4;

                        pred1D = P1[index2D-r4];
                        diff = spaceFillingValue[index] - pred1D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P1[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P1[index2D] = vce->data;
                        }

                        /* Process row-i data 1 --> data r4-1*/
                        for (j = 1; j < r4; j++)
                        {
                                index = l*r234+i*r4+j;
                                index2D = i*r4+j;

                                pred2D = P1[index2D-1] + P1[index2D-r4] - P1[index2D-r4-1];

                                diff = spaceFillingValue[index] - pred2D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P1[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[index] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P1[index2D] = vce->data;
                                }
                        }
                }


                ///////////////////////////     Process layer-1 --> layer-r2-1 ///////////////////////////

                for (k = 1; k < r2; k++)
                {
                        /* Process Row-0 data 0*/
                        index = l*r234+k*r34;
                        index2D = 0;

                        pred1D = P1[index2D];
                        diff = spaceFillingValue[index] - pred1D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[index2D] = pred1D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[index] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[index2D] = vce->data;
                        }


                        /* Process Row-0 data 1 --> data r4-1 */
                        for (j = 1; j < r4; j++)
                        {
                                index = l*r234+k*r34+j;
                                index2D = j;

                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
                                diff = spaceFillingValue[index] - pred2D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[index] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P0[index2D] = vce->data;
                                }
                        }

                        /* Process Row-1 --> Row-r3-1 */
                        for (i = 1; i < r3; i++)
                        {
                                /* Process Row-i data 0 */
                                index = l*r234+k*r34+i*r4;
                                index2D = i*r4;

                                pred2D = P0[index2D-r4] + P1[index2D] - P1[index2D-r4];
                                diff = spaceFillingValue[index] - pred2D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P0[index2D] = pred2D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[index] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P0[index2D] = vce->data;
                                }

                                /* Process Row-i data 1 --> data r4-1 */
                                for (j = 1; j < r4; j++)
                                {
                                        index = l*r234+k*r34+i*r4+j;
                                        index2D = i*r4+j;

                                        pred3D = P0[index2D-1] + P0[index2D-r4]+ P1[index2D] - P0[index2D-r4-1] - P1[index2D-r4] - P1[index2D-1] + P1[index2D-r4-1];
                                        diff = spaceFillingValue[index] - pred3D;


                                        itvNum = fabs(diff)/realPrecision + 1;

                                        if (itvNum < exe_params->intvCapacity)
                                        {
                                                if (diff < 0) itvNum = -itvNum;
                                                type[index] = (int) (itvNum/2) + exe_params->intvRadius;
                                                P0[index2D] = pred3D + 2 * (type[index] - exe_params->intvRadius) * realPrecision;
                                        }
                                        else
                                        {
                                                type[index] = 0;
                                                compressSingleDoubleValue(vce, spaceFillingValue[index], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                                memcpy(preDataBytes,vce->curBytes,8);
                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                                P0[index2D] = vce->data;
                                        }
                                }
                        }

                        double *Pt;
                        Pt = P1;
                        P1 = P0;
                        P0 = Pt;
                }
        }

        free(P0);
        free(P1);
        size_t exactDataNum = exactLeadNumArray->size;

        TightDataPointStorageD* tdps;

        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
                        type, exactMidByteArray->array, exactMidByteArray->size,
                        exactLeadNumArray->array,
                        resiBitArray->array, resiBitArray->size,
                        resiBitsLength,
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);

        return tdps;
}


char SZ_compress_args_double_NoCkRngeNoGzip_4D(unsigned char** newByteData, double *oriData, size_t r1, size_t r2, size_t r3, size_t r4, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d)
{
        TightDataPointStorageD* tdps = SZ_compress_double_4D_MDQ(oriData, r1, r2, r3, r4, realPrecision, valueRangeSize, medianValue_d);

        convertTDPStoFlatBytes_double(tdps, newByteData, outSize);

        size_t dataLength = r1*r2*r3*r4;
        if(*outSize>dataLength*sizeof(double))
                SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);

        free_TightDataPointStorageD(tdps);
        return 0;
}

void SZ_compress_args_double_withinRange(unsigned char** newByteData, double *oriData, size_t dataLength, size_t *outSize)
{
        TightDataPointStorageD* tdps = (TightDataPointStorageD*) malloc(sizeof(TightDataPointStorageD));
        tdps->rtypeArray = NULL;
        tdps->typeArray = NULL;
        tdps->leadNumArray = NULL;
        tdps->residualMidBits = NULL;
        
        tdps->allSameData = 1;
        tdps->dataSeriesLength = dataLength;
        tdps->exactMidBytes = (unsigned char*)malloc(sizeof(unsigned char)*8);
        tdps->pwrErrBoundBytes = NULL;
        tdps->isLossless = 0;
        double value = oriData[0];
        doubleToBytes(tdps->exactMidBytes, value);
        tdps->exactMidBytes_size = 8;
        
        size_t tmpOutSize;
        //unsigned char *tmpByteData;
        convertTDPStoFlatBytes_double(tdps, newByteData, &tmpOutSize);
        //convertTDPStoFlatBytes_double(tdps, &tmpByteData, &tmpOutSize);

        //*newByteData = (unsigned char*)malloc(sizeof(unsigned char)*16); //for floating-point data (1+3+4+4)
        //memcpy(*newByteData, tmpByteData, 16);
        *outSize = tmpOutSize;//12==3+1+8(double_size)+MetaDataByteLength
        free_TightDataPointStorageD(tdps);      
}

int SZ_compress_args_double_wRngeNoGzip(unsigned char** newByteData, double *oriData, 
size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
int errBoundMode, double absErr_Bound, double relBoundRatio, double pwrErrRatio)
{
        int status = SZ_SCES;
        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
        double valueRangeSize = 0, medianValue = 0;
        
        double min = computeRangeSize_double(oriData, dataLength, &valueRangeSize, &medianValue);
        double max = min+valueRangeSize;
        double realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
                
        if(valueRangeSize <= realPrecision)
        {
                SZ_compress_args_double_withinRange(newByteData, oriData, dataLength, outSize);
        }
        else
        {
                if(r5==0&&r4==0&&r3==0&&r2==0)
                {
                        if(errBoundMode>=PW_REL)
                        {
                                //SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(newByteData, oriData, realPrecision, r1, outSize, min, max);
                                SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(newByteData, oriData, r1, absErr_Bound, relBoundRatio, pwrErrRatio, valueRangeSize, medianValue, outSize);                           
                        }
                        else
                                SZ_compress_args_double_NoCkRngeNoGzip_1D(newByteData, oriData, r1, realPrecision, outSize, valueRangeSize, medianValue);
                }
                else if(r5==0&&r4==0&&r3==0)
                {
                        if(errBoundMode>=PW_REL)
                                SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(newByteData, oriData, realPrecision, r2, r1, outSize, min, max);
                        else
                                SZ_compress_args_double_NoCkRngeNoGzip_2D(newByteData, oriData, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
                }
                else if(r5==0&&r4==0)
                {
                        if(errBoundMode>=PW_REL)
                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r3, r2, r1, outSize, min, max);
                        else
                                SZ_compress_args_double_NoCkRngeNoGzip_3D(newByteData, oriData, r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
                }
                else if(r5==0)
                {
                        if(errBoundMode>=PW_REL)
                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(newByteData, oriData, realPrecision, r4*r3, r2, r1, outSize, min, max);
                        else
                                SZ_compress_args_double_NoCkRngeNoGzip_3D(newByteData, oriData, r4*r3, r2, r1, realPrecision, outSize, valueRangeSize, medianValue);
                }
        }
        return status;
}

int SZ_compress_args_double(unsigned char** newByteData, double *oriData, 
size_t r5, size_t r4, size_t r3, size_t r2, size_t r1, size_t *outSize, 
int errBoundMode, double absErr_Bound, double relBoundRatio, double pwRelBoundRatio)
{
        confparams_cpr->errorBoundMode = errBoundMode;
        if(errBoundMode==PW_REL)
        {
                confparams_cpr->pw_relBoundRatio = pwRelBoundRatio;     
                //confparams_cpr->pwr_type = SZ_PWR_MIN_TYPE;
                if(confparams_cpr->pwr_type==SZ_PWR_AVG_TYPE && r3 != 0 )
                {
                        printf("Error: Current version doesn't support 3D data compression with point-wise relative error bound being based on pwrType=AVG\n");
                        exit(0);
                        return SZ_NSCS;
                }
        }                               
                
        int status = SZ_SCES;
        size_t dataLength = computeDataLength(r5,r4,r3,r2,r1);
        
        if(dataLength <= MIN_NUM_OF_ELEMENTS)
        {
                *newByteData = SZ_skip_compress_double(oriData, dataLength, outSize);
                return status;
        }
        
        double valueRangeSize = 0, medianValue = 0;
        
        double min = computeRangeSize_double(oriData, dataLength, &valueRangeSize, &medianValue);
        double max = min+valueRangeSize;

        double realPrecision = 0; 
        
        if(confparams_cpr->errorBoundMode==PSNR)
        {
                confparams_cpr->errorBoundMode = ABS;
                realPrecision = confparams_cpr->absErrBound = computeABSErrBoundFromPSNR(confparams_cpr->psnr, (double)confparams_cpr->predThreshold, valueRangeSize);
        }
        else
                realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);
                
        if(valueRangeSize <= realPrecision)
        {
                SZ_compress_args_double_withinRange(newByteData, oriData, dataLength, outSize);
        }
        else
        {
                size_t tmpOutSize = 0;
                unsigned char* tmpByteData;
                if (r2==0)
                {
                        if(confparams_cpr->errorBoundMode>=PW_REL)
                        {
                                //SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr(&tmpByteData, oriData, realPrecision, r1, &tmpOutSize, min, max);
                                SZ_compress_args_double_NoCkRngeNoGzip_1D_pwrgroup(&tmpByteData, oriData, r1, absErr_Bound, relBoundRatio, pwRelBoundRatio, 
                                valueRangeSize, medianValue, &tmpOutSize);
                        }
                        else
#ifdef HAVE_TIMECMPR
                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                   
                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                                else
#endif
                                        SZ_compress_args_double_NoCkRngeNoGzip_1D(&tmpByteData, oriData, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                }
                else
                if (r3==0)
                {
                        if(confparams_cpr->errorBoundMode>=PW_REL)
                                SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr(&tmpByteData, oriData, realPrecision, r2, r1, &tmpOutSize, min, max);
                        else
#ifdef HAVE_TIMECMPR
                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                   
                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                                else
#endif
                                        SZ_compress_args_double_NoCkRngeNoGzip_2D(&tmpByteData, oriData, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                }
                else
                if (r4==0)
                {
                        if(confparams_cpr->errorBoundMode>=PW_REL)
                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r3, r2, r1, &tmpOutSize, min, max);
                        else
#ifdef HAVE_TIMECMPR
                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)
                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                                else
#endif
                                        SZ_compress_args_double_NoCkRngeNoGzip_3D(&tmpByteData, oriData, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                }
                else
                if (r5==0)
                {
                        if(confparams_cpr->errorBoundMode>=PW_REL)
                                SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr(&tmpByteData, oriData, realPrecision, r4*r3, r2, r1, &tmpOutSize, min, max);
                        else
#ifdef HAVE_TIMECMPR
                                if(confparams_cpr->szMode == SZ_TEMPORAL_COMPRESSION)                   
                                        multisteps->compressionType = SZ_compress_args_double_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                                else
#endif
                                        SZ_compress_args_double_NoCkRngeNoGzip_4D(&tmpByteData, oriData, r4, r3, r2, r1, realPrecision, &tmpOutSize, valueRangeSize, medianValue);
                }
                else
                {
                        printf("Error: doesn't support 5 dimensions for now.\n");
                        status = SZ_DERR;
                }
                                
                //Call Gzip to do the further compression.
                if(confparams_cpr->szMode==SZ_BEST_SPEED)
                {
                        *outSize = tmpOutSize;
                        *newByteData = tmpByteData;                     
                }
                else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
                {
                        *outSize = zlib_compress5(tmpByteData, tmpOutSize, newByteData, confparams_cpr->gzipMode);
                        free(tmpByteData);
                }
                else
                {
                        printf("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
                        status = SZ_MERR;       
                }
        }

        return status;
}

//TODO
int SZ_compress_args_double_subblock(unsigned char* compressedBytes, double *oriData,
size_t r5, size_t r4, size_t r3, size_t r2, size_t r1,
size_t s5, size_t s4, size_t s3, size_t s2, size_t s1,
size_t e5, size_t e4, size_t e3, size_t e2, size_t e1,
size_t *outSize, int errBoundMode, double absErr_Bound, double relBoundRatio)
{
        int status = SZ_SCES;
        double valueRangeSize = 0, medianValue = 0;
        computeRangeSize_double_subblock(oriData, &valueRangeSize, &medianValue, r5, r4, r3, r2, r1, s5, s4, s3, s2, s1, e5, e4, e3, e2, e1);

        double realPrecision = getRealPrecision_double(valueRangeSize, errBoundMode, absErr_Bound, relBoundRatio, &status);

        if(valueRangeSize <= realPrecision)
        {
                //TODO
                //SZ_compress_args_double_withinRange_subblock();
        }
        else
        {
                if (r2==0)
                {
                        //TODO
                        if(errBoundMode==PW_REL)
                        {
                                //TODO
                                //SZ_compress_args_double_NoCkRngeNoGzip_1D_pwr_subblock();
                                printf ("Current subblock version does not support point-wise relative error bound.\n");
                        }
                        else
                                SZ_compress_args_double_NoCkRnge_1D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r1, s1, e1);
                }
                else
                if (r3==0)
                {
                        if(errBoundMode==PW_REL)
                        {
                                //TODO
                                //SZ_compress_args_double_NoCkRngeNoGzip_2D_pwr_subblock();
                                printf ("Current subblock version does not support point-wise relative error bound.\n");
                        }
                        else
                                SZ_compress_args_double_NoCkRnge_2D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r2, r1, s2, s1, e2, e1);
                }
                else
                if (r4==0)
                {
                        if(errBoundMode==PW_REL)
                        {
                                //TODO
                                //SZ_compress_args_double_NoCkRngeNoGzip_3D_pwr_subblock();
                                printf ("Current subblock version does not support point-wise relative error bound.\n");
                        }
                        else
                                SZ_compress_args_double_NoCkRnge_3D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r3, r2, r1, s3, s2, s1, e3, e2, e1);
                }
                else
                if (r5==0)
                {
                        if(errBoundMode==PW_REL)
                        {
                                //TODO
                                //SZ_compress_args_double_NoCkRngeNoGzip_4D_pwr_subblock();
                                printf ("Current subblock version does not support point-wise relative error bound.\n");
                        }
                        else
                                SZ_compress_args_double_NoCkRnge_4D_subblock(compressedBytes, oriData, realPrecision, outSize, valueRangeSize, medianValue, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);
                }
                else
                {
                        printf("Error: doesn't support 5 dimensions for now.\n");
                        status = SZ_DERR; //dimension error
                }
        }
        return status;
}

void SZ_compress_args_double_NoCkRnge_1D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
size_t r1, size_t s1, size_t e1)
{
        TightDataPointStorageD* tdps = SZ_compress_double_1D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r1, s1, e1);

        if (confparams_cpr->szMode==SZ_BEST_SPEED)
                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
        {
                unsigned char *tmpCompBytes;
                size_t tmpOutSize;
                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
                free(tmpCompBytes);
        }
        else
        {
                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
        }

        //TODO
//      if(*outSize>dataLength*sizeof(double))
//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);

        free_TightDataPointStorageD(tdps);
}

void SZ_compress_args_double_NoCkRnge_2D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
size_t r2, size_t r1, size_t s2, size_t s1, size_t e2, size_t e1)
{
        TightDataPointStorageD* tdps = SZ_compress_double_2D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r2, r1, s2, s1, e2, e1);

        if (confparams_cpr->szMode==SZ_BEST_SPEED)
                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
        {
                unsigned char *tmpCompBytes;
                size_t tmpOutSize;
                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
                free(tmpCompBytes);
        }
        else
        {
                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
        }

        //TODO
//      if(*outSize>dataLength*sizeof(double))
//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);

        free_TightDataPointStorageD(tdps);
}

void SZ_compress_args_double_NoCkRnge_3D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
size_t r3, size_t r2, size_t r1, size_t s3, size_t s2, size_t s1, size_t e3, size_t e2, size_t e1)
{
        TightDataPointStorageD* tdps = SZ_compress_double_3D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r3, r2, r1, s3, s2, s1, e3, e2, e1);

        if (confparams_cpr->szMode==SZ_BEST_SPEED)
                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
        {
                unsigned char *tmpCompBytes;
                size_t tmpOutSize;
                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
                free(tmpCompBytes);
        }
        else
        {
                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
        }

        //TODO
//      if(*outSize>dataLength*sizeof(double))
//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);

        free_TightDataPointStorageD(tdps);
}

void SZ_compress_args_double_NoCkRnge_4D_subblock(unsigned char* compressedBytes, double *oriData, double realPrecision, size_t *outSize, double valueRangeSize, double medianValue_d,
size_t r4, size_t r3, size_t r2, size_t r1, size_t s4, size_t s3, size_t s2, size_t s1, size_t e4, size_t e3, size_t e2, size_t e1)
{
        TightDataPointStorageD* tdps = SZ_compress_double_4D_MDQ_subblock(oriData, realPrecision, valueRangeSize, medianValue_d, r4, r3, r2, r1, s4, s3, s2, s1, e4, e3, e2, e1);

        if (confparams_cpr->szMode==SZ_BEST_SPEED)
                convertTDPStoFlatBytes_double_args(tdps, compressedBytes, outSize);
        else if(confparams_cpr->szMode==SZ_BEST_COMPRESSION || confparams_cpr->szMode==SZ_DEFAULT_COMPRESSION)
        {
                unsigned char *tmpCompBytes;
                size_t tmpOutSize;
                convertTDPStoFlatBytes_double(tdps, &tmpCompBytes, &tmpOutSize);
                *outSize = zlib_compress3(tmpCompBytes, tmpOutSize, compressedBytes, confparams_cpr->gzipMode);
                free(tmpCompBytes);
        }
        else
        {
                printf ("Error: Wrong setting of confparams_cpr->szMode in the double compression.\n");
        }

        //TODO
//      if(*outSize>dataLength*sizeof(double))
//              SZ_compress_args_double_StoreOriData(oriData, dataLength, tdps, newByteData, outSize);

        free_TightDataPointStorageD(tdps);
}


unsigned int optimize_intervals_double_1D_subblock(double *oriData, double realPrecision, size_t r1, size_t s1, size_t e1)
{
        size_t dataLength = e1 - s1 + 1;
        oriData = oriData + s1;

        size_t i = 0;
        unsigned long radiusIndex;
        double pred_value = 0, pred_err;
        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
        size_t totalSampleSize = dataLength/confparams_cpr->sampleDistance;
        for(i=2;i<dataLength;i++)
        {
                if(i%confparams_cpr->sampleDistance==0)
                {
                        pred_value = 2*oriData[i-1] - oriData[i-2];
                        //pred_value = oriData[i-1];
                        pred_err = fabs(pred_value - oriData[i]);
                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
                        intervals[radiusIndex]++;
                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }

        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

unsigned int optimize_intervals_double_2D_subblock(double *oriData, double realPrecision, size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
{
        size_t R1 = e1 - s1 + 1;
        size_t R2 = e2 - s2 + 1;

        size_t i,j, index;
        unsigned long radiusIndex;
        double pred_value = 0, pred_err;
        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
        size_t totalSampleSize = R1*R2/confparams_cpr->sampleDistance;
        for(i=s1+1;i<=e1;i++)
        {
                for(j=s2+1;j<=e2;j++)
                {
                        if((i+j)%confparams_cpr->sampleDistance==0)
                        {
                                index = i*r2+j;
                                pred_value = oriData[index-1] + oriData[index-r2] - oriData[index-r2-1];
                                pred_err = fabs(pred_value - oriData[index]);
                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
                                intervals[radiusIndex]++;
                        }
                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

unsigned int optimize_intervals_double_3D_subblock(double *oriData, double realPrecision, size_t r1, size_t r2, size_t r3, size_t s1, size_t s2, size_t s3, size_t e1, size_t e2, size_t e3)
{
        size_t R1 = e1 - s1 + 1;
        size_t R2 = e2 - s2 + 1;
        size_t R3 = e3 - s3 + 1;

        size_t r23 = r2*r3;

        size_t i,j,k, index;
        unsigned long radiusIndex;
        double pred_value = 0, pred_err;
        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
        size_t totalSampleSize = R1*R2*R3/confparams_cpr->sampleDistance;
        for(i=s1+1;i<=e1;i++)
        {
                for(j=s2+1;j<=e2;j++)
                {
                        for(k=s3+1;k<=e3;k++)
                        {
                                if((i+j+k)%confparams_cpr->sampleDistance==0)
                                {
                                        index = i*r23+j*r3+k;
                                        pred_value = oriData[index-1] + oriData[index-r3] + oriData[index-r23]
                                        - oriData[index-1-r23] - oriData[index-r3-1] - oriData[index-r3-r23] + oriData[index-r3-r23-1];
                                        pred_err = fabs(pred_value - oriData[index]);
                                        radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                                        if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                                radiusIndex = confparams_cpr->maxRangeRadius - 1;
                                        intervals[radiusIndex]++;
                                }
                        }

                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;

        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

unsigned int optimize_intervals_double_4D_subblock(double *oriData, double realPrecision,
size_t r1, size_t r2, size_t r3, size_t r4, size_t s1, size_t s2, size_t s3, size_t s4, size_t e1, size_t e2, size_t e3, size_t e4)
{
        size_t R1 = e1 - s1 + 1;
        size_t R2 = e2 - s2 + 1;
        size_t R3 = e3 - s3 + 1;
        size_t R4 = e4 - s4 + 1;

        size_t r34 = r3*r4;
        size_t r234 = r2*r3*r4;

        size_t i,j,k,l, index;
        unsigned long radiusIndex;
        double pred_value = 0, pred_err;
        int *intervals = (int*)malloc(confparams_cpr->maxRangeRadius*sizeof(int));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(int));
        size_t totalSampleSize = R1*R2*R3*R4/confparams_cpr->sampleDistance;
        for(i=s1+1;i<=e1;i++)
        {
                for(j=s2+1;j<=e2;j++)
                {
                        for(k=s3+1;k<=e3;k++)
                        {
                                for(l=s4+1;l<=e4;l++)
                                {
                                        if((i+j+k+l)%confparams_cpr->sampleDistance==0)
                                        {
                                                index = i*r234+j*r34+k*r4+l;
                                                pred_value = oriData[index-1] + oriData[index-r4] + oriData[index-r34]
                                                                - oriData[index-1-r34] - oriData[index-r4-1] - oriData[index-r4-r34] + oriData[index-r4-r34-1];
                                                pred_err = fabs(pred_value - oriData[index]);
                                                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                                                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                                                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
                                                intervals[radiusIndex]++;
                                        }
                                }
                        }

                }
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;

        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

TightDataPointStorageD* SZ_compress_double_1D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
size_t r1, size_t s1, size_t e1)
{
        size_t dataLength = e1 - s1 + 1;

        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
                quantization_intervals = optimize_intervals_double_1D_subblock(oriData, realPrecision, r1, s1, e1);
        else
                quantization_intervals = exe_params->intvCapacity;
        updateQuantizationInfo(quantization_intervals);

        size_t i; 
        int reqLength;
        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);

        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);

        int* type = (int*) malloc(dataLength*sizeof(int));

        double* spaceFillingValue = oriData + s1; //

        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);

        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);

        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        type[0] = 0;

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);

        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;
        double last3CmprsData[3] = {0};

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));

        //add the first data
        compressSingleDoubleValue(vce, spaceFillingValue[0], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        listAdd_double(last3CmprsData, vce->data);

        //add the second data
        type[1] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[1], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        listAdd_double(last3CmprsData, vce->data);

        int state;
        double checkRadius;
        double curData;
        double pred;
        double predAbsErr;
        checkRadius = (exe_params->intvCapacity-1)*realPrecision;
        double interval = 2*realPrecision;

        for(i=2;i<dataLength;i++)
        {
                //printf("%.30G\n",last3CmprsData[0]);
                curData = spaceFillingValue[i];
                pred = 2*last3CmprsData[0] - last3CmprsData[1];
                //pred = last3CmprsData[0];
                predAbsErr = fabs(curData - pred);
                if(predAbsErr<=checkRadius)
                {
                        state = (predAbsErr/realPrecision+1)/2;
                        if(curData>=pred)
                        {
                                type[i] = exe_params->intvRadius+state;
                                pred = pred + state*interval;
                        }
                        else //curData<pred
                        {
                                type[i] = exe_params->intvRadius-state;
                                pred = pred - state*interval;
                        }
                        listAdd_double(last3CmprsData, pred);
                        continue;
                }

                //unpredictable data processing
                type[i] = 0;
                compressSingleDoubleValue(vce, curData, realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);

                listAdd_double(last3CmprsData, vce->data);
        }//end of for

        size_t exactDataNum = exactLeadNumArray->size;

        TightDataPointStorageD* tdps;

        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
                        type, exactMidByteArray->array, exactMidByteArray->size,
                        exactLeadNumArray->array,
                        resiBitArray->array, resiBitArray->size,
                        resiBitsLength,
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);

        return tdps;
}


TightDataPointStorageD* SZ_compress_double_2D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
size_t r1, size_t r2, size_t s1, size_t s2, size_t e1, size_t e2)
{
        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
        {
                quantization_intervals = optimize_intervals_double_2D_subblock(oriData, realPrecision, r1, r2, s1, s2, e1, e2);
                updateQuantizationInfo(quantization_intervals);
        }
        else
                quantization_intervals = exe_params->intvCapacity;

        size_t i,j; 
        int reqLength;
        double pred1D, pred2D;
        double diff = 0.0;
        double itvNum = 0;
        double *P0, *P1;

        size_t R1 = e1 - s1 + 1;
        size_t R2 = e2 - s2 + 1;
        size_t dataLength = R1*R2;

        P0 = (double*)malloc(R2*sizeof(double));
        memset(P0, 0, R2*sizeof(double));
        P1 = (double*)malloc(R2*sizeof(double));
        memset(P1, 0, R2*sizeof(double));

        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);
        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);

        int* type = (int*) malloc(dataLength*sizeof(int));

        double* spaceFillingValue = oriData; //

        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);

        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);

        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);

        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));

        /* Process Row-s1 data s2*/
        size_t gIndex;
        size_t lIndex;

        gIndex = s1*r2+s2;
        lIndex = 0;

        type[lIndex] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        P1[0] = vce->data;

        /* Process Row-s1 data s2+1*/
        gIndex = s1*r2+(s2+1);
        lIndex = 1;

        pred1D = P1[0];
        diff = spaceFillingValue[gIndex] - pred1D;

        itvNum =  fabs(diff)/realPrecision + 1;

        if (itvNum < exe_params->intvCapacity)
        {
                if (diff < 0) itvNum = -itvNum;
                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                P1[1] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
        }
        else
        {
                type[lIndex] = 0;
                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                P1[1] = vce->data;
        }

    /* Process Row-s1 data s2+2 --> data e2 */
        for (j = 2; j < R2; j++)
        {
                gIndex = s1*r2+(s2+j);
                lIndex = j;

                pred1D = 2*P1[j-1] - P1[j-2];
                diff = spaceFillingValue[gIndex] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[j] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[lIndex] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[j] = vce->data;
                }
        }

        /* Process Row-s1+1 --> Row-e1 */
        for (i = 1; i < R1; i++)
        {
                /* Process row-s1+i data s2 */
                gIndex = (s1+i)*r2+s2;
                lIndex = i*R2;

                pred1D = P1[0];
                diff = spaceFillingValue[gIndex] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                        P0[0] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[lIndex] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P0[0] = vce->data;
                }

                /* Process row-s1+i data s2+1 --> e2 */
                for (j = 1; j < R2; j++)
                {
                        gIndex = (s1+i)*r2+(s2+j);
                        lIndex = i*R2+j;

                        pred2D = P0[j-1] + P1[j] - P1[j-1];
                        diff = spaceFillingValue[gIndex] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[j] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[j] = vce->data;
                        }
                }

                double *Pt;
                Pt = P1;
                P1 = P0;
                P0 = Pt;
        }

        free(P0);
        free(P1);
        size_t exactDataNum = exactLeadNumArray->size;

        TightDataPointStorageD* tdps;

        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
                        type, exactMidByteArray->array, exactMidByteArray->size,
                        exactLeadNumArray->array,
                        resiBitArray->array, resiBitArray->size,
                        resiBitsLength,
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);

        return tdps;
}

TightDataPointStorageD* SZ_compress_double_3D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
size_t r1, size_t r2, size_t r3, size_t s1, size_t s2, size_t s3, size_t e1, size_t e2, size_t e3)
{
        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
        {
                quantization_intervals = optimize_intervals_double_3D_subblock(oriData, realPrecision, r1, r2, r3, s1, s2, s3, e1, e2, e3);
                updateQuantizationInfo(quantization_intervals);
        }
        else
                quantization_intervals = exe_params->intvCapacity;

        size_t i,j,k; 
        int reqLength;
        double pred1D, pred2D, pred3D;
        double diff = 0.0;
        double itvNum = 0;
        double *P0, *P1;

        size_t R1 = e1 - s1 + 1;
        size_t R2 = e2 - s2 + 1;
        size_t R3 = e3 - s3 + 1;
        size_t dataLength = R1*R2*R3;

        size_t r23 = r2*r3;
        size_t R23 = R2*R3;

        P0 = (double*)malloc(R23*sizeof(double));
        P1 = (double*)malloc(R23*sizeof(double));

        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);
        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);

        int* type = (int*) malloc(dataLength*sizeof(int));

        double* spaceFillingValue = oriData; //

        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);

        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);

        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);

        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));


        ///////////////////////////     Process layer-s1 ///////////////////////////
        /* Process Row-s2 data s3*/
        size_t gIndex;  //global index
        size_t lIndex;  //local index
        size_t index2D;         //local 2D index

        gIndex = s1*r23+s2*r3+s3;
        lIndex = 0;
        index2D = 0;

        type[lIndex] = 0;
        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
        memcpy(preDataBytes,vce->curBytes,8);
        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
        P1[index2D] = vce->data;

        /* Process Row-s2 data s3+1*/
        gIndex = s1*r23+s2*r3+s3+1;
        lIndex = 1;
        index2D = 1;

        pred1D = P1[index2D-1];
        diff = spaceFillingValue[gIndex] - pred1D;

        itvNum = fabs(diff)/realPrecision + 1;

        if (itvNum < exe_params->intvCapacity)
        {
                if (diff < 0) itvNum = -itvNum;
                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
        }
        else
        {
                type[lIndex] = 0;
                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                P1[index2D] = vce->data;
        }

    /* Process Row-s2 data s3+2 --> data e3 */
        for (j = 2; j < R3; j++)
        {
                gIndex = s1*r23+s2*r3+s3+j;
                lIndex = j;
                index2D = j;

                pred1D = 2*P1[index2D-1] - P1[index2D-2];
                diff = spaceFillingValue[gIndex] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[lIndex] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[index2D] = vce->data;
                }
        }

        /* Process Row-s2+1 --> Row-e2 */
        for (i = 1; i < R2; i++)
        {
                /* Process row-s2+i data s3 */
                gIndex = s1*r23+(s2+i)*r3+s3;
                lIndex = i*R3;
                index2D = i*R3;

                pred1D  = P1[index2D-R3];
                diff    = spaceFillingValue[gIndex] - pred1D;

                itvNum  = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[lIndex] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[index2D] = vce->data;
                }

                /* Process row-s2+i data s3+1 --> data e3*/
                for (j = 1; j < R3; j++)
                {
                        gIndex = s1*r23+(s2+i)*r3+s3+j;
                        lIndex = i*R3+j;
                        index2D = i*R3+j;

                        pred2D  = P1[index2D-1] + P1[index2D-R3] - P1[index2D-R3-1];
                        diff = spaceFillingValue[gIndex] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P1[index2D] = vce->data;
                        }
                }
        }


        ///////////////////////////     Process layer-s1+1 --> layer-e1 ///////////////////////////

        for (k = 1; k < R1; k++)
        {
                /* Process Row-s2 data s3*/
                gIndex = (s1+k)*r23+s2*r3+s3;
                lIndex = k*R23;
                index2D = 0;

                pred1D = P1[index2D];
                diff = spaceFillingValue[gIndex] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                        P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[lIndex] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P0[index2D] = vce->data;
                }


            /* Process Row-s2 data s3+1 --> data e3 */
                for (j = 1; j < R3; j++)
                {
                        gIndex = (s1+k)*r23+s2*r3+s3+j;
                        lIndex = k*R23+j;
                        index2D = j;

                        pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
                        diff = spaceFillingValue[gIndex] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[index2D] = vce->data;
                        }
                }

            /* Process Row-s2+1 --> Row-e2 */
                for (i = 1; i < R2; i++)
                {
                        /* Process Row-s2+i data s3 */
                        gIndex = (s1+k)*r23+(s2+i)*r3+s3;
                        lIndex = k*R23+i*R3;
                        index2D = i*R3;

                        pred2D = P0[index2D-R3] + P1[index2D] - P1[index2D-R3];
                        diff = spaceFillingValue[gIndex] - pred2D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[index2D] = vce->data;
                        }

                        /* Process Row-s2+i data s3+1 --> data e3 */
                        for (j = 1; j < R3; j++)
                        {
                                gIndex = (s1+k)*r23+(s2+i)*r3+s3+j;
                                lIndex = k*R23+i*R3+j;
                                index2D = i*R3+j;

                                pred3D = P0[index2D-1] + P0[index2D-R3]+ P1[index2D] - P0[index2D-R3-1] - P1[index2D-R3] - P1[index2D-1] + P1[index2D-R3-1];
                                diff = spaceFillingValue[gIndex] - pred3D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[lIndex] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P0[index2D] = vce->data;
                                }
                        }
                }

                double *Pt;
                Pt = P1;
                P1 = P0;
                P0 = Pt;
        }

        free(P0);
        free(P1);
        size_t exactDataNum = exactLeadNumArray->size;

        TightDataPointStorageD* tdps;

        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
                        type, exactMidByteArray->array, exactMidByteArray->size,
                        exactLeadNumArray->array,
                        resiBitArray->array, resiBitArray->size,
                        resiBitsLength,
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);

        return tdps;
}

TightDataPointStorageD* SZ_compress_double_4D_MDQ_subblock(double *oriData, double realPrecision, double valueRangeSize, double medianValue_d,
size_t r1, size_t r2, size_t r3, size_t r4, size_t s1, size_t s2, size_t s3, size_t s4, size_t e1, size_t e2, size_t e3, size_t e4)
{
        unsigned int quantization_intervals;
        if(exe_params->optQuantMode==1)
        {
                quantization_intervals = optimize_intervals_double_4D_subblock(oriData, realPrecision, r1, r2, r3, r4, s1, s2, s3, s4, e1, e2, e3, e4);
                updateQuantizationInfo(quantization_intervals);
        }
        else
                quantization_intervals = exe_params->intvCapacity;

        size_t i,j,k; 
        int reqLength;
        double pred1D, pred2D, pred3D;
        double diff = 0.0;
        double itvNum = 0;
        double *P0, *P1;

        size_t R1 = e1 - s1 + 1;
        size_t R2 = e2 - s2 + 1;
        size_t R3 = e3 - s3 + 1;
        size_t R4 = e4 - s4 + 1;

        size_t dataLength = R1*R2*R3*R4;

        size_t r34 = r3*r4;
        size_t r234 = r2*r3*r4;
        size_t R34 = R3*R4;
        size_t R234 = R2*R3*R4;

        P0 = (double*)malloc(R34*sizeof(double));
        P1 = (double*)malloc(R34*sizeof(double));

        double medianValue = medianValue_d;
        short radExpo = getExponent_double(valueRangeSize/2);
        computeReqLength_double(realPrecision, radExpo, &reqLength, &medianValue);

        int* type = (int*) malloc(dataLength*sizeof(int));

        double* spaceFillingValue = oriData; //

        DynamicIntArray *exactLeadNumArray;
        new_DIA(&exactLeadNumArray, DynArrayInitLen);

        DynamicByteArray *exactMidByteArray;
        new_DBA(&exactMidByteArray, DynArrayInitLen);

        DynamicIntArray *resiBitArray;
        new_DIA(&resiBitArray, DynArrayInitLen);

        unsigned char preDataBytes[8];
        longToBytes_bigEndian(preDataBytes, 0);

        int reqBytesLength = reqLength/8;
        int resiBitsLength = reqLength%8;

        DoubleValueCompressElement *vce = (DoubleValueCompressElement*)malloc(sizeof(DoubleValueCompressElement));
        LossyCompressionElement *lce = (LossyCompressionElement*)malloc(sizeof(LossyCompressionElement));

        size_t l;
        for (l = 0; l < R1; l++)
        {

                ///////////////////////////     Process layer-s2 ///////////////////////////
                /* Process Row-s3 data s4*/
                size_t gIndex;  //global index
                size_t lIndex;  //local index
                size_t index2D;         //local 2D index

                gIndex = (s1+l)*r234+s2*r34+s3*r4+s4;
                lIndex = l*R234;
                index2D = 0;

                type[lIndex] = 0;
                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                memcpy(preDataBytes,vce->curBytes,8);
                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                P1[index2D] = vce->data;

                /* Process Row-s3 data s4+1*/
                gIndex = (s1+l)*r234+s2*r34+s3*r4+s4+1;
                lIndex = l*R234+1;
                index2D = 1;

                pred1D = P1[index2D-1];
                diff = spaceFillingValue[gIndex] - pred1D;

                itvNum = fabs(diff)/realPrecision + 1;

                if (itvNum < exe_params->intvCapacity)
                {
                        if (diff < 0) itvNum = -itvNum;
                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                        P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                }
                else
                {
                        type[lIndex] = 0;
                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                        memcpy(preDataBytes,vce->curBytes,8);
                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                        P1[index2D] = vce->data;
                }

                /* Process Row-s3 data s4+2 --> data e4 */
                for (j = 2; j < R4; j++)
                {
                        gIndex = (s1+l)*r234+s2*r34+s3*r4+s4+j;
                        lIndex = l*R234+j;
                        index2D = j;

                        pred1D = 2*P1[index2D-1] - P1[index2D-2];
                        diff = spaceFillingValue[gIndex] - pred1D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P1[index2D] = vce->data;
                        }
                }

                /* Process Row-s3+1 --> Row-e3 */
                for (i = 1; i < R3; i++)
                {
                        /* Process row-s2+i data s3 */
                        gIndex = (s1+l)*r234+s2*r34+(s3+i)*r4+s4;
                        lIndex = l*R234+i*R4;
                        index2D = i*R4;

                        pred1D  = P1[index2D-R4];
                        diff    = spaceFillingValue[gIndex] - pred1D;

                        itvNum  = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P1[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P1[index2D] = vce->data;
                        }

                        /* Process row-s3+i data s4+1 --> data e4*/
                        for (j = 1; j < R4; j++)
                        {
                                gIndex = (s1+l)*r234+s2*r34+(s3+i)*r4+s4+j;
                                lIndex = l*R234+i*R4+j;
                                index2D = i*R4+j;

                                pred2D  = P1[index2D-1] + P1[index2D-R4] - P1[index2D-R4-1];
                                diff = spaceFillingValue[gIndex] - pred2D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P1[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[lIndex] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P1[index2D] = vce->data;
                                }
                        }
                }


                ///////////////////////////     Process layer-s2+1 --> layer-e2 ///////////////////////////

                for (k = 1; k < R2; k++)
                {
                        /* Process Row-s3 data s4*/
                        gIndex = (s1+l)*r234+(s2+k)*r34+s3*r4+s4;
                        lIndex = l*R234+k*R34;
                        index2D = 0;

                        pred1D = P1[index2D];
                        diff = spaceFillingValue[gIndex] - pred1D;

                        itvNum = fabs(diff)/realPrecision + 1;

                        if (itvNum < exe_params->intvCapacity)
                        {
                                if (diff < 0) itvNum = -itvNum;
                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                P0[index2D] = pred1D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                        }
                        else
                        {
                                type[lIndex] = 0;
                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                memcpy(preDataBytes,vce->curBytes,8);
                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                P0[index2D] = vce->data;
                        }


                        /* Process Row-s3 data s4+1 --> data e4 */
                        for (j = 1; j < R4; j++)
                        {
                                gIndex = (s1+l)*r234+(s2+k)*r34+s3*r4+s4+j;
                                lIndex = l*R234+k*R34+j;
                                index2D = j;

                                pred2D = P0[index2D-1] + P1[index2D] - P1[index2D-1];
                                diff = spaceFillingValue[gIndex] - pred2D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[lIndex] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P0[index2D] = vce->data;
                                }
                        }

                        /* Process Row-s3+1 --> Row-e3 */
                        for (i = 1; i < R3; i++)
                        {
                                /* Process Row-s3+i data s4 */
                                gIndex = (s1+l)*r234+(s2+k)*r34+(s3+i)*r4+s4;
                                lIndex = l*R234+k*R34+i*R4;
                                index2D = i*R4;

                                pred2D = P0[index2D-R4] + P1[index2D] - P1[index2D-R4];
                                diff = spaceFillingValue[gIndex] - pred2D;

                                itvNum = fabs(diff)/realPrecision + 1;

                                if (itvNum < exe_params->intvCapacity)
                                {
                                        if (diff < 0) itvNum = -itvNum;
                                        type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                        P0[index2D] = pred2D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                                }
                                else
                                {
                                        type[lIndex] = 0;
                                        compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                        updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                        memcpy(preDataBytes,vce->curBytes,8);
                                        addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                        P0[index2D] = vce->data;
                                }

                                /* Process Row-s3+i data s4+1 --> data e4 */
                                for (j = 1; j < R4; j++)
                                {
                                        gIndex = (s1+l)*r234+(s2+k)*r34+(s3+i)*r4+s4+j;
                                        lIndex = l*R234+k*R34+i*R4+j;
                                        index2D = i*R4+j;

//                                      printf ("global index = %d, local index = %d\n", gIndex, lIndex);

                                        pred3D = P0[index2D-1] + P0[index2D-R4]+ P1[index2D] - P0[index2D-R4-1] - P1[index2D-R4] - P1[index2D-1] + P1[index2D-R4-1];
                                        diff = spaceFillingValue[gIndex] - pred3D;

                                        itvNum = fabs(diff)/realPrecision + 1;

                                        if (itvNum < exe_params->intvCapacity)
                                        {
                                                if (diff < 0) itvNum = -itvNum;
                                                type[lIndex] = (int) (itvNum/2) + exe_params->intvRadius;
                                                P0[index2D] = pred3D + 2 * (type[lIndex] - exe_params->intvRadius) * realPrecision;
                                        }
                                        else
                                        {
                                                type[lIndex] = 0;
                                                compressSingleDoubleValue(vce, spaceFillingValue[gIndex], realPrecision, medianValue, reqLength, reqBytesLength, resiBitsLength);
                                                updateLossyCompElement_Double(vce->curBytes, preDataBytes, reqBytesLength, resiBitsLength, lce);
                                                memcpy(preDataBytes,vce->curBytes,8);
                                                addExactData(exactMidByteArray, exactLeadNumArray, resiBitArray, lce);
                                                P0[index2D] = vce->data;
                                        }
                                }
                        }

                        double *Pt;
                        Pt = P1;
                        P1 = P0;
                        P0 = Pt;
                }
        }

        free(P0);
        free(P1);
        size_t exactDataNum = exactLeadNumArray->size;

        TightDataPointStorageD* tdps;

        new_TightDataPointStorageD(&tdps, dataLength, exactDataNum,
                        type, exactMidByteArray->array, exactMidByteArray->size,
                        exactLeadNumArray->array,
                        resiBitArray->array, resiBitArray->size,
                        resiBitsLength,
                        realPrecision, medianValue, (char)reqLength, quantization_intervals, NULL, 0, 0);

        //free memory
        free_DIA(exactLeadNumArray);
        free_DIA(resiBitArray);
        free(type);
        free(vce);
        free(lce);
        free(exactMidByteArray); //exactMidByteArray->array has been released in free_TightDataPointStorageF(tdps);

        return tdps;
}

/**
 * 
 * This is a fast implementation for optimize_intervals_double_3D()
 * */
unsigned int optimize_intervals_double_3D_opt(double *oriData, size_t r1, size_t r2, size_t r3, double realPrecision){  
        size_t i;
        size_t radiusIndex;
        size_t r23=r2*r3;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = 0;

        size_t offset_count = confparams_cpr->sampleDistance - 2; // count r3 offset
        size_t offset_count_2;
        double * data_pos = oriData + r23 + r3 + offset_count;
        size_t n1_count = 1, n2_count = 1; // count i,j sum
        size_t len = r1 * r2 * r3;
        while(data_pos - oriData < len){
                totalSampleSize++;
                pred_value = data_pos[-1] + data_pos[-r3] + data_pos[-r23] - data_pos[-1-r23] - data_pos[-r3-1] - data_pos[-r3-r23] + data_pos[-r3-r23-1];
                pred_err = fabs(pred_value - *data_pos);
                radiusIndex = (pred_err/realPrecision+1)/2;
                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                {
                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
                        //printf("radiusIndex=%d\n", radiusIndex);
                }
                intervals[radiusIndex]++;
                // printf("TEST: %ld, i: %ld\tj: %ld\tk: %ld\n", data_pos - oriData);
                // fflush(stdout);
                offset_count += confparams_cpr->sampleDistance;
                if(offset_count >= r3){
                        n2_count ++;
                        if(n2_count == r2){
                                n1_count ++;
                                n2_count = 1;
                                data_pos += r3;
                        }
                        offset_count_2 = (n1_count + n2_count) % confparams_cpr->sampleDistance;
                        data_pos += (r3 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
                        offset_count = (confparams_cpr->sampleDistance - offset_count_2);
                        if(offset_count == 0) offset_count ++;
                }
                else data_pos += confparams_cpr->sampleDistance;
        }       
        // printf("sample_count: %ld\n", sample_count);
        // fflush(stdout);
        // if(*max_freq < 0.15) *max_freq *= 2;
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;
        free(intervals);
        //printf("targetCount=%d, sum=%d, totalSampleSize=%d, ratio=%f, accIntervals=%d, powerOf2=%d\n", targetCount, sum, totalSampleSize, (double)sum/(double)totalSampleSize, accIntervals, powerOf2);
        return powerOf2;
}

unsigned int optimize_intervals_double_2D_opt(double *oriData, size_t r1, size_t r2, double realPrecision)
{       
        size_t i;
        size_t radiusIndex;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = 0;//(r1-1)*(r2-1)/confparams_cpr->sampleDistance;

        size_t offset_count = confparams_cpr->sampleDistance - 1; // count r2 offset
        size_t offset_count_2;
        double * data_pos = oriData + r2 + offset_count;
        size_t n1_count = 1; // count i sum
        size_t len = r1 * r2;
        while(data_pos - oriData < len){
                totalSampleSize++;
                pred_value = data_pos[-1] + data_pos[-r2] - data_pos[-r2-1];
                pred_err = fabs(pred_value - *data_pos);
                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                        radiusIndex = confparams_cpr->maxRangeRadius - 1;
                intervals[radiusIndex]++;

                offset_count += confparams_cpr->sampleDistance;
                if(offset_count >= r2){
                        n1_count ++;
                        offset_count_2 = n1_count % confparams_cpr->sampleDistance;
                        data_pos += (r2 + confparams_cpr->sampleDistance - offset_count) + (confparams_cpr->sampleDistance - offset_count_2);
                        offset_count = (confparams_cpr->sampleDistance - offset_count_2);
                        if(offset_count == 0) offset_count ++;
                }
                else data_pos += confparams_cpr->sampleDistance;
        }

        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);

        if(powerOf2<32)
                powerOf2 = 32;

        free(intervals);
        return powerOf2;
}

unsigned int optimize_intervals_double_1D_opt(double *oriData, size_t dataLength, double realPrecision)
{       
        size_t i = 0, radiusIndex;
        double pred_value = 0, pred_err;
        size_t *intervals = (size_t*)malloc(confparams_cpr->maxRangeRadius*sizeof(size_t));
        memset(intervals, 0, confparams_cpr->maxRangeRadius*sizeof(size_t));
        size_t totalSampleSize = 0;//dataLength/confparams_cpr->sampleDistance;

        double * data_pos = oriData + 2;
        while(data_pos - oriData < dataLength){
                totalSampleSize++;
                //pred_value = 2*data_pos[-1] - data_pos[-2];
                pred_value = data_pos[-1];
                pred_err = fabs(pred_value - *data_pos);
                radiusIndex = (unsigned long)((pred_err/realPrecision+1)/2);
                if(radiusIndex>=confparams_cpr->maxRangeRadius)
                        radiusIndex = confparams_cpr->maxRangeRadius - 1;                       
                intervals[radiusIndex]++;

                data_pos += confparams_cpr->sampleDistance;
        }
        //compute the appropriate number
        size_t targetCount = totalSampleSize*confparams_cpr->predThreshold;
        size_t sum = 0;
        for(i=0;i<confparams_cpr->maxRangeRadius;i++)
        {
                sum += intervals[i];
                if(sum>targetCount)
                        break;
        }
        if(i>=confparams_cpr->maxRangeRadius)
                i = confparams_cpr->maxRangeRadius-1;
                
        unsigned int accIntervals = 2*(i+1);
        unsigned int powerOf2 = roundUpToPowerOf2(accIntervals);
        
        if(powerOf2<32)
                powerOf2 = 32;
        
        free(intervals);
        //printf("accIntervals=%d, powerOf2=%d\n", accIntervals, powerOf2);
        return powerOf2;
}