/[escript]/trunk/escript/src/Utils.cpp
ViewVC logotype

Contents of /trunk/escript/src/Utils.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2637 - (show annotations)
Fri Aug 28 04:12:24 2009 UTC (9 years, 10 months ago) by jfenwick
File size: 8891 byte(s)
saveDataCSV looks good but does not do MPI yet.
Also not unit tested.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include <string.h>
16
17 // added for saveCSV
18 #include <boost/python.hpp>
19 #include "Data.h"
20
21 #include "Utils.h"
22 #include "DataVector.h"
23
24 #ifdef _OPENMP
25 #include <omp.h>
26 #endif
27
28 #ifdef PASO_MPI
29 #include <mpi.h>
30 #endif
31
32 #ifdef _WIN32
33 #include <WinSock2.h>
34 #else
35 #include <unistd.h>
36 #endif
37
38 namespace escript {
39
40 int getSvnVersion()
41 {
42 #ifdef SVN_VERSION
43 return SVN_VERSION;
44 #else
45 return 0;
46 #endif
47 }
48
49 /* This is probably not very robust, but it works on Savanna today and is useful for performance analysis */
50 int get_core_id() {
51 int processor_num=-1;
52 #ifdef CORE_ID1
53 FILE *fp;
54 int i, count_spaces=0;
55 char fname[100];
56 char buf[1000];
57
58 sprintf(fname, "/proc/%d/stat", getpid());
59 fp = fopen(fname, "r");
60 if (fp == NULL) return(-1);
61 fgets(buf, 1000, fp);
62 fclose(fp);
63
64 for (i=strlen(buf)-1; i>=0; i--) {
65 if (buf[i] == ' ') count_spaces++;
66 if (count_spaces == 4) break;
67 }
68 processor_num = atoi(&buf[i+1]);
69 #endif
70 return(processor_num);
71 }
72
73
74 void printParallelThreadCnt()
75 {
76 int mpi_iam=0, mpi_num=1;
77 char hname[64];
78
79 #ifdef HAVE_GETHOSTNAME
80 gethostname(hname, 64);
81 hname[63] = '\0';
82 #else
83 strcpy(hname, "unknown host");
84 #endif
85
86 #ifdef PASO_MPI
87 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
88 MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
89 #endif
90
91 #pragma omp parallel
92 {
93 int omp_iam=0, omp_num=1;
94 #ifdef _OPENMP
95 omp_iam = omp_get_thread_num(); /* Call in a parallel region */
96 omp_num = omp_get_num_threads();
97 #endif
98 #pragma omp critical (printthrdcount)
99 printf("printParallelThreadCounts: MPI=%03d/%03d OpenMP=%03d/%03d running on %s core %d\n",
100 mpi_iam, mpi_num, omp_iam, omp_num, hname, get_core_id());
101 }
102 }
103
104 void setNumberOfThreads(const int num_threads)
105 {
106
107 #ifdef _OPENMP
108 omp_set_num_threads(num_threads);
109 #endif
110
111 }
112
113 int getNumberOfThreads()
114 {
115 #ifdef _OPENMP
116 return omp_get_max_threads();
117 #else
118 return 1;
119 #endif
120
121 }
122
123 ESCRIPT_DLL_API int getMPISizeWorld() {
124 int mpi_num = 1;
125 #ifdef PASO_MPI
126 MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
127 #endif
128 return mpi_num;
129 }
130
131 ESCRIPT_DLL_API int getMPIRankWorld() {
132 int mpi_iam = 0;
133 #ifdef PASO_MPI
134 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
135 #endif
136 return mpi_iam;
137 }
138
139 ESCRIPT_DLL_API int getMPIWorldMax(const int val) {
140 #ifdef PASO_MPI
141 int val2 = val;
142 int out = val;
143 MPI_Allreduce( &val2, &out, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD );
144 #else
145 int out = val;
146 #endif
147 return out;
148 }
149
150 ESCRIPT_DLL_API int getMPIWorldSum(const int val) {
151 #ifdef PASO_MPI
152 int val2 = val;
153 int out = 0;
154 MPI_Allreduce( &val2, &out, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
155 #else
156 int out = val;
157 #endif
158 return out;
159 }
160
161 ESCRIPT_DLL_API double getMachinePrecision() {
162 return DBL_EPSILON;
163 }
164 ESCRIPT_DLL_API double getMaxFloat() {
165 return DBL_MAX;
166 }
167 ESCRIPT_DLL_API void MPIBarrierWorld() {
168 #ifdef PASO_MPI
169 MPI_Barrier(MPI_COMM_WORLD );
170 #endif
171 }
172
173 ESCRIPT_DLL_API
174 void
175 saveDataCSV(const std::string& filename, boost::python::dict arg, const std::string& sep, const std::string& csep,
176 bool append)
177 {
178 using std::cout;
179 using std::endl;
180 boost::python::list keys=arg.keys();
181 int numdata = boost::python::extract<int>(arg.attr("__len__")());
182 bool hasmask=arg.has_key("mask");
183 Data mask;
184 if (hasmask)
185 {
186 mask=boost::python::extract<escript::Data>(arg["mask"]);
187 keys.remove("mask");
188 numdata--;
189 if (mask.getDataPointRank()!=0)
190 {
191 throw DataException("saveDataCSVcpp: masks must be scalar.");
192 }
193 }
194 if (numdata<1)
195 {
196 throw DataException("saveDataCSVcpp: no data to save specified.");
197 }
198 std::vector<int> step(numdata);
199 std::vector<std::string> names(numdata);
200 std::vector<Data> data(numdata);
201 std::vector<const DataAbstract::ValueType::value_type*> samples(numdata);
202 std::vector<int> offset(numdata);
203 std::vector<int> fstypes(numdata); // FunctionSpace types for each data
204
205 // We need to interpret the samples correctly even if they are different types
206 // for this reason, we should interate over samples
207 for (int i=0;i<numdata;++i)
208 {
209 names[i]=boost::python::extract<std::string>(keys[i]);
210 data[i]=boost::python::extract<escript::Data>(arg[keys[i]]);
211 step[i]=(data[i].actsExpanded()?DataTypes::noValues(data[i].getDataPointShape()):0);
212 fstypes[i]=data[i].getFunctionSpace().getTypeCode();
213 if (i>0)
214 {
215 if (data[i].getDomain()!=data[i-1].getDomain())
216 {
217 throw DataException("saveDataCSVcpp: all data must be on the same domain.");
218 }
219 }
220 }
221 int bestfnspace=0;
222 if (!data[0].getDomain()->commonFunctionSpace(fstypes, bestfnspace))
223 {
224 throw DataException("saveDataCSVcpp: FunctionSpaces of data are incompatible");
225 }
226 // now we interpolate all data to the same type
227 FunctionSpace best(data[0].getDomain(),bestfnspace);
228 for (int i=0;i<numdata;++i)
229 {
230 data[i]=data[i].interpolate(best);
231 }
232 int numsamples=data[0].getNumSamples(); // these must be the same for all data
233 int dpps=data[0].getNumDataPointsPerSample();
234
235
236 std::ofstream os;
237 if (append)
238 {
239 os.open(filename.c_str(), std::ios_base::app);
240 }
241 else
242 {
243 os.open(filename.c_str());
244 }
245 if (!os.is_open())
246 {
247 throw DataException("saveDataCSVcpp: unable to open file for writing");
248 }
249
250 bool first=true;
251 for (int i=0;i<numdata;++i)
252 {
253 const DataTypes::ShapeType& s=data[i].getDataPointShape();
254 switch (data[i].getDataPointRank())
255 {
256 case 0: if (!first)
257 {
258 os << sep;
259 }
260 else
261 {
262 first=false;
263 }
264 os << names[i]; break;
265 case 1: for (int j=0;j<s[0];++j)
266 {
267 if (!first)
268 {
269 os << sep;
270 }
271 else
272 {
273 first=false;
274 }
275 os << names[i] << csep << j;
276 }
277 break;
278 case 2: for (int j=0;j<s[0];++j)
279 {
280 for (int k=0;k<s[1];++k)
281 {
282 if (!first)
283 {
284 os << sep;
285 }
286 else
287 {
288 first=false;
289 }
290 os << names[i] << csep << k << csep << j;
291 }
292 }
293 break;
294 case 3: for (int j=0;j<s[0];++j)
295 {
296 for (int k=0;k<s[1];++k)
297 {
298 for (int l=0;l<s[2];++l)
299 {
300 if (!first)
301 {
302 os << sep;
303 }
304 else
305 {
306 first=false;
307 }
308 os << names[i] << csep << k << csep << j << csep << l;
309 }
310 }
311 }
312 break;
313 case 4: for (int j=0;j<s[0];++j)
314 {
315 for (int k=0;k<s[1];++k)
316 {
317 for (int l=0;l<s[2];++l)
318 {
319 for (int m=0;m<s[3];++m)
320 {
321 if (!first)
322 {
323 os << sep;
324 }
325 else
326 {
327 first=false;
328 }
329 os << names[i] << csep << k << csep << j << csep << l << csep << m;
330 }
331 }
332 }
333 }
334 break;
335 default:
336 throw DataException("saveDataCSV: Illegal rank");
337 }
338 }
339 os << endl;
340
341 boost::scoped_ptr<BufferGroup> maskbuffer; // sample buffer for the mask [if we have one]
342 const double* masksample=0;
343 int maskoffset=0;
344 //the use of shared_ptr here is just to ensure the buffer group is freed
345 //I would have used scoped_ptr but they don't work in vectors
346 std::vector<boost::shared_ptr<BufferGroup> > bg(numdata);
347 for (int d=0;d<numdata;++d)
348 {
349 bg[d].reset(data[d].allocSampleBuffer());
350 }
351
352 bool expandedmask=false; // does the mask act expanded. Are there mask value for each point in the sample
353 bool wantrow=true; // do we output this row?
354 if (hasmask)
355 {
356 maskbuffer.reset(mask.allocSampleBuffer());
357 if (mask.actsExpanded())
358 {
359 maskoffset=DataTypes::noValues(mask.getDataPointShape());
360 expandedmask=true;
361 }
362 }
363 try{
364 for (int i=0;i<numsamples;++i)
365 {
366 wantrow=true;
367 for (int d=0;d<numdata;++d)
368 {
369 samples[d]=data[d].getSampleDataRO(i,bg[d].get());
370 }
371 if (hasmask)
372 {
373 masksample=mask.getSampleDataRO(i, maskbuffer.get());
374 if (!expandedmask) // mask controls whole sample
375 {
376 if (masksample[0]<=0) // masks are scalar
377 {
378 wantrow=false;
379 }
380 }
381 }
382 for (int j=0;j<dpps;++j)
383 {
384 // now we need to check if this point is masked off
385 if (expandedmask)
386 {
387 wantrow=(masksample[j]>0); // masks are scalar to the relevant value is at [j]
388 }
389 if (wantrow)
390 {
391 bool needsep=false;
392 for (int d=0;d<numdata;++d)
393 {
394 DataTypes::pointToStream(os, samples[d], data[d].getDataPointShape(), offset[d], needsep, sep);
395 needsep=true;
396 offset[d]+=step[d];
397 }
398 os << endl;
399 }
400 }
401 for (int d=0;d<numdata;++d)
402 {
403 offset[d]=0;
404 }
405 }
406 } catch (...)
407 {
408 os.close();
409 throw;
410 }
411 os.close();
412
413 cout << "This method is not MPI safe" << endl;
414
415 }
416
417 } // end of namespace

  ViewVC Help
Powered by ViewVC 1.1.26