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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

temp_trunk_copy/escript/src/Utils.cpp revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC trunk/escript/src/Utils.cpp revision 2637 by jfenwick, Fri Aug 28 04:12:24 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * 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"  #include "Utils.h"
22  #include "DataVector.h"  #include "DataVector.h"
# Line 20  Line 25 
25  #include <omp.h>  #include <omp.h>
26  #endif  #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 {  namespace escript {
39    
40  int getSvnVersion()  int getSvnVersion()
# Line 31  int getSvnVersion() Line 46  int getSvnVersion()
46  #endif  #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)  void setNumberOfThreads(const int num_threads)
105  {  {
106    
# Line 50  int getNumberOfThreads() Line 120  int getNumberOfThreads()
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  }  // end of namespace

Legend:
Removed from v.1384  
changed lines
  Added in v.2637

  ViewVC Help
Powered by ViewVC 1.1.26