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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1827 - (show annotations)
Thu Oct 2 04:28:07 2008 UTC (10 years, 10 months ago) by ksteube
File size: 31591 byte(s)
MPI parallelized dump/load for Constant and Tagged data.
run_escriptOnFinley.py runs hybrid MPI/OpenMP.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26