Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/Container/Observations.cpp (revision 13260) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/Container/Observations.cpp (revision 13261) @@ -35,14 +35,14 @@ return; } /*}}}*/ -/*FUNCTION Observations::Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options){{{*/ -Observations::Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options){ +/*FUNCTION Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){{{*/ +Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){ /*Intermediaries*/ int i,j,maxdepth,level,counter,index; int xi,yi; - IssmDouble xmin,xmax,ymin,ymax; - IssmDouble offset,minlength,minspacing,mintrimming,maxtrimming; + IssmPDouble xmin,xmax,ymin,ymax; + IssmPDouble offset,minlength,minspacing,mintrimming,maxtrimming; Observation *observation = NULL; /*Get extrema*/ @@ -65,11 +65,11 @@ if(options->GetOption("boxlength")){ options->Get(&minlength,"boxlength"); if(minlength<=0)_error_("boxlength should be a positive number"); - maxdepth=reCast(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0)); + maxdepth=reCast(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0)); } else{ maxdepth = 30; - minlength=max(xmax-xmin,ymax-ymin)/IssmDouble((1L<(indices); }/*}}}*/ -/*FUNCTION Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata){{{*/ -void Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata){ +/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){{{*/ +void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){ /*Output and Intermediaries*/ bool stop; int nobs,tempnobs,i,j,k,n,counter; - IssmDouble h2,radius2; + IssmPDouble h2,radius2; int *indices = NULL; int *tempindices = NULL; - IssmDouble *dists = NULL; - IssmDouble *x = NULL; - IssmDouble *y = NULL; - IssmDouble *obs = NULL; + IssmPDouble *dists = NULL; + IssmPDouble *x = NULL; + IssmPDouble *y = NULL; + IssmPDouble *obs = NULL; Observation *observation = NULL; /*If radius is not provided or is 0, return all observations*/ @@ -190,7 +190,7 @@ this->quadtree->RangeSearch(&tempindices,&tempnobs,x_interp,y_interp,radius); if(tempnobs){ indices = xNew(tempnobs); - dists = xNew(tempnobs); + dists = xNew(tempnobs); } nobs = 0; for (i=0;i(dists); + xDelete(dists); xDelete(tempindices); if(nobs){ /*Allocate vectors*/ - x = xNew(nobs); - y = xNew(nobs); - obs = xNew(nobs); + x = xNew(nobs); + y = xNew(nobs); + obs = xNew(nobs); /*Loop over all observations and fill in x, y and obs*/ for (i=0;iSize(); if(nobs){ - x = xNew(nobs); - y = xNew(nobs); - obs = xNew(nobs); + x = xNew(nobs); + y = xNew(nobs); + obs = xNew(nobs); for(int i=0;iSize();i++){ observation=(Observation*)this->GetObjectByOffset(i); observation->WriteXYObs(&x[i],&y[i],&obs[i]); @@ -277,15 +277,15 @@ *pnobs=nobs; }/*}}}*/ /*FUNCTION Observations::InterpolationIDW{{{*/ -void Observations::InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power){ +void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){ /*Intermediaries*/ - int i,n_obs; - IssmDouble prediction; - IssmDouble numerator,denominator,h,weight; - IssmDouble *x = NULL; - IssmDouble *y = NULL; - IssmDouble *obs = NULL; + int i,n_obs; + IssmPDouble prediction; + IssmPDouble numerator,denominator,h,weight; + IssmPDouble *x = NULL; + IssmPDouble *y = NULL; + IssmPDouble *obs = NULL; /*Some checks*/ _assert_(maxdata>0); @@ -321,26 +321,26 @@ /*clean-up*/ *pprediction = prediction; - xDelete(x); - xDelete(y); - xDelete(obs); + xDelete(x); + xDelete(y); + xDelete(obs); }/*}}}*/ /*FUNCTION Observations::InterpolationKriging{{{*/ -void Observations::InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram){ +void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){ /*Intermediaries*/ int i,j,n_obs; - IssmDouble prediction,error; - IssmDouble numerator,denominator,ratio; - IssmDouble *x = NULL; - IssmDouble *y = NULL; - IssmDouble *obs = NULL; - IssmDouble *Gamma = NULL; - IssmDouble *GinvG0 = NULL; - IssmDouble *Ginv1 = NULL; - IssmDouble *GinvZ = NULL; - IssmDouble *gamma0 = NULL; - IssmDouble *ones = NULL; + IssmPDouble prediction,error; + IssmPDouble numerator,denominator,ratio; + IssmPDouble *x = NULL; + IssmPDouble *y = NULL; + IssmPDouble *obs = NULL; + IssmPDouble *Gamma = NULL; + IssmPDouble *GinvG0 = NULL; + IssmPDouble *Ginv1 = NULL; + IssmPDouble *GinvZ = NULL; + IssmPDouble *gamma0 = NULL; + IssmPDouble *ones = NULL; /*Some checks*/ _assert_(mindata>0 && maxdata>0); @@ -360,9 +360,9 @@ } /*Allocate intermediary matrix and vectors*/ - Gamma = xNew(n_obs*n_obs); - gamma0 = xNew(n_obs); - ones = xNew(n_obs); + Gamma = xNew(n_obs*n_obs); + gamma0 = xNew(n_obs); + ones = xNew(n_obs); /*First: Create semivariogram matrix for observations*/ for(i=0;i(x); - xDelete(y); - xDelete(obs); - xDelete(Gamma); - xDelete(gamma0); - xDelete(ones); - xDelete(GinvG0); - xDelete(Ginv1); - xDelete(GinvZ); + xDelete(x); + xDelete(y); + xDelete(obs); + xDelete(Gamma); + xDelete(gamma0); + xDelete(ones); + xDelete(GinvG0); + xDelete(Ginv1); + xDelete(GinvZ); }/*}}}*/ /*FUNCTION Observations::InterpolationNearestNeighbor{{{*/ -void Observations::InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius){ +void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){ /*Intermediaries*/ - IssmDouble x,y,obs; + IssmPDouble x,y,obs; /*Get clostest observation*/ this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius); @@ -425,27 +425,27 @@ *pprediction = obs; }/*}}}*/ /*FUNCTION Observations::QuadtreeColoring{{{*/ -void Observations::QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n){ +void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){ int xi,yi,level; for(int i=0;iquadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]); this->quadtree->QuadtreeDepth(&level,xi,yi); - A[i]=(IssmDouble)level; + A[i]=(IssmPDouble)level; } }/*}}}*/ /*FUNCTION Observations::Variomap{{{*/ -void Observations::Variomap(IssmDouble* gamma,IssmDouble *x,int n){ +void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){ /*Output and Intermediaries*/ int i,j,k; - IssmDouble distance; + IssmPDouble distance; Observation *observation1 = NULL; Observation *observation2 = NULL; - IssmDouble *counter = xNew(n); + IssmPDouble *counter = xNew(n); for(j=0;j(counter); + xDelete(counter); }/*}}}*/ Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/modules/Solverx/Solverx.h =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 13260) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 13261) @@ -23,6 +23,7 @@ #endif void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf,Parameters* parameters); +void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n); void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n); #ifdef _HAVE_ADOLC_ Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 13260) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 13261) @@ -51,21 +51,32 @@ #endif }/*}}}*/ +void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/ + + /*Allocate output*/ + double* X = xNew(n); + + /*Solve*/ + SolverxSeq(X,A,B,n); + + /*Assign output pointer*/ + *pX=X; +} +/*}}}*/ + #ifdef _HAVE_ADOLC_ - -int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y) { +int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){/*{{{*/ if(m*(m+1)!=n)_error_("Stiffness matrix should be square!"); SolverxSeq(y,x, x+m*m, m); return 0; -} - -void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){//{{{ +}/*}}}*/ +void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/ // pack inputs to conform to the EDF-prescribed interface - IssmDouble* adoubleEDFin=xNew(n*(n+1)); // packed inputs, i.e. matrix and right hand side - for(int i=0; i(n*(n+1)); // packed inputs, i.e. matrix and right hand side + for(int i=0; i(n*(n+1)); // provide space to transfer inputs during call_ext_fct - IssmPDouble* pdoubleEDFout=xNew(n); // provide space to transfer outputs during call_ext_fct + IssmPDouble* pdoubleEDFout=xNew(n); // provide space to transfer outputs during call_ext_fct // call the wrapped solver through the registry entry we retrieve from parameters call_ext_fct(dynamic_cast * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p, n*(n+1), pdoubleEDFin, adoubleEDFin, @@ -76,8 +87,8 @@ } /*}}}*/ #endif -void SolverxSeq(IssmPDouble * X, IssmPDouble * A, IssmPDouble * B,int n){ //{{{ - #ifdef _HAVE_GSL_ +void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/ +#ifdef _HAVE_GSL_ /*GSL Matrices and vectors: */ int s; gsl_matrix_view a; @@ -108,6 +119,6 @@ xDelete(Acopy); gsl_permutation_free(p); gsl_vector_free(x); - #endif +#endif } /*}}}*/