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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1808 - (show annotations)
Thu Sep 25 03:14:56 2008 UTC (11 years, 1 month ago) by jfenwick
File size: 31595 byte(s)
Removed some commented out lines.
Modified DataExpanded to not throw when creating objects with zero 
samples.
Modified toString() to report "(data contains no samples)" rather than 
printing a blank line.
Modified DataExpanded::dump() and load so that they do not attempt so 
save/load the ids and data fields if the data object contains no 
samples.



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26