1 |
// $Id$ |
|
2 |
/* |
/* $Id$ */ |
3 |
************************************************************ |
|
4 |
* Copyright 2006 by ACcESS MNRF * |
/******************************************************* |
5 |
* * |
* |
6 |
* http://www.access.edu.au * |
* Copyright 2003-2007 by ACceSS MNRF |
7 |
* Primary Business: Queensland, Australia * |
* Copyright 2007 by University of Queensland |
8 |
* Licensed under the Open Software License version 3.0 * |
* |
9 |
* http://www.opensource.org/licenses/osl-3.0.php * |
* http://esscc.uq.edu.au |
10 |
* * |
* Primary Business: Queensland, Australia |
11 |
************************************************************ |
* Licensed under the Open Software License version 3.0 |
12 |
*/ |
* http://www.opensource.org/licenses/osl-3.0.php |
13 |
|
* |
14 |
|
*******************************************************/ |
15 |
|
|
16 |
#include "DataExpanded.h" |
#include "DataExpanded.h" |
17 |
#include "DataException.h" |
#include "DataException.h" |
18 |
#include "DataConstant.h" |
#include "DataConstant.h" |
19 |
#include "DataTagged.h" |
#include "DataTagged.h" |
20 |
|
#ifdef USE_NETCDF |
21 |
|
#include <netcdfcpp.h> |
22 |
|
#endif |
23 |
|
|
24 |
#include <boost/python/extract.hpp> |
#include <boost/python/extract.hpp> |
25 |
|
|
159 |
{ |
{ |
160 |
} |
} |
161 |
|
|
|
void |
|
|
DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape) |
|
|
{ |
|
|
if (getPointDataView().getRank()!=0) { |
|
|
stringstream temp; |
|
|
temp << "Error - Can only reshape Data with data points of rank 0. " |
|
|
<< "This Data has data points with rank: " |
|
|
<< getPointDataView().getRank(); |
|
|
throw DataException(temp.str()); |
|
|
} |
|
|
// |
|
|
// create the new DataBlocks2D data array, and a corresponding DataArrayView |
|
|
DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape)); |
|
|
DataArrayView newView(newData.getData(),shape); |
|
|
// |
|
|
// Copy the original data to every value for the new shape |
|
|
int i,j; |
|
|
int nRows=m_data.getNumRows(); |
|
|
int nCols=m_data.getNumCols(); |
|
|
#pragma omp parallel for private(i,j) schedule(static) |
|
|
for (i=0;i<nRows;i++) { |
|
|
for (j=0;j<nCols;j++) { |
|
|
// NOTE: An exception may be thown from this call if |
|
|
// DOASSERT is on which of course will play |
|
|
// havoc with the omp threads. Run single threaded |
|
|
// if using DOASSERT. |
|
|
newView.copy(newData.index(i,j),m_data(i,j)); |
|
|
} |
|
|
} |
|
|
// swap the new data array for the original |
|
|
m_data.Swap(newData); |
|
|
// set the corresponding DataArrayView |
|
|
DataArrayView temp(m_data.getData(),shape); |
|
|
setPointDataView(temp); |
|
|
} |
|
|
|
|
162 |
DataAbstract* |
DataAbstract* |
163 |
DataExpanded::getSlice(const DataArrayView::RegionType& region) const |
DataExpanded::getSlice(const DataArrayView::RegionType& region) const |
164 |
{ |
{ |
165 |
return new DataExpanded(*this,region); |
return new DataExpanded(*this,region); |
166 |
} |
} |
167 |
|
|
168 |
void |
void |
169 |
DataExpanded::setSlice(const DataAbstract* value, |
DataExpanded::setSlice(const DataAbstract* value, |
170 |
const DataArrayView::RegionType& region) |
const DataArrayView::RegionType& region) |
171 |
{ |
{ |
172 |
const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value); |
const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value); |
173 |
if (tempDataExp==0) { |
if (tempDataExp==0) { |
203 |
} |
} |
204 |
|
|
205 |
void |
void |
206 |
DataExpanded::copy(const DataArrayView& value) |
DataExpanded::copy(const DataArrayView& value) |
207 |
{ |
{ |
208 |
// |
// |
209 |
// copy a single value to every data point in this object |
// copy a single value to every data point in this object |
213 |
#pragma omp parallel for private(i,j) schedule(static) |
#pragma omp parallel for private(i,j) schedule(static) |
214 |
for (i=0;i<nRows;i++) { |
for (i=0;i<nRows;i++) { |
215 |
for (j=0;j<nCols;j++) { |
for (j=0;j<nCols;j++) { |
216 |
// NOTE: An exception may be thown from this call if |
// NOTE: An exception may be thown from this call if |
217 |
// DOASSERT is on which of course will play |
// DOASSERT is on which of course will play |
218 |
// havoc with the omp threads. Run single threaded |
// havoc with the omp threads. Run single threaded |
219 |
// if using DOASSERT. |
// if using DOASSERT. |
220 |
getPointDataView().copy(m_data.index(i,j),value); |
getPointDataView().copy(getPointOffset(i,j),value); |
221 |
} |
} |
222 |
} |
} |
223 |
} |
} |
224 |
|
|
225 |
void |
void |
226 |
DataExpanded::copy(const boost::python::numeric::array& value) |
DataExpanded::copy(const boost::python::numeric::array& value) |
227 |
{ |
{ |
228 |
// |
|
229 |
// first convert the numarray into a DataArray object |
// extract the shape of the numarray |
230 |
DataArray temp(value); |
DataArrayView::ShapeType tempShape; |
231 |
|
for (int i=0; i < value.getrank(); i++) { |
232 |
|
tempShape.push_back(extract<int>(value.getshape()[i])); |
233 |
|
} |
234 |
|
|
235 |
|
// get the space for the data vector |
236 |
|
int len = DataArrayView::noValues(tempShape); |
237 |
|
DataVector temp_data(len, 0.0, len); |
238 |
|
DataArrayView temp_dataView(temp_data, tempShape); |
239 |
|
temp_dataView.copy(value); |
240 |
|
|
241 |
// |
// |
242 |
// check the input shape matches this shape |
// check the input shape matches this shape |
243 |
if (!getPointDataView().checkShape(temp.getView().getShape())) { |
if (!getPointDataView().checkShape(temp_dataView.getShape())) { |
244 |
throw DataException(getPointDataView().createShapeErrorMessage( |
throw DataException(getPointDataView().createShapeErrorMessage( |
245 |
"Error - (DataExpanded) Cannot copy due to shape mismatch.", |
"Error - (DataExpanded) Cannot copy due to shape mismatch.", |
246 |
temp.getView().getShape())); |
temp_dataView.getShape())); |
247 |
} |
} |
248 |
// |
// |
249 |
// now copy over the data |
// now copy over the data |
250 |
copy(temp.getView()); |
copy(temp_dataView); |
251 |
} |
} |
252 |
|
|
253 |
void |
void |
268 |
DataExpanded::toString() const |
DataExpanded::toString() const |
269 |
{ |
{ |
270 |
stringstream temp; |
stringstream temp; |
271 |
|
FunctionSpace fs=getFunctionSpace(); |
272 |
// |
// |
273 |
// create a temporary view as the offset will be changed |
// create a temporary view as the offset will be changed |
274 |
DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset()); |
DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset()); |
275 |
for (int i=0;i<m_data.getNumRows();i++) { |
for (int i=0;i<m_data.getNumRows();i++) { |
276 |
for (int j=0;j<m_data.getNumCols();j++) { |
for (int j=0;j<m_data.getNumCols();j++) { |
277 |
tempView.setOffset(m_data.index(i,j)); |
tempView.setOffset(getPointOffset(i,j)); |
278 |
stringstream suffix; |
stringstream suffix; |
279 |
suffix << "(" << i << "," << j << ")"; |
suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")"; |
280 |
temp << tempView.toString(suffix.str()); |
temp << tempView.toString(suffix.str()); |
281 |
if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) { |
if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) { |
282 |
temp << endl; |
temp << endl; |
297 |
DataExpanded::getDataPoint(int sampleNo, |
DataExpanded::getDataPoint(int sampleNo, |
298 |
int dataPointNo) |
int dataPointNo) |
299 |
{ |
{ |
300 |
DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo)); |
DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo)); |
301 |
return temp; |
return temp; |
302 |
} |
} |
303 |
|
|
307 |
return m_data.size(); |
return m_data.size(); |
308 |
} |
} |
309 |
|
|
|
void |
|
|
DataExpanded::setRefValue(int ref, |
|
|
const DataArray& value) |
|
|
{ |
|
|
// |
|
|
// Get the number of samples and data-points per sample. |
|
|
int numSamples = getNumSamples(); |
|
|
int numDPPSample = getNumDPPSample(); |
|
|
|
|
|
// |
|
|
// Determine the sample number which corresponds to this reference number. |
|
|
int sampleNo = -1; |
|
|
int tempRef = -1; |
|
|
for (int n=0; n<numSamples; n++) { |
|
|
tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n); |
|
|
if (tempRef == ref) { |
|
|
sampleNo = n; |
|
|
break; |
|
|
} |
|
|
} |
|
|
if (sampleNo == -1) { |
|
|
throw DataException("DataExpanded::setRefValue error: invalid ref number supplied."); |
|
|
} |
|
|
|
|
|
for (int n=0; n<numDPPSample; n++) { |
|
|
// |
|
|
// Get *each* data-point in the sample in turn. |
|
|
DataArrayView pointView = getDataPoint(sampleNo, n); |
|
|
// |
|
|
// Assign the values in the DataArray to this data-point. |
|
|
pointView.copy(value.getView()); |
|
|
} |
|
|
} |
|
|
|
|
|
void |
|
|
DataExpanded::getRefValue(int ref, |
|
|
DataArray& value) |
|
|
{ |
|
|
// |
|
|
// Get the number of samples and data-points per sample. |
|
|
int numSamples = getNumSamples(); |
|
|
int numDPPSample = getNumDPPSample(); |
|
|
|
|
|
// |
|
|
// Determine the sample number which corresponds to this reference number |
|
|
int sampleNo = -1; |
|
|
int tempRef = -1; |
|
|
for (int n=0; n<numSamples; n++) { |
|
|
tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n); |
|
|
if (tempRef == ref) { |
|
|
sampleNo = n; |
|
|
break; |
|
|
} |
|
|
} |
|
|
if (sampleNo == -1) { |
|
|
throw DataException("DataExpanded::getRefValue error: invalid ref number supplied."); |
|
|
} |
|
|
|
|
|
// |
|
|
// Get the *first* data-point associated with this sample number. |
|
|
DataArrayView pointView = getDataPoint(sampleNo, 0); |
|
|
|
|
|
// |
|
|
// Load the values from this data-point into the DataArray |
|
|
value.getView().copy(pointView); |
|
|
} |
|
|
|
|
310 |
int |
int |
311 |
DataExpanded::archiveData(ofstream& archiveFile, |
DataExpanded::archiveData(ofstream& archiveFile, |
312 |
const DataArrayView::ValueType::size_type noValues) const |
const DataArrayView::ValueType::size_type noValues) const |
322 |
} |
} |
323 |
|
|
324 |
void |
void |
325 |
|
DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) { |
326 |
|
// |
327 |
|
// Get the number of samples and data-points per sample. |
328 |
|
int numSamples = getNumSamples(); |
329 |
|
int numDataPointsPerSample = getNumDPPSample(); |
330 |
|
int dataPointRank = getPointDataView().getRank(); |
331 |
|
ShapeType dataPointShape = getPointDataView().getShape(); |
332 |
|
if (numSamples*numDataPointsPerSample > 0) { |
333 |
|
//TODO: global error handling |
334 |
|
if ((sampleNo >= numSamples) || (sampleNo < 0 )) { |
335 |
|
throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo."); |
336 |
|
} |
337 |
|
if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) { |
338 |
|
throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample."); |
339 |
|
} |
340 |
|
DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo); |
341 |
|
if (dataPointRank==0) { |
342 |
|
dataPointView()=value; |
343 |
|
} else if (dataPointRank==1) { |
344 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
345 |
|
dataPointView(i)=value; |
346 |
|
} |
347 |
|
} else if (dataPointRank==2) { |
348 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
349 |
|
for (int j=0; j<dataPointShape[1]; j++) { |
350 |
|
dataPointView(i,j)=value; |
351 |
|
} |
352 |
|
} |
353 |
|
} else if (dataPointRank==3) { |
354 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
355 |
|
for (int j=0; j<dataPointShape[1]; j++) { |
356 |
|
for (int k=0; k<dataPointShape[2]; k++) { |
357 |
|
dataPointView(i,j,k)=value; |
358 |
|
} |
359 |
|
} |
360 |
|
} |
361 |
|
} else if (dataPointRank==4) { |
362 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
363 |
|
for (int j=0; j<dataPointShape[1]; j++) { |
364 |
|
for (int k=0; k<dataPointShape[2]; k++) { |
365 |
|
for (int l=0; l<dataPointShape[3]; l++) { |
366 |
|
dataPointView(i,j,k,l)=value; |
367 |
|
} |
368 |
|
} |
369 |
|
} |
370 |
|
} |
371 |
|
} |
372 |
|
} |
373 |
|
} |
374 |
|
void |
375 |
|
DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) { |
376 |
|
// |
377 |
|
// Get the number of samples and data-points per sample. |
378 |
|
int numSamples = getNumSamples(); |
379 |
|
int numDataPointsPerSample = getNumDPPSample(); |
380 |
|
int dataPointRank = getPointDataView().getRank(); |
381 |
|
ShapeType dataPointShape = getPointDataView().getShape(); |
382 |
|
// |
383 |
|
// check rank: |
384 |
|
if (value.getrank()!=dataPointRank) |
385 |
|
throw DataException("Rank of numarray does not match Data object rank"); |
386 |
|
if (numSamples*numDataPointsPerSample > 0) { |
387 |
|
//TODO: global error handling |
388 |
|
if ((sampleNo >= numSamples) || (sampleNo < 0 )) { |
389 |
|
throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo."); |
390 |
|
} |
391 |
|
if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) { |
392 |
|
throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample."); |
393 |
|
} |
394 |
|
DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo); |
395 |
|
if (dataPointRank==0) { |
396 |
|
dataPointView()=extract<double>(value[0]); |
397 |
|
} else if (dataPointRank==1) { |
398 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
399 |
|
dataPointView(i)=extract<double>(value[i]); |
400 |
|
} |
401 |
|
} else if (dataPointRank==2) { |
402 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
403 |
|
for (int j=0; j<dataPointShape[1]; j++) { |
404 |
|
dataPointView(i,j)=extract<double>(value[i][j]); |
405 |
|
} |
406 |
|
} |
407 |
|
} else if (dataPointRank==3) { |
408 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
409 |
|
for (int j=0; j<dataPointShape[1]; j++) { |
410 |
|
for (int k=0; k<dataPointShape[2]; k++) { |
411 |
|
dataPointView(i,j,k)=extract<double>(value[i][j][k]); |
412 |
|
} |
413 |
|
} |
414 |
|
} |
415 |
|
} else if (dataPointRank==4) { |
416 |
|
for (int i=0; i<dataPointShape[0]; i++) { |
417 |
|
for (int j=0; j<dataPointShape[1]; j++) { |
418 |
|
for (int k=0; k<dataPointShape[2]; k++) { |
419 |
|
for (int l=0; l<dataPointShape[3]; l++) { |
420 |
|
dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]); |
421 |
|
} |
422 |
|
} |
423 |
|
} |
424 |
|
} |
425 |
|
} |
426 |
|
} |
427 |
|
} |
428 |
|
void |
429 |
DataExpanded::copyAll(const boost::python::numeric::array& value) { |
DataExpanded::copyAll(const boost::python::numeric::array& value) { |
430 |
// |
// |
431 |
// Get the number of samples and data-points per sample. |
// Get the number of samples and data-points per sample. |
480 |
} |
} |
481 |
} |
} |
482 |
void |
void |
483 |
|
DataExpanded::symmetric(DataAbstract* ev) |
484 |
|
{ |
485 |
|
int sampleNo,dataPointNo; |
486 |
|
int numSamples = getNumSamples(); |
487 |
|
int numDataPointsPerSample = getNumDPPSample(); |
488 |
|
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
489 |
|
if (temp_ev==0) { |
490 |
|
throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error)."); |
491 |
|
} |
492 |
|
DataArrayView& thisView=getPointDataView(); |
493 |
|
DataArrayView& evView=ev->getPointDataView(); |
494 |
|
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
495 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
496 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
497 |
|
DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo), |
498 |
|
evView,ev->getPointOffset(sampleNo,dataPointNo)); |
499 |
|
} |
500 |
|
} |
501 |
|
} |
502 |
|
void |
503 |
|
DataExpanded::nonsymmetric(DataAbstract* ev) |
504 |
|
{ |
505 |
|
int sampleNo,dataPointNo; |
506 |
|
int numSamples = getNumSamples(); |
507 |
|
int numDataPointsPerSample = getNumDPPSample(); |
508 |
|
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
509 |
|
if (temp_ev==0) { |
510 |
|
throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error)."); |
511 |
|
} |
512 |
|
DataArrayView& thisView=getPointDataView(); |
513 |
|
DataArrayView& evView=ev->getPointDataView(); |
514 |
|
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
515 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
516 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
517 |
|
DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo), |
518 |
|
evView,ev->getPointOffset(sampleNo,dataPointNo)); |
519 |
|
} |
520 |
|
} |
521 |
|
} |
522 |
|
void |
523 |
|
DataExpanded::trace(DataAbstract* ev, int axis_offset) |
524 |
|
{ |
525 |
|
int sampleNo,dataPointNo; |
526 |
|
int numSamples = getNumSamples(); |
527 |
|
int numDataPointsPerSample = getNumDPPSample(); |
528 |
|
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
529 |
|
if (temp_ev==0) { |
530 |
|
throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error)."); |
531 |
|
} |
532 |
|
DataArrayView& thisView=getPointDataView(); |
533 |
|
DataArrayView& evView=ev->getPointDataView(); |
534 |
|
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
535 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
536 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
537 |
|
DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo), |
538 |
|
evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset); |
539 |
|
} |
540 |
|
} |
541 |
|
} |
542 |
|
|
543 |
|
void |
544 |
|
DataExpanded::transpose(DataAbstract* ev, int axis_offset) |
545 |
|
{ |
546 |
|
int sampleNo,dataPointNo; |
547 |
|
int numSamples = getNumSamples(); |
548 |
|
int numDataPointsPerSample = getNumDPPSample(); |
549 |
|
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
550 |
|
if (temp_ev==0) { |
551 |
|
throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error)."); |
552 |
|
} |
553 |
|
DataArrayView& thisView=getPointDataView(); |
554 |
|
DataArrayView& evView=ev->getPointDataView(); |
555 |
|
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
556 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
557 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
558 |
|
DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo), |
559 |
|
evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset); |
560 |
|
} |
561 |
|
} |
562 |
|
} |
563 |
|
|
564 |
|
void |
565 |
|
DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1) |
566 |
|
{ |
567 |
|
int sampleNo,dataPointNo; |
568 |
|
int numSamples = getNumSamples(); |
569 |
|
int numDataPointsPerSample = getNumDPPSample(); |
570 |
|
DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev); |
571 |
|
if (temp_ev==0) { |
572 |
|
throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error)."); |
573 |
|
} |
574 |
|
DataArrayView& thisView=getPointDataView(); |
575 |
|
DataArrayView& evView=ev->getPointDataView(); |
576 |
|
#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static) |
577 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
578 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
579 |
|
DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo), |
580 |
|
evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1); |
581 |
|
} |
582 |
|
} |
583 |
|
} |
584 |
|
void |
585 |
DataExpanded::eigenvalues(DataAbstract* ev) |
DataExpanded::eigenvalues(DataAbstract* ev) |
586 |
{ |
{ |
587 |
int sampleNo,dataPointNo; |
int sampleNo,dataPointNo; |
629 |
} |
} |
630 |
} |
} |
631 |
|
|
632 |
|
void |
633 |
|
DataExpanded::setToZero(){ |
634 |
|
int numSamples = getNumSamples(); |
635 |
|
int numDataPointsPerSample = getNumDPPSample(); |
636 |
|
DataArrayView& thisView=getPointDataView(); |
637 |
|
DataArrayView::ValueType::size_type n = thisView.noValues(); |
638 |
|
double* p; |
639 |
|
int sampleNo,dataPointNo, i; |
640 |
|
#pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static) |
641 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
642 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
643 |
|
p=&(m_data[getPointOffset(sampleNo,dataPointNo)]); |
644 |
|
for (int i=0; i<n ;++i) p[i]=0.; |
645 |
|
} |
646 |
|
} |
647 |
|
} |
648 |
|
|
649 |
|
|
650 |
|
void |
651 |
|
DataExpanded::dump(const std::string fileName) const |
652 |
|
{ |
653 |
|
#ifdef PASO_MPI |
654 |
|
throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet."); |
655 |
|
#endif |
656 |
|
#ifdef USE_NETCDF |
657 |
|
const int ldims=2+DataArrayView::maxRank; |
658 |
|
const NcDim* ncdims[ldims]; |
659 |
|
NcVar *var, *ids; |
660 |
|
int rank = getPointDataView().getRank(); |
661 |
|
int type= getFunctionSpace().getTypeCode(); |
662 |
|
int ndims =0; |
663 |
|
long dims[ldims]; |
664 |
|
const double* d_ptr=&(m_data[0]); |
665 |
|
DataArrayView::ShapeType shape = getPointDataView().getShape(); |
666 |
|
|
667 |
|
// netCDF error handler |
668 |
|
NcError err(NcError::verbose_nonfatal); |
669 |
|
// Create the file. |
670 |
|
NcFile dataFile(fileName.c_str(), NcFile::Replace); |
671 |
|
// check if writing was successful |
672 |
|
if (!dataFile.is_valid()) |
673 |
|
throw DataException("Error - DataExpanded:: opening of netCDF file for output failed."); |
674 |
|
if (!dataFile.add_att("type_id",2) ) |
675 |
|
throw DataException("Error - DataExpanded:: appending data type to netCDF file failed."); |
676 |
|
if (!dataFile.add_att("rank",rank) ) |
677 |
|
throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed."); |
678 |
|
if (!dataFile.add_att("function_space_type",type)) |
679 |
|
throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed."); |
680 |
|
ndims=rank+2; |
681 |
|
if ( rank >0 ) { |
682 |
|
dims[0]=shape[0]; |
683 |
|
if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) ) |
684 |
|
throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed."); |
685 |
|
} |
686 |
|
if ( rank >1 ) { |
687 |
|
dims[1]=shape[1]; |
688 |
|
if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) ) |
689 |
|
throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed."); |
690 |
|
} |
691 |
|
if ( rank >2 ) { |
692 |
|
dims[2]=shape[2]; |
693 |
|
if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) ) |
694 |
|
throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed."); |
695 |
|
} |
696 |
|
if ( rank >3 ) { |
697 |
|
dims[3]=shape[3]; |
698 |
|
if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) ) |
699 |
|
throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed."); |
700 |
|
} |
701 |
|
dims[rank]=getFunctionSpace().getNumDataPointsPerSample(); |
702 |
|
if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) ) |
703 |
|
throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed."); |
704 |
|
dims[rank+1]=getFunctionSpace().getNumSamples(); |
705 |
|
if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) ) |
706 |
|
throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed."); |
707 |
|
|
708 |
|
if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) ) |
709 |
|
throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed."); |
710 |
|
const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs(); |
711 |
|
if (! (ids->put(ids_p,dims[rank+1])) ) |
712 |
|
throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed."); |
713 |
|
|
714 |
|
if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) ) |
715 |
|
throw DataException("Error - DataExpanded:: appending variable to netCDF file failed."); |
716 |
|
if (! (var->put(d_ptr,dims)) ) |
717 |
|
throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed."); |
718 |
|
#else |
719 |
|
throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager."); |
720 |
|
#endif |
721 |
|
} |
722 |
|
|
723 |
|
void |
724 |
|
DataExpanded::setTaggedValue(int tagKey, |
725 |
|
const DataArrayView& value) |
726 |
|
{ |
727 |
|
int numSamples = getNumSamples(); |
728 |
|
int numDataPointsPerSample = getNumDPPSample(); |
729 |
|
int sampleNo,dataPointNo, i; |
730 |
|
DataArrayView& thisView=getPointDataView(); |
731 |
|
DataArrayView::ValueType::size_type n = thisView.noValues(); |
732 |
|
double* p,*in=&(value.getData()[0]); |
733 |
|
|
734 |
|
if (value.noValues() != n) { |
735 |
|
throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points."); |
736 |
|
} |
737 |
|
|
738 |
|
#pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static) |
739 |
|
for (sampleNo = 0; sampleNo < numSamples; sampleNo++) { |
740 |
|
if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) { |
741 |
|
for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) { |
742 |
|
p=&(m_data[getPointOffset(sampleNo,dataPointNo)]); |
743 |
|
for (int i=0; i<n ;++i) p[i]=in[i]; |
744 |
|
} |
745 |
|
} |
746 |
|
} |
747 |
|
} |
748 |
|
|
749 |
|
|
750 |
} // end of namespace |
} // end of namespace |