1 |
// $Id$ |
2 |
/* |
3 |
************************************************************ |
4 |
* Copyright 2006 by ACcESS MNRF * |
5 |
* * |
6 |
* http://www.access.edu.au * |
7 |
* Primary Business: Queensland, Australia * |
8 |
* Licensed under the Open Software License version 3.0 * |
9 |
* http://www.opensource.org/licenses/osl-3.0.php * |
10 |
* * |
11 |
************************************************************ |
12 |
*/ |
13 |
|
14 |
#include "DataExpanded.h" |
15 |
#include "DataException.h" |
16 |
#include "DataConstant.h" |
17 |
#include "DataTagged.h" |
18 |
|
19 |
#include <boost/python/extract.hpp> |
20 |
|
21 |
using namespace std; |
22 |
using namespace boost::python; |
23 |
using namespace boost; |
24 |
|
25 |
namespace escript { |
26 |
|
27 |
DataExpanded::DataExpanded(const boost::python::numeric::array& value, |
28 |
const FunctionSpace& what) |
29 |
: DataAbstract(what) |
30 |
{ |
31 |
DataArrayView::ShapeType tempShape; |
32 |
// |
33 |
// extract the shape of the python numarray |
34 |
for (int i=0; i<value.getrank(); i++) { |
35 |
tempShape.push_back(extract<int>(value.getshape()[i])); |
36 |
} |
37 |
// |
38 |
// initialise the data array for this object |
39 |
initialise(tempShape,what.getNumSamples(),what.getNumDPPSample()); |
40 |
// |
41 |
// copy the given value to every data point |
42 |
copy(value); |
43 |
} |
44 |
|
45 |
DataExpanded::DataExpanded(const DataExpanded& other) |
46 |
: DataAbstract(other.getFunctionSpace()), |
47 |
m_data(other.m_data) |
48 |
{ |
49 |
// |
50 |
// create the view for the data |
51 |
DataArrayView temp(m_data.getData(),other.getPointDataView().getShape()); |
52 |
setPointDataView(temp); |
53 |
} |
54 |
|
55 |
DataExpanded::DataExpanded(const DataConstant& other) |
56 |
: DataAbstract(other.getFunctionSpace()) |
57 |
{ |
58 |
// |
59 |
// initialise the data array for this object |
60 |
initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample()); |
61 |
// |
62 |
// DataConstant only has one value, copy this to every data point |
63 |
copy(other.getPointDataView()); |
64 |
} |
65 |
|
66 |
DataExpanded::DataExpanded(const DataTagged& other) |
67 |
: DataAbstract(other.getFunctionSpace()) |
68 |
{ |
69 |
// |
70 |
// initialise the data array for this object |
71 |
initialise(other.getPointDataView().getShape(),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 |
DataArrayView::ValueType::size_type numRows=m_data.getNumRows(); |
77 |
DataArrayView::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 |
getPointDataView().copy(getPointOffset(i,j), |
83 |
other.getPointDataView(), |
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 DataArrayView::RegionType& region) |
95 |
: DataAbstract(other.getFunctionSpace()) |
96 |
{ |
97 |
// |
98 |
// get the shape of the slice |
99 |
DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region)); |
100 |
// |
101 |
// initialise this Data object to the shape of the slice |
102 |
initialise(shape,other.getNumSamples(),other.getNumDPPSample()); |
103 |
// |
104 |
// copy the data |
105 |
DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region); |
106 |
DataArrayView::ValueType::size_type numRows=m_data.getNumRows(); |
107 |
DataArrayView::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 |
} |
118 |
catch (std::exception& e) { |
119 |
cout << e.what() << endl; |
120 |
} |
121 |
} |
122 |
} |
123 |
} |
124 |
|
125 |
DataExpanded::DataExpanded(const DataArrayView& value, |
126 |
const FunctionSpace& what) |
127 |
: DataAbstract(what) |
128 |
{ |
129 |
// |
130 |
// get the shape of the given data value |
131 |
DataArrayView::ShapeType tempShape=value.getShape(); |
132 |
// |
133 |
// initialise this Data object to the shape of the given data value |
134 |
initialise(tempShape,what.getNumSamples(),what.getNumDPPSample()); |
135 |
// |
136 |
// copy the given value to every data point |
137 |
copy(value); |
138 |
} |
139 |
|
140 |
DataExpanded::DataExpanded(const FunctionSpace& what, |
141 |
const DataArrayView::ShapeType &shape, |
142 |
const DataArrayView::ValueType &data) |
143 |
: DataAbstract(what) |
144 |
{ |
145 |
// |
146 |
// create the view of the data |
147 |
initialise(shape,what.getNumSamples(),what.getNumDPPSample()); |
148 |
// |
149 |
// copy the data in the correct format |
150 |
m_data.getData()=data; |
151 |
} |
152 |
|
153 |
DataExpanded::~DataExpanded() |
154 |
{ |
155 |
} |
156 |
|
157 |
DataAbstract* |
158 |
DataExpanded::getSlice(const DataArrayView::RegionType& region) const |
159 |
{ |
160 |
return new DataExpanded(*this,region); |
161 |
} |
162 |
|
163 |
void |
164 |
DataExpanded::setSlice(const DataAbstract* value, |
165 |
const DataArrayView::RegionType& region) |
166 |
{ |
167 |
const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value); |
168 |
if (tempDataExp==0) { |
169 |
throw DataException("Programming error - casting to DataExpanded."); |
170 |
} |
171 |
// |
172 |
// get shape of slice |
173 |
DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region)); |
174 |
DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region); |
175 |
// |
176 |
// check shape |
177 |
if (getPointDataView().getRank()!=region.size()) { |
178 |
throw DataException("Error - Invalid slice region."); |
179 |
} |
180 |
if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) { |
181 |
throw DataException (value->getPointDataView().createShapeErrorMessage( |
182 |
"Error - Couldn't copy slice due to shape mismatch.",shape)); |
183 |
} |
184 |
// |
185 |
// copy the data from the slice into this object |
186 |
DataArrayView::ValueType::size_type numRows=m_data.getNumRows(); |
187 |
DataArrayView::ValueType::size_type numCols=m_data.getNumCols(); |
188 |
int i, j; |
189 |
#pragma omp parallel for private(i,j) schedule(static) |
190 |
for (i=0;i<numRows;i++) { |
191 |
for (j=0;j<numCols;j++) { |
192 |
getPointDataView().copySliceFrom(getPointOffset(i,j), |
193 |
tempDataExp->getPointDataView(), |
194 |
tempDataExp->getPointOffset(i,j), |
195 |
region_loop_range); |
196 |
} |
197 |
} |
198 |
} |
199 |
|
200 |
void |
201 |
DataExpanded::copy(const DataArrayView& value) |
202 |
{ |
203 |
// |
204 |
// copy a single value to every data point in this object |
205 |
int nRows=m_data.getNumRows(); |
206 |
int nCols=m_data.getNumCols(); |
207 |
int i,j; |
208 |
#pragma omp parallel for private(i,j) schedule(static) |
209 |
for (i=0;i<nRows;i++) { |
210 |
for (j=0;j<nCols;j++) { |
211 |
// NOTE: An exception may be thown from this call if |
212 |
// DOASSERT is on which of course will play |
213 |
// havoc with the omp threads. Run single threaded |
214 |
// if using DOASSERT. |
215 |
getPointDataView().copy(m_data.index(i,j),value); |
216 |
} |
217 |
} |
218 |
} |
219 |
|
220 |
void |
221 |
DataExpanded::copy(const boost::python::numeric::array& value) |
222 |
{ |
223 |
// |
224 |
// first convert the numarray into a DataArray object |
225 |
DataArray temp(value); |
226 |
// |
227 |
// check the input shape matches this shape |
228 |
if (!getPointDataView().checkShape(temp.getView().getShape())) { |
229 |
throw DataException(getPointDataView().createShapeErrorMessage( |
230 |
"Error - (DataExpanded) Cannot copy due to shape mismatch.", |
231 |
temp.getView().getShape())); |
232 |
} |
233 |
// |
234 |
// now copy over the data |
235 |
copy(temp.getView()); |
236 |
} |
237 |
|
238 |
void |
239 |
DataExpanded::initialise(const DataArrayView::ShapeType& shape, |
240 |
int noSamples, |
241 |
int noDataPointsPerSample) |
242 |
{ |
243 |
// |
244 |
// resize data array to the required size |
245 |
m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape)); |
246 |
// |
247 |
// create the data view of the data array |
248 |
DataArrayView temp(m_data.getData(),shape); |
249 |
setPointDataView(temp); |
250 |
} |
251 |
|
252 |
string |
253 |
DataExpanded::toString() const |
254 |
{ |
255 |
stringstream temp; |
256 |
FunctionSpace fs=getFunctionSpace(); |
257 |
// |
258 |
// create a temporary view as the offset will be changed |
259 |
DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset()); |
260 |
for (int i=0;i<m_data.getNumRows();i++) { |
261 |
for (int j=0;j<m_data.getNumCols();j++) { |
262 |
tempView.setOffset(m_data.index(i,j)); |
263 |
stringstream suffix; |
264 |
suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")"; |
265 |
temp << tempView.toString(suffix.str()); |
266 |
if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) { |
267 |
temp << endl; |
268 |
} |
269 |
} |
270 |
} |
271 |
return temp.str(); |
272 |
} |
273 |
|
274 |
DataArrayView::ValueType::size_type |
275 |
DataExpanded::getPointOffset(int sampleNo, |
276 |
int dataPointNo) const |
277 |
{ |
278 |
return m_data.index(sampleNo,dataPointNo); |
279 |
} |
280 |
|
281 |
DataArrayView |
282 |
DataExpanded::getDataPoint(int sampleNo, |
283 |
int dataPointNo) |
284 |
{ |
285 |
DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo)); |
286 |
return temp; |
287 |
} |
288 |
|
289 |
DataArrayView::ValueType::size_type |
290 |
DataExpanded::getLength() const |
291 |
{ |
292 |
return m_data.size(); |
293 |
} |
294 |
|
295 |
int |
296 |
DataExpanded::archiveData(ofstream& archiveFile, |
297 |
const DataArrayView::ValueType::size_type noValues) const |
298 |
{ |
299 |
return(m_data.archiveData(archiveFile, noValues)); |
300 |
} |
301 |
|
302 |
int |
303 |
DataExpanded::extractData(ifstream& archiveFile, |
304 |
const DataArrayView::ValueType::size_type noValues) |
305 |
{ |
306 |
return(m_data.extractData(archiveFile, noValues)); |
307 |
} |
308 |
|
309 |
void |
310 |
DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) { |
311 |
// |
312 |
// Get the number of samples and data-points per sample. |
313 |
int numSamples = getNumSamples(); |
314 |
int numDataPointsPerSample = getNumDPPSample(); |
315 |
int dataPointRank = getPointDataView().getRank(); |
316 |
ShapeType dataPointShape = getPointDataView().getShape(); |
317 |
if (numSamples*numDataPointsPerSample > 0) { |
318 |
//TODO: global error handling |
319 |
if ((sampleNo >= numSamples) || (sampleNo < 0 )) { |
320 |
throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo."); |
321 |
} |
322 |
if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) { |
323 |
throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample."); |
324 |
} |
325 |
DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo); |
326 |
if (dataPointRank==0) { |
327 |
dataPointView()=value; |
328 |
} else if (dataPointRank==1) { |
329 |
for (int i=0; i<dataPointShape[0]; i++) { |
330 |
dataPointView(i)=value; |
331 |
} |
332 |
} else if (dataPointRank==2) { |
333 |
for (int i=0; i<dataPointShape[0]; i++) { |
334 |
for (int j=0; j<dataPointShape[1]; j++) { |
335 |
dataPointView(i,j)=value; |
336 |
} |
337 |
} |
338 |
} else if (dataPointRank==3) { |
339 |
for (int i=0; i<dataPointShape[0]; i++) { |
340 |
for (int j=0; j<dataPointShape[1]; j++) { |
341 |
for (int k=0; k<dataPointShape[2]; k++) { |
342 |
dataPointView(i,j,k)=value; |
343 |
} |
344 |
} |
345 |
} |
346 |
} else if (dataPointRank==4) { |
347 |
for (int i=0; i<dataPointShape[0]; i++) { |
348 |
for (int j=0; j<dataPointShape[1]; j++) { |
349 |
for (int k=0; k<dataPointShape[2]; k++) { |
350 |
for (int l=0; l<dataPointShape[3]; l++) { |
351 |
dataPointView(i,j,k,l)=value; |
352 |
} |
353 |
} |
354 |
} |
355 |
} |
356 |
} |
357 |
} |
358 |
} |
359 |
void |
360 |
DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) { |
361 |
// |
362 |
// Get the number of samples and data-points per sample. |
363 |
int numSamples = getNumSamples(); |
364 |
int numDataPointsPerSample = getNumDPPSample(); |
365 |
int dataPointRank = getPointDataView().getRank(); |
366 |
ShapeType dataPointShape = getPointDataView().getShape(); |
367 |
// |
368 |
// check rank: |
369 |
if (value.getrank()!=dataPointRank) |
370 |
throw DataException("Rank of numarray does not match Data object rank"); |
371 |
if (numSamples*numDataPointsPerSample > 0) { |
372 |
//TODO: global error handling |
373 |
if ((sampleNo >= numSamples) || (sampleNo < 0 )) { |
374 |
throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo."); |
375 |
} |
376 |
if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) { |
377 |
throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample."); |
378 |
} |
379 |
DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo); |
380 |
if (dataPointRank==0) { |
381 |
dataPointView()=extract<double>(value[0]); |
382 |
} else if (dataPointRank==1) { |
383 |
for (int i=0; i<dataPointShape[0]; i++) { |
384 |
dataPointView(i)=extract<double>(value[i]); |
385 |
} |
386 |
} else if (dataPointRank==2) { |
387 |
for (int i=0; i<dataPointShape[0]; i++) { |
388 |
for (int j=0; j<dataPointShape[1]; j++) { |
389 |
dataPointView(i,j)=extract<double>(value[i][j]); |
390 |
} |
391 |
} |
392 |
} else if (dataPointRank==3) { |
393 |
for (int i=0; i<dataPointShape[0]; i++) { |
394 |
for (int j=0; j<dataPointShape[1]; j++) { |
395 |
for (int k=0; k<dataPointShape[2]; k++) { |
396 |
dataPointView(i,j,k)=extract<double>(value[i][j][k]); |
397 |
} |
398 |
} |
399 |
} |
400 |
} else if (dataPointRank==4) { |
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 |
for (int l=0; l<dataPointShape[3]; l++) { |
405 |
dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]); |
406 |
} |
407 |
} |
408 |
} |
409 |
} |
410 |
} |
411 |
} |
412 |
} |
413 |
void |
414 |
DataExpanded::copyAll(const boost::python::numeric::array& value) { |
415 |
// |
416 |
// Get the number of samples and data-points per sample. |
417 |
int numSamples = getNumSamples(); |
418 |
int numDataPointsPerSample = getNumDPPSample(); |
419 |
int dataPointRank = getPointDataView().getRank(); |
420 |
ShapeType dataPointShape = getPointDataView().getShape(); |
421 |
// |
422 |
// check rank: |
423 |
if (value.getrank()!=dataPointRank+1) |
424 |
throw DataException("Rank of numarray does not match Data object rank"); |
425 |
if (value.getshape()[0]!=numSamples*numDataPointsPerSample) |
426 |
throw DataException("leading dimension of numarray is too small"); |
427 |
// |
428 |
int dataPoint = 0; |
429 |
for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
430 |
for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
431 |
DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo); |
432 |
if (dataPointRank==0) { |
433 |
dataPointView()=extract<double>(value[dataPoint]); |
434 |
} else if (dataPointRank==1) { |
435 |
for (int i=0; i<dataPointShape[0]; i++) { |
436 |
dataPointView(i)=extract<double>(value[dataPoint][i]); |
437 |
} |
438 |
} else if (dataPointRank==2) { |
439 |
for (int i=0; i<dataPointShape[0]; i++) { |
440 |
for (int j=0; j<dataPointShape[1]; j++) { |
441 |
dataPointView(i,j)=extract<double>(value[dataPoint][i][j]); |
442 |
} |
443 |
} |
444 |
} else if (dataPointRank==3) { |
445 |
for (int i=0; i<dataPointShape[0]; i++) { |
446 |
for (int j=0; j<dataPointShape[1]; j++) { |
447 |
for (int k=0; k<dataPointShape[2]; k++) { |
448 |
dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]); |
449 |
} |
450 |
} |
451 |
} |
452 |
} else if (dataPointRank==4) { |
453 |
for (int i=0; i<dataPointShape[0]; i++) { |
454 |
for (int j=0; j<dataPointShape[1]; j++) { |
455 |
for (int k=0; k<dataPointShape[2]; k++) { |
456 |
for (int l=0; l<dataPointShape[3]; l++) { |
457 |
dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]); |
458 |
} |
459 |
} |
460 |
} |
461 |
} |
462 |
} |
463 |
dataPoint++; |
464 |
} |
465 |
} |
466 |
} |
467 |
void |
468 |
DataExpanded::symmetric(DataAbstract* ev) |
469 |
{ |
470 |
int sampleNo,dataPointNo; |
471 |
int numSamples = getNumSamples(); |
472 |
int numDataPointsPerSample = getNumDPPSample(); |
473 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
474 |
if (temp_ev==0) { |
475 |
throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error)."); |
476 |
} |
477 |
DataArrayView& thisView=getPointDataView(); |
478 |
DataArrayView& evView=ev->getPointDataView(); |
479 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
480 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
481 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
482 |
DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo), |
483 |
evView,ev->getPointOffset(sampleNo,dataPointNo)); |
484 |
} |
485 |
} |
486 |
} |
487 |
void |
488 |
DataExpanded::nonsymmetric(DataAbstract* ev) |
489 |
{ |
490 |
int sampleNo,dataPointNo; |
491 |
int numSamples = getNumSamples(); |
492 |
int numDataPointsPerSample = getNumDPPSample(); |
493 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
494 |
if (temp_ev==0) { |
495 |
throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error)."); |
496 |
} |
497 |
DataArrayView& thisView=getPointDataView(); |
498 |
DataArrayView& evView=ev->getPointDataView(); |
499 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
500 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
501 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
502 |
DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo), |
503 |
evView,ev->getPointOffset(sampleNo,dataPointNo)); |
504 |
} |
505 |
} |
506 |
} |
507 |
void |
508 |
DataExpanded::trace(DataAbstract* ev, int axis_offset) |
509 |
{ |
510 |
int sampleNo,dataPointNo; |
511 |
int numSamples = getNumSamples(); |
512 |
int numDataPointsPerSample = getNumDPPSample(); |
513 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
514 |
if (temp_ev==0) { |
515 |
throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error)."); |
516 |
} |
517 |
DataArrayView& thisView=getPointDataView(); |
518 |
DataArrayView& evView=ev->getPointDataView(); |
519 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
520 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
521 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
522 |
DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo), |
523 |
evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset); |
524 |
} |
525 |
} |
526 |
} |
527 |
|
528 |
void |
529 |
DataExpanded::transpose(DataAbstract* ev, int axis_offset) |
530 |
{ |
531 |
int sampleNo,dataPointNo; |
532 |
int numSamples = getNumSamples(); |
533 |
int numDataPointsPerSample = getNumDPPSample(); |
534 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
535 |
if (temp_ev==0) { |
536 |
throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error)."); |
537 |
} |
538 |
DataArrayView& thisView=getPointDataView(); |
539 |
DataArrayView& evView=ev->getPointDataView(); |
540 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
541 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
542 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
543 |
DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo), |
544 |
evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset); |
545 |
} |
546 |
} |
547 |
} |
548 |
|
549 |
void |
550 |
DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1) |
551 |
{ |
552 |
int sampleNo,dataPointNo; |
553 |
int numSamples = getNumSamples(); |
554 |
int numDataPointsPerSample = getNumDPPSample(); |
555 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
556 |
if (temp_ev==0) { |
557 |
throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error)."); |
558 |
} |
559 |
DataArrayView& thisView=getPointDataView(); |
560 |
DataArrayView& evView=ev->getPointDataView(); |
561 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
562 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
563 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
564 |
DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo), |
565 |
evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1); |
566 |
} |
567 |
} |
568 |
} |
569 |
void |
570 |
DataExpanded::eigenvalues(DataAbstract* ev) |
571 |
{ |
572 |
int sampleNo,dataPointNo; |
573 |
int numSamples = getNumSamples(); |
574 |
int numDataPointsPerSample = getNumDPPSample(); |
575 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
576 |
if (temp_ev==0) { |
577 |
throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error)."); |
578 |
} |
579 |
DataArrayView& thisView=getPointDataView(); |
580 |
DataArrayView& evView=ev->getPointDataView(); |
581 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
582 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
583 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
584 |
DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo), |
585 |
evView,ev->getPointOffset(sampleNo,dataPointNo)); |
586 |
} |
587 |
} |
588 |
} |
589 |
void |
590 |
DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol) |
591 |
{ |
592 |
int numSamples = getNumSamples(); |
593 |
int numDataPointsPerSample = getNumDPPSample(); |
594 |
int sampleNo,dataPointNo; |
595 |
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
596 |
if (temp_ev==0) { |
597 |
throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error)."); |
598 |
} |
599 |
DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V); |
600 |
if (temp_V==0) { |
601 |
throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error)."); |
602 |
} |
603 |
DataArrayView& thisView=getPointDataView(); |
604 |
DataArrayView& evView=ev->getPointDataView(); |
605 |
DataArrayView& VView=V->getPointDataView(); |
606 |
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
607 |
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
608 |
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
609 |
DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo), |
610 |
evView,ev->getPointOffset(sampleNo,dataPointNo), |
611 |
VView,V->getPointOffset(sampleNo,dataPointNo), |
612 |
tol); |
613 |
} |
614 |
} |
615 |
} |
616 |
|
617 |
void |
618 |
DataExpanded::dump(const std::string fileName) const |
619 |
{ |
620 |
#ifdef PASO_MPI |
621 |
throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.") |
622 |
#endif |
623 |
const NcDim* ncdims[2+DataArrayView::maxRank]; |
624 |
NcVar* var; |
625 |
int rank = getPointDataView().getRank(); |
626 |
int type= getFunctionSpace().getTypeCode(); |
627 |
int ndims =0; |
628 |
long dims[2+DataArrayView::maxRank]; |
629 |
DataArrayView::ShapeType shape = getPointDataView().getShape(); |
630 |
|
631 |
// netCDF error handler |
632 |
NcError err(NcError::verbose_nonfatal); |
633 |
// Create the file. |
634 |
NcFile dataFile(fileName.c_str(), NcFile::Replace); |
635 |
// check if writing was successful |
636 |
if (!dataFile.is_valid()) |
637 |
throw DataException("Error - DataExpanded:: opening of netCDF file for output failed."); |
638 |
if (!dataFile.add_att("type","constant") ) |
639 |
throw DataException("Error - DataExpanded:: appending data type to netCDF file failed."); |
640 |
if (!dataFile.add_att("rank",rank) ) |
641 |
throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed."); |
642 |
if (!dataFile.add_att("function_space_type",type)) |
643 |
throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed."); |
644 |
ndims=rank+2; |
645 |
if ( rank >1 ) { |
646 |
dims[0]=shape[0]; |
647 |
if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) ) |
648 |
throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed."); |
649 |
} |
650 |
if ( rank >1 ) { |
651 |
dims[1]=shape[1]; |
652 |
if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) ) |
653 |
throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed."); |
654 |
} |
655 |
if ( rank >2 ) { |
656 |
dims[2]=shape[2]; |
657 |
if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) ) |
658 |
throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed."); |
659 |
} |
660 |
if ( rank >3 ) { |
661 |
dims[3]=shape[3]; |
662 |
if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) ) |
663 |
throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed."); |
664 |
} |
665 |
dims[rank]= getNumDataPointsPerSample(); |
666 |
if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) ) |
667 |
throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed."); |
668 |
dims[rank+1]= getNumSamples(); |
669 |
if (! (ncdims[rank+1] = dataFile.add_dim("num_sample", dims[rank+1])) ) |
670 |
throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed."); |
671 |
if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) ) |
672 |
throw DataException("Error - DataExpanded:: appending variable to netCDF file failed."); |
673 |
if (! (var->put(&m_data[0],dims)) ) |
674 |
throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed."); |
675 |
} |
676 |
|
677 |
|
678 |
} // end of namespace |