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