Index: /issm/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12376) +++ /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12377) @@ -21,4 +21,5 @@ #include "../include/include.h" #include "../EnumDefinitions/EnumDefinitions.h" +#include "../io/io.h" using namespace std; @@ -71,5 +72,5 @@ /*Initialize Quadtree*/ - printf("Generating quadtree with a maximum box size %g (depth=%i)... ",minlength,maxdepth); + _printf_(true,"Generating quadtree with a maximum box size %g (depth=%i)... ",minlength,maxdepth); this->quadtree = new Quadtree(xmin,xmax,ymin,ymax,maxdepth); @@ -101,7 +102,7 @@ } } - printf("done\n"); - printf("Initial number of observations: %i\n",n); - printf(" Final number of observations: %i\n",this->quadtree->NbObs); + _printf_(true,"done\n"); + _printf_(true,"Initial number of observations: %i\n",n); + _printf_(true," Final number of observations: %i\n",this->quadtree->NbObs); } /*}}}*/ Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12376) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12377) @@ -633,4 +633,23 @@ ./modules/Krigingx/Krigingx.h +#For parallel kriging, only difference is ./modules/Krigingx/pKrigingx.cpp with no multithreading +pkriging_sources = ./Container/Observations.h\ + ./Container/Observations.cpp\ + ./objects/Kriging/Variogram.h \ + ./objects/Kriging/GaussianVariogram.h\ + ./objects/Kriging/GaussianVariogram.cpp\ + ./objects/Kriging/ExponentialVariogram.h\ + ./objects/Kriging/ExponentialVariogram.cpp\ + ./objects/Kriging/SphericalVariogram.h\ + ./objects/Kriging/SphericalVariogram.cpp\ + ./objects/Kriging/PowerVariogram.h\ + ./objects/Kriging/PowerVariogram.cpp\ + ./objects/Kriging/Quadtree.h\ + ./objects/Kriging/Quadtree.cpp\ + ./objects/Kriging/Observation.h\ + ./objects/Kriging/Observation.cpp\ + ./modules/Krigingx/pKrigingx.cpp\ + ./modules/Krigingx/Krigingx.h + #}}} #Kml sources {{{1 @@ -964,4 +983,5 @@ libISSMCore_a_SOURCES = $(issm_sources) +libISSMCore_a_SOURCES += $(pkriging_sources) libISSMCore_a_CXXFLAGS = $(ALLCXXFLAGS) @@ -1005,4 +1025,8 @@ issm_SOURCES = solutions/issm.cpp issm_CXXFLAGS= -fPIC $(CXXFLAGS) $(CXXOPTFLAGS) $(COPTFLAGS) + +bin_PROGRAMS += kriging +kriging_SOURCES = solutions/kriging.cpp +kriging_CXXFLAGS= -fPIC $(CXXFLAGS) $(CXXOPTFLAGS) $(COPTFLAGS) #}}} #Automatic differentiation: append this fold to the end of the src/c/Makefile.am to get this Makefile.am {{{ Index: /issm/trunk-jpl/src/c/io/Disk/diskio.h =================================================================== --- /issm/trunk-jpl/src/c/io/Disk/diskio.h (revision 12376) +++ /issm/trunk-jpl/src/c/io/Disk/diskio.h (revision 12377) @@ -10,7 +10,4 @@ #include "../../include/include.h" -class DataSet; -class Parameters; - FILE* pfopen(char* filename,const char* format); void pfclose(FILE* fid,char* filename); Index: /issm/trunk-jpl/src/c/io/io.h =================================================================== --- /issm/trunk-jpl/src/c/io/io.h (revision 12376) +++ /issm/trunk-jpl/src/c/io/io.h (revision 12377) @@ -6,11 +6,9 @@ #define _ISSM_IO_H_ -#ifdef HAVE_CONFIG_H //config.h {{{ +#ifdef HAVE_CONFIG_H #include #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif -//}}} - #include "./Disk/diskio.h" Index: /issm/trunk-jpl/src/c/issm.h =================================================================== --- /issm/trunk-jpl/src/c/issm.h (revision 12376) +++ /issm/trunk-jpl/src/c/issm.h (revision 12377) @@ -22,4 +22,3 @@ #include "./modules/modules.h" - #endif //ifndef _ISSM_H_ Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12376) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12377) @@ -10,10 +10,7 @@ #include "../../Container/Observations.h" #include "../modules.h" - #ifdef _HAVE_GSL_ #include #endif - -#include "../../objects/Kriging/GaussianVariogram.h" /*FUNCTION Krigingx{{{*/ int Krigingx(double** ppredictions,double **perror,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){ Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12376) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12377) @@ -13,4 +13,5 @@ int Krigingx(double** ppredictions,double **perror,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options); +int pKrigingx(double** ppredictions,double **perror,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options); void ProcessVariogram(Variogram **pvariogram,Options* options); void GslSolve(double** pX,double* A,double* B,int n); Index: /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp (revision 12377) +++ /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp (revision 12377) @@ -0,0 +1,212 @@ +/*!\file: Kriging.cpp + * \brief "c" core code for Kriging + */ + +#include "./Krigingx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../objects/objects.h" +#include "../../Container/Observations.h" +#include "../modules.h" + +#ifdef _HAVE_GSL_ +#include +#endif + +/*FUNCTION pKrigingx{{{*/ +int pKrigingx(double** ppredictions,double **perror,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){ + + extern int num_procs; + extern int my_rank; + + /*output*/ + double *predictions = NULL; + double *error = NULL; + + /*Intermediaries*/ + int mindata,maxdata; + double radius; + char *output = NULL; + Variogram *variogram = NULL; + Observations *observations = NULL; + + /*Get Variogram from Options*/ + ProcessVariogram(&variogram,options); + options->Get(&radius,"searchradius",0.); + options->Get(&mindata,"mindata",1); + options->Get(&maxdata,"maxdata",50); + + /*Process observation dataset*/ + observations=new Observations(obs_list,obs_x,obs_y,obs_length,options); + + /*Allocate output*/ + predictions =(double*)xcalloc(n_interp,sizeof(double)); + error =(double*)xcalloc(n_interp,sizeof(double)); + + /*Get output*/ + options->Get(&output,"output","prediction"); + + if(strcmp(output,"quadtree")==0){ + observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp); + } + else if(strcmp(output,"variomap")==0){ + observations->Variomap(predictions,x_interp,n_interp); + } + else if(strcmp(output,"prediction")==0){ + + /*Intermediaries*/ + int i,j,n_obs; + double numerator,denominator,ratio,localpercent; + double *x = NULL; + double *y = NULL; + double *obs = NULL; + double *Gamma = NULL; + double *GinvG0 = NULL; + double *Ginv1 = NULL; + double *GinvZ = NULL; + double *gamma0 = NULL; + double *ones = NULL; + + /*partition loop across threads: */ + for(int idx=my_rank;idxObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata); + if(n_obsSemiVariogram(x[i]-x[j],y[i]-y[j]); + Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]); + Gamma[j*n_obs+i] = Gamma[i*n_obs+j]; + } + } + for(i=0;iSemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]); + for(i=0;iCovariance(x[i]-x_interp[idx],y[i]-y_interp[idx]); + + /*Solve the three linear systems*/ + GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0 + GslSolve(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones + GslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z + + /*Prepare predictor*/ + numerator=-1.; denominator=0.; + for(i=0;iGetOption("model")){ + options->Get(&model,"model"); + if (strcmp(model,"gaussian")==0) variogram = new GaussianVariogram(options); + else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options); + else if(strcmp(model,"spherical")==0) variogram = new SphericalVariogram(options); + else if(strcmp(model,"power")==0) variogram = new PowerVariogram(options); + else _error_("variogram %s not supported yet (list of supported variogram: gaussian, exponential, spherical and power)",model); + } + else variogram = new GaussianVariogram(options); + + /*Assign output pointer*/ + xfree((void**)&model); + *pvariogram = variogram; +}/*}}}*/ +void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/ +#ifdef _HAVE_GSL_ + + /*GSL Matrices and vectors: */ + int s; + gsl_matrix_view a; + gsl_vector_view b; + gsl_vector *x = NULL; + gsl_permutation *p = NULL; + + /*A will be modified by LU decomposition. Use copy*/ + double* Acopy = (double*)xmalloc(n*n*sizeof(double)); + memcpy(Acopy,A,n*n*sizeof(double)); + + /*Initialize gsl matrices and vectors: */ + a = gsl_matrix_view_array (Acopy,n,n); + b = gsl_vector_view_array (B,n); + x = gsl_vector_alloc (n); + + /*Run LU and solve: */ + p = gsl_permutation_alloc (n); + gsl_linalg_LU_decomp (&a.matrix, p, &s); + gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x); + + //printf ("x = \n"); + //gsl_vector_fprintf (stdout, x, "%g"); + + /*Copy result*/ + double* X = (double*)xmalloc(n*sizeof(double)); + memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double)); + + /*Clean up and assign output pointer*/ + xfree((void**)&Acopy); + gsl_permutation_free(p); + gsl_vector_free(x); + *pX=X; +#else + _error_("GSL support required"); +#endif + }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 12376) +++ /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 12377) @@ -16,5 +16,5 @@ #include "../../objects/objects.h" -void OutputResultsx( Elements* elements, Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Results* results){ +void OutputResultsx(Elements* elements, Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Results* results){ extern int my_rank; Index: /issm/trunk-jpl/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/IoModel.cpp (revision 12376) +++ /issm/trunk-jpl/src/c/objects/IoModel.cpp (revision 12377) @@ -71,11 +71,13 @@ IoModel::~IoModel(){ - delete this->constants; + if(this->constants) delete this->constants; /*Some checks in debugging mode*/ #ifdef _ISSM_DEBUG_ - for(int i=0;idata[i]){ - _printf_("Info: previous pointer of %s has not been freed (DeleteData has not been called)",EnumToStringx(i)); + if(this->data){ + for(int i=0;idata[i]){ + _printf_("Info: previous pointer of %s has not been freed (DeleteData has not been called)",EnumToStringx(i)); + } } } @@ -83,5 +85,4 @@ xfree((void**)&this->data); - xfree((void**)&this->my_elements); xfree((void**)&this->my_nodes); @@ -98,8 +99,6 @@ int record_enum = 0; - /*Check that some fields have been allocated*/ _assert_(this->fid || my_rank); - /*Go find in the binary file, the position of the data we want to fetch: */ Index: /issm/trunk-jpl/src/c/objects/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/objects/IoModel.h (revision 12376) +++ /issm/trunk-jpl/src/c/objects/IoModel.h (revision 12377) @@ -18,5 +18,4 @@ private: - FILE *fid; //pointer to input file IssmDouble **data; //this dataset holds temporary data, memory intensive. Parameters *constants; //this dataset holds all IssmDouble, int, bool and char *parameters read in from the input file.* @@ -24,4 +23,5 @@ public: /*This data needs to stay memory resident at all time, even if it's memory intensive: */ + FILE *fid; //pointer to input file bool *my_elements; bool *my_nodes; Index: /issm/trunk-jpl/src/c/solutions/kriging.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/kriging.cpp (revision 12377) +++ /issm/trunk-jpl/src/c/solutions/kriging.cpp (revision 12377) @@ -0,0 +1,257 @@ +/*!\file: kriging.cpp + * \brief: kriging main parallel program + */ + +#include "../issm.h" +#include "../include/globals.h" + +/*Local prototypes*/ +void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,int argc,char **argv); +void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid); + +int main(int argc,char **argv){ + + /*I/O: */ + FILE *output_fid = NULL; + FILE *input_fid = NULL; + bool waitonlock = false; + + /*File names*/ + char *lockfilename = NULL; + char *binfilename = NULL; + char *outbinfilename = NULL; + + /*Input*/ + int ninterp,nobs; + double *x = NULL; + double *y = NULL; + double *data = NULL; + double *x_interp = NULL; + double *y_interp = NULL; + Options *options = NULL; + + /*Output*/ + double *predictions = NULL; + double *error = NULL; + + ISSMBOOT(); + + /*Initialize environments: Petsc, MPI, etc...: */ +#ifdef _HAVE_PETSC_ + int ierr=PetscInitialize(&argc,&argv,(char*)0,""); + if(ierr) _error_("Could not initialize Petsc"); +#else +#ifdef _HAVE_MPI_ + MPI_Init(&argc,&argv); +#endif +#endif + + /*Size and rank: */ +#ifdef _HAVE_MPI_ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); +#endif + + /*First process inputs*/ + _printf_(true,"\n"); + _printf_(true,"Ice Sheet System Model (%s) version %s\n",PACKAGE_NAME,PACKAGE_VERSION); + _printf_(true,"(website: %s contact: %s)\n",PACKAGE_URL,PACKAGE_BUGREPORT); + _printf_(true,"\n"); + ProcessArguments2(&binfilename,&outbinfilename,&lockfilename,argc,argv); + + /*Process input files*/ + input_fid=pfopen(binfilename,"rb"); + ProcessInputfile(&x,&y,&data,&nobs,&x_interp,&y_interp,&ninterp,&options,input_fid); + pfclose(input_fid,binfilename); + + _printf_(true,"call computational core:\n"); + pKrigingx(&predictions,&error,x,y,data,nobs,x_interp,y_interp,ninterp,options); + + _printf_(true,"write results to disk:\n"); + output_fid=pfopen(outbinfilename,"wb"); + Results *results = new Results(); + results->AddObject(new DoubleVecExternalResult(results->Size()+1,0,predictions,ninterp,1,0)); + results->AddObject(new DoubleVecExternalResult(results->Size()+1,0,error,ninterp,1,0)); + for(int i=0;iSize();i++){ + ExternalResult* result=(ExternalResult*)results->GetObjectByOffset(i); + result->WriteData(output_fid,1); + } + pfclose(output_fid,outbinfilename); + + /*Close output and petsc options file and write lock file if requested*/ + _printf_(true,"write lock file:\n"); + WriteLockFile(lockfilename); + + /*Free ressources */ + xfree((void**)&lockfilename); + xfree((void**)&binfilename); + xfree((void**)&outbinfilename); + xfree((void**)&x); + xfree((void**)&y); + xfree((void**)&data); + xfree((void**)&x_interp); + xfree((void**)&y_interp); + xfree((void**)&predictions); + xfree((void**)&error); + delete options; + delete results; + +#ifdef _HAVE_PETSC_ + _printf_(true,"closing MPI and Petsc\n"); + PetscFinalize(); +#else +#ifdef _HAVE_MPI_ + _printf_(true,"closing MPI and Petsc\n"); + MPI_Finalize(); +#endif +#endif + + /*end module: */ + ISSMEND(); + + return 0; //unix success return; +} + +void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,int argc,char **argv){ + + char *modelname = NULL; + char *binfilename = NULL; + char *outbinfilename = NULL; + char *lockfilename = NULL; + + if(argc<2)_error_("Usage error: missing model name"); + modelname=argv[2]; + binfilename = (char*)xmalloc((strlen(modelname)+strlen(".bin") +1)*sizeof(char)); sprintf(binfilename, "%s%s",modelname,".bin"); + outbinfilename = (char*)xmalloc((strlen(modelname)+strlen(".outbin")+1)*sizeof(char)); sprintf(outbinfilename,"%s%s",modelname,".outbin"); + lockfilename = (char*)xmalloc((strlen(modelname)+strlen(".lock") +1)*sizeof(char)); sprintf(lockfilename, "%s%s",modelname,".lock"); + + /*Clean up and assign output pointer*/ + *pbinfilename=binfilename; + *poutbinfilename=outbinfilename; + *plockfilename=lockfilename; +} + +void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){ + + int ninterp,nobs; + double *x = NULL; + double *y = NULL; + double *data = NULL; + double *x_interp = NULL; + double *y_interp = NULL; + Options *options = NULL; + + int dummy,M,N; + IoModel* iomodel = new IoModel(); + iomodel->fid=fid; + iomodel->CheckEnumSync(); + + /*Read x*/ + fid=iomodel->SetFilePointerToData(&dummy,&dummy,0); + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); + if(M*N){ + x=(double*)xmalloc(M*N*sizeof(double)); + if(my_rank==0){ + if(fread(x,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix "); + } + MPI_Bcast(x,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); + } + nobs=M*N; + + /*Read y*/ + fid=iomodel->SetFilePointerToData(&dummy,&dummy,1); + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); + if(M*N){ + y=(double*)xmalloc(M*N*sizeof(double)); + if(my_rank==0){ + if(fread(y,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix "); + } + MPI_Bcast(y,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); + } + _assert_(M*N==nobs); + + /*Read data*/ + fid=iomodel->SetFilePointerToData(&dummy,&dummy,2); + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); + if(M*N){ + data=(double*)xmalloc(M*N*sizeof(double)); + if(my_rank==0){ + if(fread(data,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix "); + } + MPI_Bcast(data,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); + } + _assert_(M*N==nobs); + + /*Read x_interp*/ + fid=iomodel->SetFilePointerToData(&dummy,&dummy,3); + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); + if(M*N){ + x_interp=(double*)xmalloc(M*N*sizeof(double)); + if(my_rank==0){ + if(fread(x_interp,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix "); + } + MPI_Bcast(x_interp,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); + } + ninterp=N*M; + + /*Read y_interp*/ + fid=iomodel->SetFilePointerToData(&dummy,&dummy,4); + if(my_rank==0){ + if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); + } + MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); + if(my_rank==0){ + if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); + } + MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); + if(M*N){ + y_interp=(double*)xmalloc(M*N*sizeof(double)); + if(my_rank==0){ + if(fread(y_interp,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix "); + } + MPI_Bcast(y_interp,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); + } + _assert_(N*M==ninterp); + + /*Read options*/ + options = new Options(); + + /*Assign output pointer*/ + *px = x; + *py = y; + *pdata = data; + *pnobs = nobs; + *px_interp = x_interp; + *py_interp = y_interp; + *pninterp = ninterp; + *poptions = options; + delete iomodel; +}