source: issm/branches/trunk-jpl-damage/src/c/solutions/kriging.cpp@ 13101

Last change on this file since 13101 was 13101, checked in by cborstad, 13 years ago

merged trunk-jpl through revision 13099 into branch
File size: 5.1 KB
Line 
1/*!\file: kriging.cpp
2 * \brief: kriging main parallel program
3 */
4
5#include "../issm.h"
6#include "../include/globals.h"
7
8/*Local prototypes*/
9void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,int argc,char **argv);
10void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid);
11
12int main(int argc,char **argv){
13
14 /*I/O: */
15 FILE *output_fid = NULL;
16 FILE *input_fid = NULL;
17 bool waitonlock = false;
18
19 /*File names*/
20 char *lockfilename = NULL;
21 char *binfilename = NULL;
22 char *outbinfilename = 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 ISSMBOOT();
38
39 /*Initialize environments: Petsc, MPI, etc...: */
40#ifdef _HAVE_PETSC_
41 int ierr=PetscInitialize(&argc,&argv,(char*)0,"");
42 if(ierr) _error_("Could not initialize Petsc");
43#else
44#ifdef _HAVE_MPI_
45 MPI_Init(&argc,&argv);
46#endif
47#endif
48
49 /*Size and rank: */
50#ifdef _HAVE_MPI_
51 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
52 MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
53#endif
54
55 /*First process inputs*/
56 _pprintLine_("");
57 _pprintLine_("Ice Sheet System Model (" << PACKAGE_NAME << ") version " << PACKAGE_VERSION);
58 _pprintLine_("(website: " << PACKAGE_URL << " contact: " << PACKAGE_BUGREPORT << ")");
59 _pprintLine_("");
60 ProcessArguments2(&binfilename,&outbinfilename,&lockfilename,argc,argv);
61
62 /*Process input files*/
63 input_fid=pfopen(binfilename,"rb");
64 ProcessInputfile(&x,&y,&data,&nobs,&x_interp,&y_interp,&ninterp,&options,input_fid);
65 pfclose(input_fid,binfilename);
66
67 _pprintLine_("call computational core:");
68 pKrigingx(&predictions,&error,x,y,data,nobs,x_interp,y_interp,ninterp,options);
69
70 _pprintLine_("write results to disk:");
71 Results *results = new Results();
72 if(my_rank==0){
73 output_fid=pfopen(outbinfilename,"wb");
74 results->AddObject(new DoubleVecExternalResult(results->Size()+1,0,predictions,ninterp,1,0));
75 results->AddObject(new DoubleVecExternalResult(results->Size()+1,1,error,ninterp,1,0));
76 for(int i=0;i<results->Size();i++){
77 ExternalResult* result=(ExternalResult*)results->GetObjectByOffset(i);
78 result->WriteData(output_fid,1);
79 }
80 pfclose(output_fid,outbinfilename);
81 }
82
83 /*Close output and petsc options file and write lock file if requested*/
84 _pprintLine_("write lock file:");
85 WriteLockFile(lockfilename);
86
87 /*Free ressources */
88 xDelete<char>(lockfilename);
89 xDelete<char>(binfilename);
90 xDelete<char>(outbinfilename);
91 xDelete<IssmDouble>(x);
92 xDelete<IssmDouble>(y);
93 xDelete<IssmDouble>(data);
94 xDelete<IssmDouble>(x_interp);
95 xDelete<IssmDouble>(y_interp);
96 xDelete<IssmDouble>(predictions);
97 xDelete<IssmDouble>(error);
98 delete options;
99 delete results;
100
101#ifdef _HAVE_PETSC_
102 _pprintLine_("closing MPI and Petsc");
103 PetscFinalize();
104#else
105#ifdef _HAVE_MPI_
106 _pprintLine_("closing MPI and Petsc");
107 MPI_Finalize();
108#endif
109#endif
110
111 /*end module: */
112 ISSMEND();
113
114 return 0; //unix success return;
115}
116
117void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,int argc,char **argv){
118
119 char *modelname = NULL;
120 char *binfilename = NULL;
121 char *outbinfilename = NULL;
122 char *lockfilename = NULL;
123
124 if(argc<2)_error_("Usage error: missing model name");
125 modelname=argv[2];
126 binfilename = xNew<char>((strlen(modelname)+strlen(".bin") +1)); sprintf(binfilename, "%s%s",modelname,".bin");
127 outbinfilename = xNew<char>((strlen(modelname)+strlen(".outbin")+1)); sprintf(outbinfilename,"%s%s",modelname,".outbin");
128 lockfilename = xNew<char>((strlen(modelname)+strlen(".lock") +1)); sprintf(lockfilename, "%s%s",modelname,".lock");
129
130 /*Clean up and assign output pointer*/
131 *pbinfilename=binfilename;
132 *poutbinfilename=outbinfilename;
133 *plockfilename=lockfilename;
134}
135
136void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid){
137
138 int ninterp,nobs,numoptions;
139 IssmDouble *x = NULL;
140 IssmDouble *y = NULL;
141 IssmDouble *data = NULL;
142 IssmDouble *x_interp = NULL;
143 IssmDouble *y_interp = NULL;
144 Options *options = NULL;
145 Option *option = NULL;
146
147 int dummy,M,N;
148 IoModel* iomodel = new IoModel();
149 iomodel->fid=fid;
150 iomodel->CheckEnumSync();
151 iomodel->FetchData(&x,&M,&N,0); nobs=M*N;
152 iomodel->FetchData(&y,&M,&N,1); _assert_(M*N==nobs);
153 iomodel->FetchData(&data,&M,&N,2); _assert_(M*N==nobs);
154 iomodel->FetchData(&x_interp,&M,&N,3); ninterp=M*N;
155 iomodel->FetchData(&y_interp,&M,&N,4); _assert_(M*N==ninterp);
156
157 /*Read options*/
158 options = new Options();
159 iomodel->LastIndex(&N);
160 numoptions=(int)((N-4)/2);
161 for(int i=0;i<numoptions;i++){
162 iomodel->FetchData(&option,5+2*i);
163 options->AddOption(option);
164 option=NULL;
165 }
166
167 /*Assign output pointer*/
168 *px = x;
169 *py = y;
170 *pdata = data;
171 *pnobs = nobs;
172 *px_interp = x_interp;
173 *py_interp = y_interp;
174 *pninterp = ninterp;
175 *poptions = options;
176 delete iomodel;
177}
Note: See TracBrowser for help on using the repository browser.