Changeset 25425
- Timestamp:
- 08/18/20 08:18:44 (5 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/IoModel.cpp
r25379 r25425 1263 1263 } 1264 1264 /*}}}*/ 1265 #ifdef _HAVE_AD_ 1266 void IoModel::FetchData(IssmPDouble** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/ 1267 1268 /*First, look if has already been loaded (might be an independent variable)*/ 1269 vector<IoData*>::iterator iter; 1270 for(iter=data.begin();iter<data.end();iter++){ 1271 IoData* iodata=*iter; 1272 if(strcmp(iodata->name,data_name)==0){ 1273 _error_(data_name <<" has already been loaded as in IssmDouble and cannot be converted to IssmPDouble"); 1274 } 1275 } 1276 1277 /*output: */ 1278 int M,N; 1279 IssmPDouble *matrix = NULL; 1280 int code = 0; 1281 1282 /*recover my_rank:*/ 1283 int my_rank=IssmComm::GetRank(); 1284 1285 /*Set file pointer to beginning of the data: */ 1286 fid=this->SetFilePointerToData(&code,NULL,data_name); 1287 if(code!=5 && code!=6 && code!=7 && code!=10)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")"); 1288 1289 /*Now fetch: */ 1290 1291 /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ 1292 /*numberofelements: */ 1293 if(my_rank==0){ 1294 if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); 1295 } 1296 ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 1297 1298 if(my_rank==0){ 1299 if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); 1300 } 1301 ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 1302 1303 /*Now allocate matrix: */ 1304 if(M*N){ 1305 if(code==10){ 1306 /*Special case for Compressed mat*/ 1307 IssmPDouble offset,range; 1308 if(my_rank==0) if(fread(&offset,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read offset"); 1309 ISSM_MPI_Bcast(&offset,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 1310 1311 if(my_rank==0) if(fread(&range,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read range"); 1312 ISSM_MPI_Bcast(&range,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 1313 1314 *pmatrix=xNew<IssmPDouble>(M*N); 1315 1316 /*Read matrix*/ 1317 uint8_t* rawmatrix=xNew<uint8_t>((M-1)*N); 1318 if(my_rank==0) if(fread(rawmatrix,(M-1)*N*sizeof(char),1,fid)!=1) _error_("could not read matrix "); 1319 ISSM_MPI_Bcast(rawmatrix,(M-1)*N,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 1320 1321 for(int i=0;i<(M-1)*N;++i) (*pmatrix)[i]=offset+range*reCast<IssmPDouble>(rawmatrix[i])/255.; 1322 xDelete<uint8_t>(rawmatrix); 1323 1324 /*read time now*/ 1325 IssmPDouble* timematrix=xNew<IssmPDouble>(N); 1326 if(my_rank==0) if(fread(timematrix,N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read time in compressed matrix"); 1327 ISSM_MPI_Bcast(timematrix,N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 1328 1329 for(int i=0;i<N;++i) (*pmatrix)[(M-1)*N+i]=timematrix[i]; 1330 xDelete<IssmPDouble>(timematrix); 1331 1332 } 1333 else{ 1334 /*Read matrix on node 0, then broadcast: */ 1335 matrix=xNew<IssmPDouble>(M*N); 1336 if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); 1337 ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 1338 *pmatrix=matrix; 1339 } 1340 } 1341 else{ 1342 *pmatrix=NULL; 1343 } 1344 /*Assign output pointers: */ 1345 if(pM) *pM=M; 1346 if(pN) *pN=N; 1347 } 1348 /*}}}*/ 1349 #endif 1265 1350 void IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,const char* data_name){/*{{{*/ 1266 1351 -
issm/trunk-jpl/src/c/classes/IoModel.h
r25379 r25425 136 136 void FetchData(int** pmatrix,int* pM,int* pN,const char* data_name); 137 137 void FetchData(IssmDouble** pscalarmatrix,int* pM,int* pN,const char* data_name); 138 #if _HAVE_AD_ && !defined(_WRAPPERS_) 139 void FetchData(IssmPDouble** pscalarmatrix,int* pM,int* pN,const char* data_name); 140 #endif 138 141 void FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name); 139 142 void FetchData(Options *options,const char* data_name); -
issm/trunk-jpl/src/c/classes/kriging/Observations.cpp
r24506 r25425 44 44 45 45 /*Get tree type (FIXME)*/ 46 IssmDouble dtree = 0.;46 double dtree = 0.; 47 47 options->Get(&dtree,"treetype",1.); 48 48 this->treetype = reCast<int>(dtree); … … 552 552 IssmPDouble* B = xNew<IssmPDouble>(n_obs+1); 553 553 554 Issm Double unbias = variogram->Covariance(0.,0.);554 IssmPDouble unbias = variogram->Covariance(0.,0.); 555 555 /*First: Create semivariogram matrix for observations*/ 556 556 for(i=0;i<n_obs;i++){ -
issm/trunk-jpl/src/c/classes/kriging/Observations.h
r20827 r25425 25 25 /*constructors, destructors*/ 26 26 Observations(); 27 Observations(Issm Double* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);27 Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options); 28 28 ~Observations(); 29 29 30 30 /*Initialize data structures*/ 31 void InitCovertree(Issm Double* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);32 void InitQuadtree(Issm Double* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);31 void InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options); 32 void InitQuadtree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options); 33 33 34 34 /*Methods*/ 35 void ClosestObservation(Issm Double *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);36 void ClosestObservationCovertree(Issm Double *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);37 void ClosestObservationQuadtree(Issm Double *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);35 void ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius); 36 void ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius); 37 void ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius); 38 38 void Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius); 39 void InterpolationIDW(Issm Double *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);40 void InterpolationV4(Issm Double *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);41 void InterpolationKriging(Issm Double *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);42 void InterpolationNearestNeighbor(Issm Double *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);43 void ObservationList(Issm Double **px,IssmDouble **py,IssmDouble **pobs,int* pnobs);44 void ObservationList(Issm Double **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);45 void ObservationListCovertree(Issm Double **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);46 void ObservationListQuadtree(Issm Double **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);47 void QuadtreeColoring(Issm Double* A,IssmDouble *x,IssmDouble *y,int n);48 void Variomap(Issm Double* gamma,IssmDouble *x,int n);39 void InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power); 40 void InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata); 41 void InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram); 42 void InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius); 43 void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs); 44 void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata); 45 void ObservationListCovertree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata); 46 void ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata); 47 void QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n); 48 void Variomap(IssmPDouble* gamma,IssmPDouble *x,int n); 49 49 50 50 }; -
issm/trunk-jpl/src/c/main/kriging.cpp
r21068 r25425 7 7 /*Local prototypes*/ 8 8 void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv); 9 void ProcessInputfile( IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid);9 void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid); 10 10 11 11 int main(int argc,char **argv){ … … 23 23 /*Input*/ 24 24 int ninterp,nobs; 25 IssmDouble *x = NULL;26 IssmDouble *y = NULL;27 IssmDouble *data = NULL;28 IssmDouble *x_interp = NULL;29 IssmDouble *y_interp = NULL;25 double *x = NULL; 26 double *y = NULL; 27 double *data = NULL; 28 double *x_interp = NULL; 29 double *y_interp = NULL; 30 30 Options *options = NULL; 31 31 32 32 /*Output*/ 33 IssmDouble *predictions = NULL;34 IssmDouble *error = NULL;33 double *predictions = NULL; 34 double *error = NULL; 35 35 36 36 /*Initialize exception trapping: */ … … 73 73 xDelete<char>(outbinfilename); 74 74 xDelete<char>(rootpath); 75 xDelete< IssmDouble>(x);76 xDelete< IssmDouble>(y);77 xDelete< IssmDouble>(data);78 xDelete< IssmDouble>(x_interp);79 xDelete< IssmDouble>(y_interp);80 xDelete< IssmDouble>(predictions);81 xDelete< IssmDouble>(error);75 xDelete<double>(x); 76 xDelete<double>(y); 77 xDelete<double>(data); 78 xDelete<double>(x_interp); 79 xDelete<double>(y_interp); 80 xDelete<double>(predictions); 81 xDelete<double>(error); 82 82 delete options; 83 83 delete results; … … 130 130 } 131 131 132 void ProcessInputfile( IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid){132 void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){ 133 133 134 int ninterp,nobs;135 IssmDouble *x = NULL;136 IssmDouble *y = NULL;137 IssmDouble *data = NULL;138 IssmDouble *x_interp = NULL;139 IssmDouble *y_interp = NULL;134 int ninterp,nobs; 135 double *x = NULL; 136 double *y = NULL; 137 double *data = NULL; 138 double *x_interp = NULL; 139 double *y_interp = NULL; 140 140 Options *options = NULL; 141 141 … … 144 144 iomodel->fid=fid; 145 145 iomodel->CheckFile(); 146 iomodel->FetchData(&x,&M,&N,"md.x"); nobs=M*N;147 iomodel->FetchData(&y,&M,&N,"md.y"); _assert_(M*N==nobs);148 iomodel->FetchData(&data,&M,&N,"md.data"); _assert_(M*N==nobs);146 iomodel->FetchData(&x,&M,&N,"md.x"); nobs=M*N; 147 iomodel->FetchData(&y,&M,&N,"md.y"); _assert_(M*N==nobs); 148 iomodel->FetchData(&data,&M,&N,"md.data"); _assert_(M*N==nobs); 149 149 iomodel->FetchData(&x_interp,&M,&N,"md.x_interp"); ninterp=M*N; 150 150 iomodel->FetchData(&y_interp,&M,&N,"md.y_interp"); _assert_(M*N==ninterp);
Note:
See TracChangeset
for help on using the changeset viewer.