Changeset 13261


Ignore:
Timestamp:
09/05/12 07:38:13 (13 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixed Kriging with new solver prototype

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

Legend:

Unmodified
Added
Removed

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

    r13216 r13261  
    3636}
    3737/*}}}*/
    38 /*FUNCTION Observations::Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options){{{*/
    39 Observations::Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options){
     38/*FUNCTION Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){{{*/
     39Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){
    4040
    4141        /*Intermediaries*/
    4242        int          i,j,maxdepth,level,counter,index;
    4343        int          xi,yi;
    44         IssmDouble       xmin,xmax,ymin,ymax;
    45         IssmDouble       offset,minlength,minspacing,mintrimming,maxtrimming;
     44        IssmPDouble       xmin,xmax,ymin,ymax;
     45        IssmPDouble       offset,minlength,minspacing,mintrimming,maxtrimming;
    4646        Observation *observation = NULL;
    4747
    … …  
    6666                options->Get(&minlength,"boxlength");
    6767                if(minlength<=0)_error_("boxlength should be a positive number");
    68                 maxdepth=reCast<int,IssmDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
     68                maxdepth=reCast<int,IssmPDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
    6969        }
    7070        else{
    7171                maxdepth = 30;
    72                 minlength=max(xmax-xmin,ymax-ymin)/IssmDouble((1L<<maxdepth)-1);
     72                minlength=max(xmax-xmin,ymax-ymin)/IssmPDouble((1L<<maxdepth)-1);
    7373        }
    7474
    … …  
    118118/*Methods*/
    119119/*FUNCTION Observations::ClosestObservation{{{*/
    120 void Observations::ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius){
     120void Observations::ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
    121121
    122122        /*Output and Intermediaries*/
    123123        bool         stop;
    124124        int          nobs,i,index;
    125         IssmDouble       h2,hmin2,radius2;
     125        IssmPDouble  h2,hmin2,radius2;
    126126        int         *indices      = NULL;
    127127        Observation *observation  = NULL;
    … …  
    166166
    167167}/*}}}*/
    168 /*FUNCTION Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata){{{*/
    169 void Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata){
     168/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){{{*/
     169void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){
    170170
    171171        /*Output and Intermediaries*/
    172172        bool         stop;
    173173        int          nobs,tempnobs,i,j,k,n,counter;
    174         IssmDouble       h2,radius2;
     174        IssmPDouble  h2,radius2;
    175175        int         *indices      = NULL;
    176176        int         *tempindices  = NULL;
    177         IssmDouble      *dists        = NULL;
    178         IssmDouble      *x            = NULL;
    179         IssmDouble      *y            = NULL;
    180         IssmDouble      *obs          = NULL;
     177        IssmPDouble *dists        = NULL;
     178        IssmPDouble *x            = NULL;
     179        IssmPDouble *y            = NULL;
     180        IssmPDouble *obs          = NULL;
    181181        Observation *observation  = NULL;
    182182
    … …  
    191191        if(tempnobs){
    192192                indices = xNew<int>(tempnobs);
    193                 dists   = xNew<IssmDouble>(tempnobs);
     193                dists   = xNew<IssmPDouble>(tempnobs);
    194194        }
    195195        nobs = 0;
    … …  
    226226                }
    227227        } 
    228         xDelete<IssmDouble>(dists);
     228        xDelete<IssmPDouble>(dists);
    229229        xDelete<int>(tempindices);
    230230
    231231        if(nobs){
    232232                /*Allocate vectors*/
    233                 x   = xNew<IssmDouble>(nobs);
    234                 y   = xNew<IssmDouble>(nobs);
    235                 obs = xNew<IssmDouble>(nobs);
     233                x   = xNew<IssmPDouble>(nobs);
     234                y   = xNew<IssmPDouble>(nobs);
     235                obs = xNew<IssmPDouble>(nobs);
    236236
    237237                /*Loop over all observations and fill in x, y and obs*/
    … …  
    249249        *pnobs=nobs;
    250250}/*}}}*/
    251 /*FUNCTION Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs){{{*/
    252 void Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs){
     251/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){{{*/
     252void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){
    253253
    254254        /*Output and Intermediaries*/
    255255        int          nobs;
    256         IssmDouble      *x            = NULL;
    257         IssmDouble      *y            = NULL;
    258         IssmDouble      *obs          = NULL;
     256        IssmPDouble *x            = NULL;
     257        IssmPDouble *y            = NULL;
     258        IssmPDouble *obs          = NULL;
    259259        Observation *observation  = NULL;
    260260
    … …  
    262262
    263263        if(nobs){
    264                 x   = xNew<IssmDouble>(nobs);
    265                 y   = xNew<IssmDouble>(nobs);
    266                 obs = xNew<IssmDouble>(nobs);
     264                x   = xNew<IssmPDouble>(nobs);
     265                y   = xNew<IssmPDouble>(nobs);
     266                obs = xNew<IssmPDouble>(nobs);
    267267                for(int i=0;i<this->Size();i++){
    268268                        observation=(Observation*)this->GetObjectByOffset(i);
    … …  
    278278}/*}}}*/
    279279/*FUNCTION Observations::InterpolationIDW{{{*/
    280 void Observations::InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power){
     280void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){
    281281
    282282        /*Intermediaries*/
    283         int    i,n_obs;
    284         IssmDouble prediction;
    285         IssmDouble numerator,denominator,h,weight;
    286         IssmDouble *x   = NULL;
    287         IssmDouble *y   = NULL;
    288         IssmDouble *obs = NULL;
     283        int         i,n_obs;
     284        IssmPDouble prediction;
     285        IssmPDouble numerator,denominator,h,weight;
     286        IssmPDouble *x   = NULL;
     287        IssmPDouble *y   = NULL;
     288        IssmPDouble *obs = NULL;
    289289
    290290        /*Some checks*/
    … …  
    322322        /*clean-up*/
    323323        *pprediction = prediction;
    324         xDelete<IssmDouble>(x);
    325         xDelete<IssmDouble>(y);
    326         xDelete<IssmDouble>(obs);
     324        xDelete<IssmPDouble>(x);
     325        xDelete<IssmPDouble>(y);
     326        xDelete<IssmPDouble>(obs);
    327327}/*}}}*/
    328328/*FUNCTION Observations::InterpolationKriging{{{*/
    329 void Observations::InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram){
     329void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){
    330330
    331331        /*Intermediaries*/
    332332        int           i,j,n_obs;
    333         IssmDouble        prediction,error;
    334         IssmDouble        numerator,denominator,ratio;
    335         IssmDouble       *x            = NULL;
    336         IssmDouble       *y            = NULL;
    337         IssmDouble       *obs          = NULL;
    338         IssmDouble       *Gamma        = NULL;
    339         IssmDouble       *GinvG0       = NULL;
    340         IssmDouble       *Ginv1        = NULL;
    341         IssmDouble       *GinvZ        = NULL;
    342         IssmDouble       *gamma0       = NULL;
    343         IssmDouble       *ones         = NULL;
     333        IssmPDouble   prediction,error;
     334        IssmPDouble   numerator,denominator,ratio;
     335        IssmPDouble  *x            = NULL;
     336        IssmPDouble  *y            = NULL;
     337        IssmPDouble  *obs          = NULL;
     338        IssmPDouble  *Gamma        = NULL;
     339        IssmPDouble  *GinvG0       = NULL;
     340        IssmPDouble  *Ginv1        = NULL;
     341        IssmPDouble  *GinvZ        = NULL;
     342        IssmPDouble  *gamma0       = NULL;
     343        IssmPDouble  *ones         = NULL;
    344344
    345345        /*Some checks*/
    … …  
    361361
    362362        /*Allocate intermediary matrix and vectors*/
    363         Gamma  = xNew<IssmDouble>(n_obs*n_obs);
    364         gamma0 = xNew<IssmDouble>(n_obs);
    365         ones   = xNew<IssmDouble>(n_obs);
     363        Gamma  = xNew<IssmPDouble>(n_obs*n_obs);
     364        gamma0 = xNew<IssmPDouble>(n_obs);
     365        ones   = xNew<IssmPDouble>(n_obs);
    366366
    367367        /*First: Create semivariogram matrix for observations*/
    … …  
    402402        *pprediction = prediction;
    403403        *perror = error;
    404         xDelete<IssmDouble>(x);
    405         xDelete<IssmDouble>(y);
    406         xDelete<IssmDouble>(obs);
    407         xDelete<IssmDouble>(Gamma);
    408         xDelete<IssmDouble>(gamma0);
    409         xDelete<IssmDouble>(ones);
    410         xDelete<IssmDouble>(GinvG0);
    411         xDelete<IssmDouble>(Ginv1);
    412         xDelete<IssmDouble>(GinvZ);
     404        xDelete<IssmPDouble>(x);
     405        xDelete<IssmPDouble>(y);
     406        xDelete<IssmPDouble>(obs);
     407        xDelete<IssmPDouble>(Gamma);
     408        xDelete<IssmPDouble>(gamma0);
     409        xDelete<IssmPDouble>(ones);
     410        xDelete<IssmPDouble>(GinvG0);
     411        xDelete<IssmPDouble>(Ginv1);
     412        xDelete<IssmPDouble>(GinvZ);
    413413
    414414}/*}}}*/
    415415/*FUNCTION Observations::InterpolationNearestNeighbor{{{*/
    416 void Observations::InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius){
     416void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
    417417
    418418        /*Intermediaries*/
    419         IssmDouble        x,y,obs;
     419        IssmPDouble x,y,obs;
    420420
    421421        /*Get clostest observation*/
    … …  
    426426}/*}}}*/
    427427/*FUNCTION Observations::QuadtreeColoring{{{*/
    428 void Observations::QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n){
     428void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
    429429
    430430        int xi,yi,level;
    … …  
    433433                this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
    434434                this->quadtree->QuadtreeDepth(&level,xi,yi);
    435                 A[i]=(IssmDouble)level;
     435                A[i]=(IssmPDouble)level;
    436436        }
    437437
    438438}/*}}}*/
    439439/*FUNCTION Observations::Variomap{{{*/
    440 void Observations::Variomap(IssmDouble* gamma,IssmDouble *x,int n){
     440void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){
    441441
    442442        /*Output and Intermediaries*/
    443443        int          i,j,k;
    444         IssmDouble       distance;
     444        IssmPDouble  distance;
    445445        Observation *observation1 = NULL;
    446446        Observation *observation2 = NULL;
    447447
    448         IssmDouble *counter = xNew<IssmDouble>(n);
     448        IssmPDouble *counter = xNew<IssmPDouble>(n);
    449449        for(j=0;j<n;j++) counter[j] = 0.0;
    450450        for(j=0;j<n;j++) gamma[j]   = 0.0;
    … …  
    475475
    476476        /*Assign output pointer*/
    477         xDelete<IssmDouble>(counter);
    478 }/*}}}*/
     477        xDelete<IssmPDouble>(counter);
     478}/*}}}*/

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

    r13216 r13261  
    2424
    2525void SolverxSeq(SeqVec<IssmDouble>** puf,SeqMat<IssmDouble>* Kff, SeqVec<IssmDouble>* pf,Parameters* parameters);
     26void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
    2627void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
    2728

  • issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp

    r13216 r13261  
    5252
    5353}/*}}}*/
     54void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
     55
     56        /*Allocate output*/
     57        double* X = xNew<double>(n);
     58
     59        /*Solve*/
     60        SolverxSeq(X,A,B,n);
     61
     62        /*Assign output pointer*/
     63        *pX=X;
     64}
     65/*}}}*/
     66
    5467#ifdef _HAVE_ADOLC_
    55 
    56 int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y) {
     68int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){/*{{{*/
    5769    if(m*(m+1)!=n)_error_("Stiffness matrix should be square!");
    5870    SolverxSeq(y,x, x+m*m, m);
    5971    return 0;
    60 }
    61 
    62 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){//{{{
     72}/*}}}*/
     73void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
    6374        // pack inputs to conform to the EDF-prescribed interface
    64         IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
    65         for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i]; // pack matrix
    66         for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
     75        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1));  // packed inputs, i.e. matrix and right hand side
     76        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i];      // pack matrix
     77        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i];      // pack the right hand side
    6778        IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
    68         IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);       // provide space to transfer outputs during call_ext_fct
     79        IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);           // provide space to transfer outputs during call_ext_fct
    6980        // call the wrapped solver through the registry entry we retrieve from parameters
    7081        call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
    … …  
    7788/*}}}*/
    7889#endif
    79 void SolverxSeq(IssmPDouble * X, IssmPDouble * A, IssmPDouble * B,int n){ //{{{
    80         #ifdef _HAVE_GSL_
     90void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
     91#ifdef _HAVE_GSL_
    8192        /*GSL Matrices and vectors: */
    8293        int              s;
    … …  
    109120        gsl_permutation_free(p);
    110121        gsl_vector_free(x);
    111         #endif
     122#endif
    112123}
    113124/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.