Changeset 12231


Ignore:
Timestamp:
05/10/12 15:13:50 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added multithreading (VERY efficient for Kriging)

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Container/Observations.cpp

    r12230 r12231  
    113113        if(nobs==0){
    114114                /*No observation found, double range*/
    115                 printf("No observation found within range, doubling range\n");
     115                //printf("No observation found within range, doubling range\n");
    116116                xfree((void**)&indices);
    117117                this->ObservationList(&x,&y,&obs,&nobs,x_interp,y_interp,range*2);
    118118        }
    119119        else{
    120                 if(nobs>1000) printf("Taking more than 1000 observations\n");
     120                //if(nobs>1000) printf("Taking more than 1000 observations\n");
    121121                /*Allocate vectors*/
    122122                x   = (double*)xmalloc(nobs*sizeof(double));

  • issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

    r12230 r12231  
    1616
    1717#include "../../objects/Kriging/GaussianVariogram.h"
     18/*FUNCTION Krigingx{{{*/
    1819int Krigingx(double** ppredictions,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){
    1920
    … …  
    2223
    2324        /*Intermediaries*/
     25        double        range;
     26        char         *output       = NULL;
     27        Variogram    *variogram    = NULL;
     28        Observations *observations = NULL;
     29
     30        /*threading: */
     31        KrigingxThreadStruct gate;
     32        int num=1;
     33#ifdef _MULTITHREADING_
     34        num=_NUMTHREADS_;
     35#endif
     36
     37        /*Get Variogram from Options*/
     38        ProcessVariogram(&variogram,options);
     39        options->Get(&range,"searchrange",0.);
     40
     41        /*Process observation dataset*/
     42        observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
     43
     44        /*Allocate output*/
     45        predictions =(double*)xmalloc(n_interp*sizeof(double));
     46        for(int i=0;i<n_interp;i++) predictions[i]=0;
     47
     48        /*Get output*/
     49        options->Get(&output,"output","prediction");
     50
     51        if(strcmp(output,"quadtree")==0){
     52                observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
     53        }
     54        else if(strcmp(output,"variomap")==0){
     55                observations->Variomap(predictions,x_interp,n_interp);
     56        }
     57        else if(strcmp(output,"prediction")==0){
     58
     59                /*initialize thread parameters: */
     60                gate.n_interp     = n_interp;
     61                gate.x_interp     = x_interp;
     62                gate.y_interp     = y_interp;
     63                gate.range        = range;
     64                gate.variogram    = variogram;
     65                gate.observations = observations;
     66                gate.predictions     = predictions;
     67
     68                /*launch the thread manager with Krigingxt as a core: */
     69                LaunchThread(Krigingxt,(void*)&gate,num);
     70        }
     71        else{
     72                _error_("output '%s' not supported yet",output);
     73        }
     74
     75        /*clean-up and Assign output pointer*/
     76        delete variogram;
     77        delete observations;
     78        xfree((void**)&output);
     79        *ppredictions=predictions;
     80        return 1;
     81}/*}}}*/
     82/*FUNCTION Krigingxt{{{*/
     83void* Krigingxt(void* vpthread_handle){
     84
     85        /*gate variables :*/
     86        KrigingxThreadStruct *gate        = NULL;
     87        pthread_handle       *handle      = NULL;
     88        int my_thread;
     89        int num_threads;
     90        int i0,i1;
     91
     92        /*recover handle and gate: */
     93        handle      = (pthread_handle*)vpthread_handle;
     94        gate        = (KrigingxThreadStruct*)handle->gate;
     95        my_thread   = handle->id;
     96        num_threads = handle->num;
     97
     98        /*recover parameters :*/
     99        int           n_interp     = gate->n_interp;
     100        double       *x_interp     = gate->x_interp;
     101        double       *y_interp     = gate->y_interp;
     102        double        range        = gate->range;
     103        Variogram    *variogram    = gate->variogram;
     104        Observations *observations = gate->observations;
     105        double       *predictions  = gate->predictions;
     106
     107        /*Intermediaries*/
    24108        int           i,j,n_obs;
    25         double        range;
    26109        double        numerator,denominator,ratio;
    27110        double       *x            = NULL;
    … …  
    34117        double       *gamma0       = NULL;
    35118        double       *ones         = NULL;
    36         char         *output       = NULL;
    37         Variogram    *variogram    = NULL;
    38         Observations *observations = NULL;
    39 
    40         /*Get Variogram from Options*/
    41         ProcessVariogram(&variogram,options);
    42         options->Get(&range,"searchrange",0.);
    43 
    44         /*Process observation dataset*/
    45         observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
    46 
    47         /*Allocation output*/
    48         predictions =(double*)xmalloc(n_interp*sizeof(double));
    49         for(i=0;i<n_interp;i++) predictions[i]=0;
    50 
    51         /*Get output*/
    52         options->Get(&output,"output","quadtree");
    53 
    54         if(strcmp(output,"quadtree")==0){
    55                 observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
    56                 *ppredictions=predictions;
    57                 return 1;
    58         }
    59         if(strcmp(output,"variomap")==0){
    60                 observations->Variomap(predictions,x_interp,n_interp);
    61                 *ppredictions=predictions;
    62                 return 1;
    63         }
    64 
    65         /*Loop over all interpolations*/
    66         printf("      interpolation progress:  %5.2lf %%",0.0);
    67         for(int idx=0;idx<n_interp;idx++){
    68                 if(idx%10==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)idx/n_interp*100);
     119
     120        /*partition loop across threads: */
     121        if(my_thread==0) printf("      interpolation progress:  %5.2lf %%",0.0);
     122        PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
     123        for(int idx=i0;idx<i1;idx++){
     124                if(my_thread==0 && idx%10==0) printf("\b\b\b\b\b\b\b%5.2lf %%",double(idx-i0)/double(i1-i0)*100);
    69125
    70126                /*Get list of observations for current point*/
    … …  
    89145
    90146                /*Solve the three linear systems*/
    91                 GslSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
    92                 GslSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
    93                 GslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z
     147                GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
     148                GslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
     149                GslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
    94150
    95151                /*Prepare predictor*/
    … …  
    113169                xfree((void**)&GinvZ);
    114170        }
    115         printf("\b\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
    116 
    117         /*clean-up and Assign output pointer*/
    118         delete variogram;
    119         delete observations;
    120         xfree((void**)&output);
    121         *ppredictions=predictions;
    122         return 1;
    123 }
     171        if(my_thread==0) printf("\b\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
     172
     173        return NULL;
     174}/*}}}*/
    124175
    125176void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/

  • issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h

    r12210 r12231  
    99#include "../../toolkits/toolkits.h"
    1010
     11class Observations;
     12class Variogram;
     13
    1114int  Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
    1215void ProcessVariogram(Variogram **pvariogram,Options* options);
    1316void GslSolve(double** pX,double* A,double* B,int n);
    1417
     18/*threading: */
     19typedef struct{
     20        int           n_interp;
     21        double       *x_interp;
     22        double       *y_interp;
     23        double        range;
     24        Variogram    *variogram;
     25        Observations *observations;
     26        double       *predictions;
     27}KrigingxThreadStruct;
     28
     29void* Krigingxt(void*);
    1530#endif /* _KRIGINGX_H */

  • issm/trunk-jpl/src/c/modules/Shp2Kmlx/Shp2Kmlx.h

    r10380 r12231  
    2222
    2323/* local prototypes: */
    24 int Shp2Kmlx(char* filshp,char* filkml,
    25                          int sgn);
    26 int Shp2Kmlx(char* filshp,char* filkml,
    27                          int sgn,double cm,double sp);
     24int Shp2Kmlx(char* filshp,char* filkml, int sgn);
     25int Shp2Kmlx(char* filshp,char* filkml, int sgn,double cm,double sp);
    2826
    2927#endif  /* _SHP2KMLX_H */
Note: See TracChangeset for help on using the changeset viewer.