Changeset 13577


Ignore:
Timestamp:
10/10/12 07:39:19 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: even simpler binary_search for doubles: return segment number

Location:
issm/trunk-jpl/src/c
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp

    r13574 r13577  
    426426Input* TransientInput::GetTimeInput(IssmDouble intime){
    427427
    428         IssmDouble  deltat;
    429         IssmDouble  alpha1,alpha2;
     428        IssmDouble deltat;
     429        IssmDouble alpha1,alpha2;
     430        int        found;
     431        int        offset;
    430432       
    431433        Input *input  = NULL;
     
    433435        Input *input2 = NULL;
    434436       
    435         bool    found=false;
    436         int     found_offset=0;
    437         int     offset;
    438 
    439         /*Ok, we have the time, go through the timesteps, and figure out which interval we
     437        /*go through the timesteps, and figure out which interval we
    440438         *fall within. Then interpolate the values on this interval: */
    441         found_offset=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
    442         if (found_offset==0) _error_("Input not found (is TransientInput sorted ?)");
    443 
    444         if (found_offset==-1){
     439        found=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
     440        if(!found) _error_("Input not found (is TransientInput sorted ?)");
     441
     442        if (offset==-1){
    445443                /*get values for the first time: */
    446444                _assert_(intime<this->timesteps[0]);
    447445                input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
    448                 found=true;
    449         }
    450         else if(found_offset==3){
     446        }
     447        else if(offset==(this->numtimesteps-1)){
    451448                /*get values for the last time: */
    452                 _assert_(intime>=this->timesteps[this->numtimesteps-1]);
    453                 input=(Input*)((Input*)this->inputs->GetObjectByOffset(this->numtimesteps-1))->copy();
    454                 found=true;
    455         }
    456         else if(found_offset==2){
    457                 /*get values for this time: */
    458                 _assert_(intime==this->timesteps[offset]);
     449                _assert_(intime>=this->timesteps[offset]);
    459450                input=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
    460                 found=true;
    461451        }
    462452        else{
    463                 /*get values between two times ]offset:offset+1[, Interpolate linearly*/
    464                 _assert_(intime>this->timesteps[offset] && intime<this->timesteps[offset+1]);
     453                /*get values between two times [offset:offset+1[, Interpolate linearly*/
     454                _assert_(intime>=this->timesteps[offset] && intime<this->timesteps[offset+1]);
    465455                deltat=this->timesteps[offset+1]-this->timesteps[offset];
    466456                alpha2=(intime-this->timesteps[offset])/deltat;
     
    473463                input->Scale(alpha1);
    474464                input->AXPY(input2,alpha2);
    475 
    476                 found=true;
    477465        }
    478466
  • issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp

    r13575 r13577  
    6464        return found;
    6565} /*}}}*/
    66 int binary_search(int* poffset,double target,double* list,int num_doubles){ /*{{{*/
     66int binary_search(int* poffset,double target,double* list,int length){ /*{{{*/
     67        /*
     68         *             l[0]  l[1]  l[2]        l[n]  l[n+1]   l[length-1]
     69         *     <-------+-----+-----+-----+ ... +-----+........+-------------->
     70         * offset: -1     0     1     2           n              length -1
     71         * 
     72         *  offset = -1        target < list[0]
     73         *  offset = n         list[n] <= target < list[n+1]
     74         *  offset = length-1  list[length-1] <= target
     75         */
    6776
    6877        /*output: */
    6978        int offset = 0;
    70         int found  = 0; /*found =  0: not found.
    71                                                         found = -1: found, and target is < first element
    72                                                         found =  1: found, and target is == to value at offset
    73                                                         found =  2: found, and target is > to value at offset  and < to value at offset+1
    74                                                         found =  3: found, and target is >= last value */
     79        int found  = 0;
    7580
    7681        /*intermediary: */
    7782        int n0 = 0;
    78         int n1 = int(num_doubles/2);
    79         int n2 = num_doubles-1;
     83        int n1 = int(length/2);
     84        int n2 = length-1;
    8085
    8186        if(target<list[n0]){
    82                 found  = -1;
     87                found  = 1;
    8388                offset = -1;
    8489        }
    8590        else if(target>=list[n2]){
    86                 found  = 3;
    87                 offset = n2-1;
     91                found  = 1;
     92                offset = length-1;
    8893        }
    8994        else{
    90                 while(!found){
     95                for(;;){
    9196                        /*did we find the target?*/
    9297                        if(list[n1]<=target && list[n1+1]>target){
    93                                 found = 1;
     98                                found  = 1;
    9499                                offset = n1;
    95100                                break;
    96101                        }
    97                         if(target < list[n1]){
     102                        else if(target < list[n1]){
    98103                                n2 = n1;
    99104                                n1 = n0 + int((n2-n0)/2);
     
    104109                        }
    105110                }
    106 
    107                 //did we find an exact target?
    108                 if(list[n1]==target) found = 2;
    109111        }
    110112
    111         /*Assign output pointers:*/
     113        /*Assign output pointer and return*/
    112114        *poffset=offset;
    113        
    114         /*Return result: */
    115115        return found;
    116116} /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.