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

Annotation of /trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3675 - (hide annotations)
Thu Nov 17 00:53:38 2011 UTC (7 years, 10 months ago) by jfenwick
File size: 32778 byte(s)
pasowrap joins the trunk.

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 ksteube 1800 #include "Data.h"
16 jgs 480 #include "DataExpanded.h"
17 jgs 474 #include "DataException.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20 jfenwick 3506 #include <limits>
21 caltinay 3317
22     #include "esysUtils/Esys_MPI.h"
23    
24 gross 1023 #ifdef USE_NETCDF
25 ksteube 1312 #include <netcdfcpp.h>
26 gross 1023 #endif
27 jgs 82
28     #include <boost/python/extract.hpp>
29 jfenwick 1796 #include "DataMaths.h"
30 jgs 82
31 jfenwick 3506 //#define MKLRANDOM
32    
33     #ifdef MKLRANDOM
34     #include <mkl_vsl.h>
35 jfenwick 3550 #else
36     #include <boost/random/mersenne_twister.hpp>
37 jfenwick 3506 #endif
38    
39 jgs 82 using namespace std;
40     using namespace boost::python;
41     using namespace boost;
42 jfenwick 1796 using namespace escript::DataTypes;
43 jgs 82
44 jfenwick 2271
45     // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
46    
47     #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; *((int*)0)=17;throw DataException(ss.str());}
48    
49 jgs 82 namespace escript {
50    
51 jfenwick 2271 DataExpanded::DataExpanded(const WrappedArray& value,
52 jgs 102 const FunctionSpace& what)
53 jfenwick 2271 : parent(what,value.getShape())
54 jgs 102 {
55     //
56 jgs 117 // initialise the data array for this object
57 jfenwick 1796 initialise(what.getNumSamples(),what.getNumDPPSample());
58 jgs 102 //
59     // copy the given value to every data point
60     copy(value);
61     }
62 jgs 82
63 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
64 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape()),
65 jgs 102 m_data(other.m_data)
66 jgs 117 {
67 jgs 102 }
68 jgs 82
69 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
70 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape())
71 jgs 117 {
72     //
73     // initialise the data array for this object
74 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
75 jgs 102 //
76 jgs 117 // DataConstant only has one value, copy this to every data point
77 jfenwick 1796 copy(other);
78 jgs 102 }
79 jgs 82
80 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
81 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape())
82 jgs 117 {
83     //
84 jgs 119 // initialise the data array for this object
85 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
86 jgs 117 //
87     // for each data point in this object, extract and copy the corresponding data
88     // value from the given DataTagged object
89 jgs 102 int i,j;
90 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
91     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
92 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
93 jgs 117 for (i=0;i<numRows;i++) {
94     for (j=0;j<numCols;j++) {
95 jgs 102 try {
96 jfenwick 2271 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
97     other.getVectorRO(),
98 jgs 122 other.getPointOffset(i,j));
99 jgs 82 }
100 jgs 102 catch (std::exception& e) {
101     cout << e.what() << endl;
102     }
103 jgs 82 }
104     }
105 jgs 102 }
106 jgs 82
107 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
108 jfenwick 1796 const DataTypes::RegionType& region)
109 jfenwick 2005 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
110 jgs 102 {
111     //
112     // get the shape of the slice
113 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
114 jgs 102 //
115     // initialise this Data object to the shape of the slice
116 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
117 jgs 102 //
118 jgs 117 // copy the data
119 jfenwick 1796 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
120     DataTypes::ValueType::size_type numRows=m_data.getNumRows();
121     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
122 jgs 102 int i,j;
123 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
124 jgs 117 for (i=0;i<numRows;i++) {
125     for (j=0;j<numCols;j++) {
126 jgs 102 try {
127 jfenwick 2271 DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
128     other.getVectorRO(),
129 jfenwick 1796 other.getShape(),
130 jgs 122 other.getPointOffset(i,j),
131     region_loop_range);
132 jgs 82 }
133 jgs 102 catch (std::exception& e) {
134     cout << e.what() << endl;
135     }
136 jgs 82 }
137     }
138 jgs 102 }
139 jgs 82
140 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
141 jfenwick 1796 const DataTypes::ShapeType &shape,
142     const DataTypes::ValueType &data)
143 jfenwick 2005 : parent(what,shape)
144 jgs 119 {
145 jfenwick 1796 EsysAssert(data.size()%getNoValues()==0,
146     "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
147    
148     if (data.size()==getNoValues())
149     {
150     ValueType& vec=m_data.getData();
151     //
152     // create the view of the data
153     initialise(what.getNumSamples(),what.getNumDPPSample());
154     // now we copy this value to all elements
155     for (int i=0;i<getLength();)
156     {
157 phornby 2008 for (unsigned int j=0;j<getNoValues();++j,++i)
158 jfenwick 1796 {
159     vec[i]=data[j];
160     }
161     }
162     }
163     else
164     {
165     //
166     // copy the data in the correct format
167     m_data.getData()=data;
168     }
169    
170    
171 jgs 119 }
172    
173 jfenwick 3468 DataExpanded::DataExpanded(const FunctionSpace& what,
174     const DataTypes::ShapeType &shape,
175     const double v)
176     : parent(what,shape)
177     {
178     ValueType& vec=m_data.getData();
179     //
180     // create the view of the data
181     initialise(what.getNumSamples(),what.getNumDPPSample());
182     // now we copy this value to all elements
183     const int L=getLength();
184     int i;
185     #pragma omp parallel for schedule(static) private(i)
186     for (i=0;i<L;++i)
187     {
188     vec[i]=v;
189     }
190     }
191    
192    
193    
194 jgs 102 DataExpanded::~DataExpanded()
195     {
196     }
197 jgs 82
198 jgs 102 DataAbstract*
199 jfenwick 1799 DataExpanded::deepCopy()
200     {
201     return new DataExpanded(*this);
202     }
203    
204    
205     DataAbstract*
206 jfenwick 1796 DataExpanded::getSlice(const DataTypes::RegionType& region) const
207 jgs 102 {
208     return new DataExpanded(*this,region);
209     }
210    
211     void
212     DataExpanded::setSlice(const DataAbstract* value,
213 jfenwick 1796 const DataTypes::RegionType& region)
214 jgs 102 {
215     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
216     if (tempDataExp==0) {
217     throw DataException("Programming error - casting to DataExpanded.");
218     }
219 jfenwick 2271 CHECK_FOR_EX_WRITE
220 jgs 108 //
221 jgs 117 // get shape of slice
222 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
223     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
224 jgs 108 //
225 jgs 117 // check shape
226 jfenwick 1796 if (getRank()!=region.size()) {
227 jgs 102 throw DataException("Error - Invalid slice region.");
228     }
229 jfenwick 1796 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
230     throw DataException (DataTypes::createShapeErrorMessage(
231     "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
232 jgs 102 }
233     //
234 jgs 117 // copy the data from the slice into this object
235 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
236     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
237 jgs 102 int i, j;
238 jfenwick 2271 ValueType& vec=getVectorRW();
239 jfenwick 1796 const ShapeType& mshape=getShape();
240 jfenwick 2271 const ValueType& tVec=tempDataExp->getVectorRO();
241 jfenwick 1796 const ShapeType& tShape=tempDataExp->getShape();
242 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
243 jgs 102 for (i=0;i<numRows;i++) {
244     for (j=0;j<numCols;j++) {
245 jfenwick 1796 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
246     tVec,
247     tShape,
248     tempDataExp->getPointOffset(i,j),
249 jgs 122 region_loop_range);
250 jfenwick 1796
251 jgs 82 }
252     }
253 jgs 102 }
254 jgs 82
255 jgs 102 void
256 jfenwick 1796 DataExpanded::copy(const DataConstant& value)
257 jgs 102 {
258 jfenwick 1796 EsysAssert((checkShape(getShape(), value.getShape())),
259     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
260    
261 jgs 102 //
262 jgs 117 // copy a single value to every data point in this object
263 jgs 102 int nRows=m_data.getNumRows();
264     int nCols=m_data.getNumCols();
265 jgs 117 int i,j;
266 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
267 jgs 117 for (i=0;i<nRows;i++) {
268     for (j=0;j<nCols;j++) {
269 jfenwick 2271 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
270 jgs 82 }
271     }
272 jgs 102 }
273 jgs 82
274 jfenwick 2271 void
275     DataExpanded::copy(const WrappedArray& value)
276 jgs 102 {
277 jgs 117 // check the input shape matches this shape
278 jfenwick 2271 if (!DataTypes::checkShape(getShape(),value.getShape())) {
279 jfenwick 1796 throw DataException(DataTypes::createShapeErrorMessage(
280 jgs 102 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
281 jfenwick 2271 value.getShape(),getShape()));
282 jgs 82 }
283 jfenwick 2271 getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
284 jgs 102 }
285 jgs 82
286 jfenwick 1796
287 jgs 102 void
288 jfenwick 1796 DataExpanded::initialise(int noSamples,
289 jgs 102 int noDataPointsPerSample)
290     {
291 jfenwick 1808 if (noSamples==0) //retain the default empty object
292     {
293     return;
294     }
295 jgs 102 //
296 jgs 117 // resize data array to the required size
297 jfenwick 1796 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
298 jgs 102 }
299 jgs 82
300 jfenwick 2766 bool
301     DataExpanded::hasNaN() const
302     {
303     const ValueType& v=m_data.getData();
304     for (ValueType::size_type i=0;i<v.size();++i)
305     {
306     if (nancheck(v[i]))
307     {
308     return true;
309     }
310     }
311     return false;
312     }
313    
314    
315 jgs 102 string
316     DataExpanded::toString() const
317     {
318     stringstream temp;
319 gross 856 FunctionSpace fs=getFunctionSpace();
320 jfenwick 1796
321     int offset=0;
322 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
323     for (int j=0;j<m_data.getNumCols();j++) {
324 jfenwick 1796 offset=getPointOffset(i,j);
325 jgs 102 stringstream suffix;
326 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
327 jfenwick 2271 temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
328 jgs 102 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
329     temp << endl;
330 jgs 82 }
331     }
332     }
333 jfenwick 1808 string result=temp.str();
334     if (result.empty())
335     {
336     return "(data contains no samples)\n";
337     }
338 ksteube 1312 return temp.str();
339 jgs 102 }
340 jgs 82
341 jfenwick 1796 DataTypes::ValueType::size_type
342 jgs 102 DataExpanded::getPointOffset(int sampleNo,
343     int dataPointNo) const
344     {
345     return m_data.index(sampleNo,dataPointNo);
346     }
347 jgs 82
348 jfenwick 1796 DataTypes::ValueType::size_type
349 jfenwick 2005 DataExpanded::getPointOffset(int sampleNo,
350     int dataPointNo)
351     {
352     return m_data.index(sampleNo,dataPointNo);
353     }
354    
355     DataTypes::ValueType::size_type
356 jgs 102 DataExpanded::getLength() const
357     {
358 jgs 117 return m_data.size();
359 jgs 102 }
360 jgs 82
361 jgs 123
362    
363 matt 1319 void
364 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
365 jfenwick 2271 CHECK_FOR_EX_WRITE
366 gross 922 //
367     // Get the number of samples and data-points per sample.
368     int numSamples = getNumSamples();
369     int numDataPointsPerSample = getNumDPPSample();
370 jfenwick 1796 int dataPointRank = getRank();
371     ShapeType dataPointShape = getShape();
372 gross 922 if (numSamples*numDataPointsPerSample > 0) {
373     //TODO: global error handling
374 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
375 jfenwick 3675
376     ostringstream oss;
377     oss << "sampleNo=" << sampleNo << " numSamples=" << numSamples << endl;
378     cerr << oss.str();
379    
380    
381     throw DataException("Error - DataExpanded::copyToDataPoint[scalar] invalid sampleNo.");
382 gross 922 }
383 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
384 jfenwick 3675 throw DataException("Error - DataExpanded::copyToDataPoint[scalar] invalid dataPointNoInSample.");
385 gross 922 }
386 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
387 jfenwick 2271 ValueType& vec=getVectorRW();
388 gross 922 if (dataPointRank==0) {
389 gross 1846 vec[offset]=value;
390 gross 922 } else if (dataPointRank==1) {
391     for (int i=0; i<dataPointShape[0]; i++) {
392 jfenwick 1796 vec[offset+i]=value;
393 gross 922 }
394     } else if (dataPointRank==2) {
395     for (int i=0; i<dataPointShape[0]; i++) {
396     for (int j=0; j<dataPointShape[1]; j++) {
397 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
398 gross 922 }
399     }
400     } else if (dataPointRank==3) {
401     for (int i=0; i<dataPointShape[0]; i++) {
402     for (int j=0; j<dataPointShape[1]; j++) {
403     for (int k=0; k<dataPointShape[2]; k++) {
404 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
405 gross 922 }
406     }
407     }
408     } else if (dataPointRank==4) {
409     for (int i=0; i<dataPointShape[0]; i++) {
410     for (int j=0; j<dataPointShape[1]; j++) {
411     for (int k=0; k<dataPointShape[2]; k++) {
412     for (int l=0; l<dataPointShape[3]; l++) {
413 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
414 gross 922 }
415     }
416     }
417     }
418 matt 1319 }
419 gross 922 }
420     }
421 jfenwick 2271
422 matt 1319 void
423 jfenwick 2271 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
424     CHECK_FOR_EX_WRITE
425 gross 921 //
426     // Get the number of samples and data-points per sample.
427     int numSamples = getNumSamples();
428     int numDataPointsPerSample = getNumDPPSample();
429     //
430     // check rank:
431 jfenwick 2271 if (value.getRank()!=getRank())
432 jfenwick 2458 throw DataException("Rank of value does not match Data object rank");
433 gross 921 if (numSamples*numDataPointsPerSample > 0) {
434     //TODO: global error handling
435 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
436 jfenwick 3675
437     // ostringstream oss;
438     // oss << "sampleNo=" << sampleNo << " numSamples=" << numSamples << endl;
439     // throw DataException(oss.str());
440    
441     throw DataException("Error - DataExpanded::copyToDataPoint[wrapped] invalid sampleNo.");
442 gross 921 }
443 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
444 jfenwick 3675 throw DataException("Error - DataExpanded::copyToDataPoint[wrapped] invalid dataPointNoInSample.");
445 gross 921 }
446 gross 1962 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
447 jfenwick 2271 ValueType& vec=getVectorRW();
448     vec.copyFromArrayToOffset(value,offset,1);
449 gross 921 }
450     }
451 jfenwick 2271
452 jgs 126 void
453 ksteube 775 DataExpanded::symmetric(DataAbstract* ev)
454     {
455     int sampleNo,dataPointNo;
456     int numSamples = getNumSamples();
457     int numDataPointsPerSample = getNumDPPSample();
458     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
459     if (temp_ev==0) {
460     throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
461     }
462 jfenwick 2271 const ValueType& vec=getVectorRO();
463 jfenwick 1796 const ShapeType& shape=getShape();
464 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
465 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
466 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
467     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
468     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
469 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
470     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
471 ksteube 775 }
472     }
473     }
474     void
475     DataExpanded::nonsymmetric(DataAbstract* ev)
476     {
477     int sampleNo,dataPointNo;
478     int numSamples = getNumSamples();
479     int numDataPointsPerSample = getNumDPPSample();
480     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
481     if (temp_ev==0) {
482     throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
483     }
484 jfenwick 2271 const ValueType& vec=getVectorRO();
485 jfenwick 1796 const ShapeType& shape=getShape();
486 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
487 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
488 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
489     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
490     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
491 jfenwick 1796 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
492     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
493 ksteube 775 }
494     }
495     }
496     void
497 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
498 ksteube 775 {
499     int sampleNo,dataPointNo;
500     int numSamples = getNumSamples();
501     int numDataPointsPerSample = getNumDPPSample();
502     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
503     if (temp_ev==0) {
504 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
505 ksteube 775 }
506 jfenwick 2271 const ValueType& vec=getVectorRO();
507 jfenwick 1796 const ShapeType& shape=getShape();
508 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
509 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
510 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
511     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
512     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
513 jfenwick 1796 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
514     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
515 ksteube 775 }
516     }
517     }
518 gross 800
519 ksteube 775 void
520     DataExpanded::transpose(DataAbstract* ev, int axis_offset)
521     {
522     int sampleNo,dataPointNo;
523     int numSamples = getNumSamples();
524     int numDataPointsPerSample = getNumDPPSample();
525     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
526     if (temp_ev==0) {
527     throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
528     }
529 jfenwick 2271 const ValueType& vec=getVectorRO();
530 jfenwick 1796 const ShapeType& shape=getShape();
531 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
532 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
533 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
534     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
535     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
536 jfenwick 1796 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
537     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
538 ksteube 775 }
539     }
540     }
541 gross 800
542 ksteube 775 void
543 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
544 gross 800 {
545     int sampleNo,dataPointNo;
546     int numSamples = getNumSamples();
547     int numDataPointsPerSample = getNumDPPSample();
548     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
549     if (temp_ev==0) {
550 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
551 gross 800 }
552 jfenwick 2271 const ValueType& vec=getVectorRO();
553 jfenwick 1796 const ShapeType& shape=getShape();
554 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
555 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
556 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
557     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
558     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
559 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
560     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
561 gross 800 }
562     }
563     }
564     void
565 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
566     {
567 gross 584 int sampleNo,dataPointNo;
568 gross 580 int numSamples = getNumSamples();
569     int numDataPointsPerSample = getNumDPPSample();
570     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
571     if (temp_ev==0) {
572     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
573     }
574 jfenwick 2271 const ValueType& vec=getVectorRO();
575 jfenwick 1796 const ShapeType& shape=getShape();
576 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
577 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
578 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
579 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
580     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
581 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
582     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
583 gross 580 }
584     }
585     }
586     void
587     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
588     {
589     int numSamples = getNumSamples();
590     int numDataPointsPerSample = getNumDPPSample();
591 gross 584 int sampleNo,dataPointNo;
592 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
593     if (temp_ev==0) {
594     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
595     }
596     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
597     if (temp_V==0) {
598     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
599     }
600 jfenwick 2271 const ValueType& vec=getVectorRO();
601 jfenwick 1796 const ShapeType& shape=getShape();
602 jfenwick 2271 ValueType& evVec=temp_ev->getVectorRW();
603 jfenwick 1796 const ShapeType& evShape=temp_ev->getShape();
604 jfenwick 2271 ValueType& VVec=temp_V->getVectorRW();
605 jfenwick 1796 const ShapeType& VShape=temp_V->getShape();
606 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
607 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
608     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
609 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
610     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
611     VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
612 gross 580 }
613     }
614     }
615    
616 jfenwick 2742
617 jfenwick 2792 int
618 jfenwick 2742 DataExpanded::matrixInverse(DataAbstract* out) const
619     {
620     DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
621     if (temp==0)
622     {
623     throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
624     }
625    
626     if (getRank()!=2)
627     {
628     throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
629    
630     }
631     int sampleNo;
632     const int numdpps=getNumDPPSample();
633     const int numSamples = getNumSamples();
634     const ValueType& vec=m_data.getData();
635     int errcode=0;
636     #pragma omp parallel private(sampleNo)
637     {
638     int errorcode=0;
639     LapackInverseHelper h(getShape()[0]);
640     #pragma omp for schedule(static)
641     for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
642     {
643     // not sure I like all those virtual calls to getPointOffset
644     DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
645     int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
646     if (res>errorcode)
647     {
648     errorcode=res;
649     #pragma omp critical
650     {
651     errcode=errorcode; // I'm not especially concerned which error gets reported as long as one is
652     }
653     }
654     }
655     }
656 jfenwick 2792 return errcode;
657 jfenwick 2742 if (errcode)
658     {
659     DataMaths::matrixInverseError(errcode); // throws exceptions
660     }
661     }
662    
663     void
664 gross 1118 DataExpanded::setToZero(){
665 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
666     // Why is there no memset here? Parallel issues?
667 jfenwick 2742 // A: This ensures that memory is touched by the correct thread.
668 jfenwick 2271 CHECK_FOR_EX_WRITE
669 gross 1118 int numSamples = getNumSamples();
670     int numDataPointsPerSample = getNumDPPSample();
671 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
672 gross 1118 double* p;
673     int sampleNo,dataPointNo, i;
674     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
675     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
676     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
677     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
678 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
679 gross 1118 }
680     }
681     }
682    
683 ksteube 1748 /* Append MPI rank to file name if multiple MPI processes */
684     char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
685     /* Make plenty of room for the mpi_rank number and terminating '\0' */
686     char *newFileName = (char *)malloc(strlen(fileName)+20);
687     strncpy(newFileName, fileName, strlen(fileName)+1);
688     if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
689     return(newFileName);
690     }
691 gross 1118
692     void
693 gross 950 DataExpanded::dump(const std::string fileName) const
694     {
695 gross 1023 #ifdef USE_NETCDF
696 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
697 gross 967 const NcDim* ncdims[ldims];
698     NcVar *var, *ids;
699 jfenwick 1796 int rank = getRank();
700 gross 967 int type= getFunctionSpace().getTypeCode();
701     int ndims =0;
702     long dims[ldims];
703 gross 1141 const double* d_ptr=&(m_data[0]);
704 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
705 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
706     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
707 jfenwick 3259 #ifdef ESYS_MPI
708 ksteube 1800 MPI_Status status;
709     #endif
710 gross 965
711 jfenwick 3259 #ifdef ESYS_MPI
712 ksteube 1827 /* Serialize NetCDF I/O */
713     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
714     #endif
715    
716 gross 967 // netCDF error handler
717     NcError err(NcError::verbose_nonfatal);
718     // Create the file.
719 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
720     NcFile dataFile(newFileName, NcFile::Replace);
721 gross 967 // check if writing was successful
722     if (!dataFile.is_valid())
723     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
724 gross 1141 if (!dataFile.add_att("type_id",2) )
725 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
726     if (!dataFile.add_att("rank",rank) )
727     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
728     if (!dataFile.add_att("function_space_type",type))
729     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
730     ndims=rank+2;
731     if ( rank >0 ) {
732     dims[0]=shape[0];
733     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
734 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
735 gross 967 }
736     if ( rank >1 ) {
737     dims[1]=shape[1];
738     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
739 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
740 gross 967 }
741     if ( rank >2 ) {
742     dims[2]=shape[2];
743     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
744 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
745 gross 967 }
746     if ( rank >3 ) {
747     dims[3]=shape[3];
748     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
749 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
750 gross 967 }
751     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
752     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
753     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
754     dims[rank+1]=getFunctionSpace().getNumSamples();
755     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
756     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
757    
758 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
759     {
760    
761     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
762 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
763 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
764     if (! (ids->put(ids_p,dims[rank+1])) )
765 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
766 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
767 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
768 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
769 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
770 jfenwick 1808 }
771 jfenwick 3259 #ifdef ESYS_MPI
772 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
773 ksteube 1800 #endif
774 gross 1023 #else
775     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
776     #endif
777 gross 950 }
778 gross 1358
779 jfenwick 1796 void
780 gross 1358 DataExpanded::setTaggedValue(int tagKey,
781 jfenwick 1796 const DataTypes::ShapeType& pointshape,
782     const DataTypes::ValueType& value,
783     int dataOffset)
784 gross 1358 {
785 jfenwick 2271 CHECK_FOR_EX_WRITE
786 gross 1358 int numSamples = getNumSamples();
787     int numDataPointsPerSample = getNumDPPSample();
788     int sampleNo,dataPointNo, i;
789 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
790     double* p;
791     const double* in=&value[0+dataOffset];
792 gross 1358
793 jfenwick 1796 if (value.size() != n) {
794 gross 1358 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
795     }
796    
797     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
798     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
799     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
800     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
801     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
802 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
803 gross 1358 }
804     }
805     }
806     }
807    
808 jfenwick 1796
809 gross 1487 void
810     DataExpanded::reorderByReferenceIDs(int *reference_ids)
811     {
812 jfenwick 2271 CHECK_FOR_EX_WRITE
813 gross 1487 int numSamples = getNumSamples();
814 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
815 gross 1487 int sampleNo, sampleNo2,i;
816     double* p,*p2;
817     register double rtmp;
818     FunctionSpace fs=getFunctionSpace();
819 gross 1358
820 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
821     const int id_in=reference_ids[sampleNo];
822     const int id=fs.getReferenceIDOfSample(sampleNo);
823     if (id!=id_in) {
824     bool matched=false;
825     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
826     if (id == reference_ids[sampleNo2]) {
827     p=&(m_data[getPointOffset(sampleNo,0)]);
828     p2=&(m_data[getPointOffset(sampleNo2,0)]);
829 gross 1513 for (i=0; i<n ;i++) {
830 gross 1487 rtmp=p[i];
831     p[i]=p2[i];
832     p2[i]=rtmp;
833     }
834     reference_ids[sampleNo]=id;
835     reference_ids[sampleNo2]=id_in;
836     matched=true;
837     break;
838     }
839     }
840 phornby 1628 if (! matched) {
841 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
842     }
843     }
844     }
845     }
846    
847 jfenwick 1796 DataTypes::ValueType&
848 jfenwick 2271 DataExpanded::getVectorRW()
849 jfenwick 1796 {
850 jfenwick 2271 CHECK_FOR_EX_WRITE
851 jfenwick 1796 return m_data.getData();
852     }
853 gross 1487
854 jfenwick 1796 const DataTypes::ValueType&
855 jfenwick 2271 DataExpanded::getVectorRO() const
856 jfenwick 1796 {
857     return m_data.getData();
858     }
859    
860 jfenwick 3506
861 jfenwick 3550 #ifndef MKLRANDOM
862     namespace {
863    
864     boost::mt19937 base; // used to seed all the other generators
865     vector<boost::mt19937> gens;
866    
867    
868     void seedGens(long seed)
869     {
870     #ifdef _OPENMP
871     int numthreads=omp_get_max_threads();
872     #else
873     int numthreads=1;
874     #endif
875     if (gens.size()==0) // we haven't instantiated the generators yet
876     {
877     gens.resize(numthreads); // yes this means all the generators will be owned by one thread in a NUMA sense
878     } // I don't think it is worth a more complicated solution at this point
879     if (seed!=0)
880     {
881 jfenwick 3552 base.seed((uint32_t)seed); // without this cast, icc gets confused
882 jfenwick 3550 for (int i=0;i<numthreads;++i)
883     {
884 jfenwick 3552 uint32_t b=base();
885     gens[i].seed(b); // initialise each generator with successive random values
886 jfenwick 3550 }
887     }
888     }
889    
890    
891     }
892     #endif
893    
894 jfenwick 3506 // Idea here is to create an array of seeds by feeding the original seed into the random generator
895     // The code at the beginning of the function to compute the seed if one is given is
896     // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
897     // I make no claim about how well these initial seeds are distributed
898     void DataExpanded::randomFill(long seed)
899 jfenwick 3390 {
900     CHECK_FOR_EX_WRITE
901 jfenwick 3506 static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
902     if (seed==0) // for each one
903 jfenwick 3390 {
904 jfenwick 3506 if (prevseed==0)
905     {
906     time_t s=time(0);
907     seed=s;
908     }
909     else
910     {
911     seed=prevseed+419; // these numbers are arbitrary
912     if (seed>3040101) // I want to avoid overflow on 32bit systems
913     {
914     seed=((int)(seed)%0xABCD)+1;
915     }
916     }
917 jfenwick 3390 }
918 jfenwick 3506 // now we need to consider MPI since we don't want each rank to start with the same seed
919     seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
920     prevseed=seed;
921 jfenwick 3550
922     #ifdef MKLRANDOM
923    
924 jfenwick 3506 #ifdef _OPENMP
925     int numthreads=omp_get_max_threads();
926     #else
927     int numthreads=1;
928     #endif
929     double* seeds=new double[numthreads];
930     VSLStreamStatePtr sstream;
931    
932     int status=vslNewStream(&sstream, VSL_BRNG_MT19937, seed); // use a Mersenne Twister
933     numeric_limits<double> dlim;
934     vdRngUniform(VSL_METHOD_DUNIFORM_STD, sstream , numthreads, seeds, -1, 1);
935     vslDeleteStream(&sstream);
936     DataVector& dv=getVectorRW();
937     size_t dvsize=dv.size();
938     #pragma omp parallel
939     {
940     int tnum=0;
941     #ifdef _OPENMP
942     tnum=omp_get_thread_num();
943     #endif
944     VSLStreamStatePtr stream;
945     // the 12345 is a hack to give us a better chance of getting different integer seeds.
946     int status=vslNewStream(&stream, VSL_BRNG_MT19937, seeds[tnum]*12345); // use a Mersenne Twister
947     int bigchunk=(dvsize/numthreads+1);
948     int smallchunk=dvsize-bigchunk*(numthreads-1);
949     int chunksize=(tnum<(numthreads-1))?bigchunk:smallchunk;
950     vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream, chunksize, &(dv[bigchunk*tnum]), 0,1);
951     vslDeleteStream(&stream);
952     }
953     delete[] seeds;
954     #else
955 jfenwick 3550 boost::mt19937::result_type RMAX=base.max();
956     seedGens(seed);
957 jfenwick 3390 DataVector& dv=getVectorRW();
958 jfenwick 3506 long i;
959     const size_t dvsize=dv.size();
960 jfenwick 3550
961 jfenwick 3506 #pragma omp parallel private(i)
962 jfenwick 3390 {
963 jfenwick 3506 int tnum=0;
964     #ifdef _OPENMP
965     tnum=omp_get_thread_num();
966     #endif
967 jfenwick 3550 boost::mt19937& generator=gens[tnum];
968 jfenwick 3506
969     #pragma omp for schedule(static)
970     for (i=0;i<dvsize;++i)
971     {
972 jfenwick 3550
973    
974    
975    
976    
977 jfenwick 3544 #ifdef _WIN32
978 jfenwick 3550 dv[i]=((double)generator())/RMAX;
979 jfenwick 3544 #else
980 jfenwick 3550 dv[i]=((double)generator())/RMAX;
981 jfenwick 3544 #endif
982 jfenwick 3506 }
983 jfenwick 3390 }
984 jfenwick 3506 #endif
985 jfenwick 3390 }
986 jfenwick 2271
987 jgs 82 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26