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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3317 - (hide annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 11 months ago) by caltinay
File size: 28309 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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