source: issm/trunk/src/c/main/kriging.cpp@ 25836

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

merged trunk-jpl and trunk for revision 25834
File size: 5.4 KB
Line 
1/*!\file: kriging.cpp
2 * \brief: kriging main parallel program
3 */
4
5#include "./issm.h"
6
7/*Local prototypes*/
8void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv);
9void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid);
10
11int main(int argc,char **argv){
12
13 /*I/O: */
14 FILE *output_fid = NULL;
15 FILE *input_fid = NULL;
16
17 /*File names*/
18 char *lockfilename = NULL;
19 char *binfilename = NULL;
20 char *outbinfilename = NULL;
21 char *rootpath = NULL;
22
23 /*Input*/
24 int ninterp,nobs;
25 double *x = NULL;
26 double *y = NULL;
27 double *data = NULL;
28 double *x_interp = NULL;
29 double *y_interp = NULL;
30 Options *options = NULL;
31
32 /*Output*/
33 double *predictions = NULL;
34 double *error = NULL;
35
36 /*Initialize exception trapping: */
37 ExceptionTrapBegin();
38
39 /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
40 ISSM_MPI_Comm comm=EnvironmentInit(argc,argv);
41 IssmComm::SetComm(comm);
42
43 ProcessArguments2(&binfilename,&outbinfilename,&lockfilename,&rootpath,argc,argv);
44
45 /*Process input files*/
46 input_fid=pfopen(binfilename,"rb");
47 ProcessInputfile(&x,&y,&data,&nobs,&x_interp,&y_interp,&ninterp,&options,input_fid);
48 pfclose(input_fid,binfilename);
49
50 _printf0_("call computational core:\n");
51 pKrigingx(&predictions,&error,x,y,data,nobs,x_interp,y_interp,ninterp,options);
52
53 _printf0_("write results to disk:\n");
54 Results *results = new Results();
55 if(IssmComm::GetRank()==0){
56 output_fid=pfopen0(outbinfilename,"wb");
57 results->AddObject(new GenericExternalResult<double*>(results->Size()+1,"predictions",predictions,ninterp,1,1,0));
58 results->AddObject(new GenericExternalResult<double*>(results->Size()+1,"error",error,ninterp,1,1,0));
59 for(Object* & object :results->objects){
60 ExternalResult* result=xDynamicCast<ExternalResult*>(object);
61 result->WriteData(output_fid,1);
62 }
63 pfclose(output_fid,outbinfilename);
64 }
65
66 /*Close output and toolkits options file and write lock file if requested*/
67 _printf0_("write lock file:\n");
68 WriteLockFile(lockfilename);
69
70 /*Free ressources */
71 xDelete<char>(lockfilename);
72 xDelete<char>(binfilename);
73 xDelete<char>(outbinfilename);
74 xDelete<char>(rootpath);
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 delete options;
83 delete results;
84
85 /*Finalize environment:*/
86 EnvironmentFinalize();
87
88 /*Finalize exception trapping: */
89 ExceptionTrapEnd();
90
91 return 0; //unix success return;
92}
93
94void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv){
95 char *modelname = NULL;
96 char *binfilename = NULL;
97 char *outbinfilename = NULL;
98 char *lockfilename = NULL;
99 char *rootpatharg = NULL;
100 char *rootpath = NULL;
101
102 if(argc<1)_error_("Usage error: no execution path provided");
103 if(argc<2)_error_("Usage error: missing model name");
104
105 rootpatharg=argv[1];
106 if(strcmp(strstr(rootpatharg,"/"),"/")!=0){
107 rootpath = xNew<char>(strlen(rootpatharg)+2); sprintf(rootpath,"%s/",rootpatharg);
108 }
109 else{
110 rootpath = xNew<char>(strlen(rootpatharg)+1); sprintf(rootpath,"%s",rootpatharg);
111 }
112
113 modelname=argv[2];
114 if(strstr(modelname,rootpath)==NULL){
115 binfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".bin") +1); sprintf(binfilename, "%s%s%s",rootpath,modelname,".bin");
116 outbinfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s%s",rootpath,modelname,".outbin");
117 lockfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".lock") +1); sprintf(lockfilename, "%s%s%s",rootpath,modelname,".lock");
118 }
119 else{
120 binfilename = xNew<char>(strlen(modelname)+strlen(".bin") +1); sprintf(binfilename, "%s%s",modelname,".bin");
121 outbinfilename = xNew<char>(strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s",modelname,".outbin");
122 lockfilename = xNew<char>(strlen(modelname)+strlen(".lock") +1); sprintf(lockfilename, "%s%s",modelname,".lock");
123 }
124
125 /*Clean up and assign output pointer*/
126 *pbinfilename=binfilename;
127 *poutbinfilename=outbinfilename;
128 *plockfilename=lockfilename;
129 *prootpath=rootpath;
130}
131
132void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){
133
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 Options *options = NULL;
141
142 int M,N;
143 IoModel* iomodel = new IoModel();
144 iomodel->fid=fid;
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);
149 iomodel->FetchData(&x_interp,&M,&N,"md.x_interp"); ninterp=M*N;
150 iomodel->FetchData(&y_interp,&M,&N,"md.y_interp"); _assert_(M*N==ninterp);
151
152 /*Read options*/
153 options = new Options();
154 iomodel->FetchData(options,"md.y_interp");
155
156 /*Assign output pointer*/
157 *px = x;
158 *py = y;
159 *pdata = data;
160 *pnobs = nobs;
161 *px_interp = x_interp;
162 *py_interp = y_interp;
163 *pninterp = ninterp;
164 *poptions = options;
165 delete iomodel;
166}
Note: See TracBrowser for help on using the repository browser.