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

Last change on this file since 16137 was 16137, checked in by Mathieu Morlighem, 12 years ago

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