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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (11 years ago) by jfenwick
File size: 31840 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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 : parent(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 : parent(other.getFunctionSpace(), other.getShape()),
51 m_data(other.m_data)
52 {
53 }
54
55 DataExpanded::DataExpanded(const DataConstant& other)
56 : parent(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 : parent(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 : parent(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 : parent(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::getPointOffset(int sampleNo,
365 int dataPointNo)
366 {
367 return m_data.index(sampleNo,dataPointNo);
368 }
369
370 DataTypes::ValueType::size_type
371 DataExpanded::getLength() const
372 {
373 return m_data.size();
374 }
375
376
377
378 void
379 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
380 //
381 // Get the number of samples and data-points per sample.
382 int numSamples = getNumSamples();
383 int numDataPointsPerSample = getNumDPPSample();
384 int dataPointRank = getRank();
385 ShapeType dataPointShape = getShape();
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 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
395 ValueType& vec=getVector();
396 if (dataPointRank==0) {
397 vec[offset]=value;
398 } else if (dataPointRank==1) {
399 for (int i=0; i<dataPointShape[0]; i++) {
400 vec[offset+i]=value;
401 }
402 } else if (dataPointRank==2) {
403 for (int i=0; i<dataPointShape[0]; i++) {
404 for (int j=0; j<dataPointShape[1]; j++) {
405 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
406 }
407 }
408 } else if (dataPointRank==3) {
409 for (int i=0; i<dataPointShape[0]; i++) {
410 for (int j=0; j<dataPointShape[1]; j++) {
411 for (int k=0; k<dataPointShape[2]; k++) {
412 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
413 }
414 }
415 }
416 } else if (dataPointRank==4) {
417 for (int i=0; i<dataPointShape[0]; i++) {
418 for (int j=0; j<dataPointShape[1]; j++) {
419 for (int k=0; k<dataPointShape[2]; k++) {
420 for (int l=0; l<dataPointShape[3]; l++) {
421 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
422 }
423 }
424 }
425 }
426 }
427 }
428 }
429 void
430 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
431 //
432 // Get the number of samples and data-points per sample.
433 int numSamples = getNumSamples();
434 int numDataPointsPerSample = getNumDPPSample();
435 int dataPointRank = getRank();
436 const ShapeType& shape = getShape();
437 //
438 // check rank:
439 if (value.getrank()!=dataPointRank)
440 throw DataException("Rank of numarray does not match Data object rank");
441 if (numSamples*numDataPointsPerSample > 0) {
442 //TODO: global error handling
443 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
444 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
445 }
446 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
447 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
448 }
449 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
450 ValueType& vec=getVector();
451 if (dataPointRank==0) {
452 vec[offset]=extract<double>(value[0]);
453 } else if (dataPointRank==1) {
454 for (int i=0; i<shape[0]; i++) {
455 vec[offset+i]=extract<double>(value[i]);
456 }
457 } else if (dataPointRank==2) {
458 for (int i=0; i<shape[0]; i++) {
459 for (int j=0; j<shape[1]; j++) {
460 vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
461 }
462 }
463 } else if (dataPointRank==3) {
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 vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
468 }
469 }
470 }
471 } else if (dataPointRank==4) {
472 for (int i=0; i<shape[0]; i++) {
473 for (int j=0; j<shape[1]; j++) {
474 for (int k=0; k<shape[2]; k++) {
475 for (int l=0; l<shape[3]; l++) {
476 vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
477 }
478 }
479 }
480 }
481 }
482 }
483 }
484 void
485 DataExpanded::copyAll(const boost::python::numeric::array& value) {
486 //
487 // Get the number of samples and data-points per sample.
488 int numSamples = getNumSamples();
489 int numDataPointsPerSample = getNumDPPSample();
490 int dataPointRank = getRank();
491 const ShapeType& dataPointShape = getShape();
492 //
493 // check rank:
494 if (value.getrank()!=dataPointRank+1)
495 throw DataException("Rank of numarray does not match Data object rank");
496 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
497 throw DataException("leading dimension of numarray is too small");
498 //
499 ValueType& vec=getVector();
500 int dataPoint = 0;
501 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503 ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
504 if (dataPointRank==0) {
505 vec[offset]=extract<double>(value[dataPoint]);
506 } else if (dataPointRank==1) {
507 for (int i=0; i<dataPointShape[0]; i++) {
508 vec[offset+i]=extract<double>(value[dataPoint][i]);
509 }
510 } else if (dataPointRank==2) {
511 for (int i=0; i<dataPointShape[0]; i++) {
512 for (int j=0; j<dataPointShape[1]; j++) {
513 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
514 }
515 }
516 } else if (dataPointRank==3) {
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 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
521 }
522 }
523 }
524 } else if (dataPointRank==4) {
525 for (int i=0; i<dataPointShape[0]; i++) {
526 for (int j=0; j<dataPointShape[1]; j++) {
527 for (int k=0; k<dataPointShape[2]; k++) {
528 for (int l=0; l<dataPointShape[3]; l++) {
529 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
530 }
531 }
532 }
533 }
534 }
535 dataPoint++;
536 }
537 }
538 }
539 void
540 DataExpanded::symmetric(DataAbstract* ev)
541 {
542 int sampleNo,dataPointNo;
543 int numSamples = getNumSamples();
544 int numDataPointsPerSample = getNumDPPSample();
545 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
546 if (temp_ev==0) {
547 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
548 }
549 ValueType& vec=getVector();
550 const ShapeType& shape=getShape();
551 ValueType& evVec=temp_ev->getVector();
552 const ShapeType& evShape=temp_ev->getShape();
553 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
554 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
555 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
556 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
557 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
558 }
559 }
560 }
561 void
562 DataExpanded::nonsymmetric(DataAbstract* ev)
563 {
564 int sampleNo,dataPointNo;
565 int numSamples = getNumSamples();
566 int numDataPointsPerSample = getNumDPPSample();
567 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
568 if (temp_ev==0) {
569 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
570 }
571 ValueType& vec=getVector();
572 const ShapeType& shape=getShape();
573 ValueType& evVec=temp_ev->getVector();
574 const ShapeType& evShape=temp_ev->getShape();
575 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
576 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
577 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
578 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
579 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
580 }
581 }
582 }
583 void
584 DataExpanded::trace(DataAbstract* ev, int axis_offset)
585 {
586 int sampleNo,dataPointNo;
587 int numSamples = getNumSamples();
588 int numDataPointsPerSample = getNumDPPSample();
589 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
590 if (temp_ev==0) {
591 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
592 }
593 ValueType& vec=getVector();
594 const ShapeType& shape=getShape();
595 ValueType& evVec=temp_ev->getVector();
596 const ShapeType& evShape=temp_ev->getShape();
597 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
598 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
599 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
600 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
601 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
602 }
603 }
604 }
605
606 void
607 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
608 {
609 int sampleNo,dataPointNo;
610 int numSamples = getNumSamples();
611 int numDataPointsPerSample = getNumDPPSample();
612 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
613 if (temp_ev==0) {
614 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
615 }
616 ValueType& vec=getVector();
617 const ShapeType& shape=getShape();
618 ValueType& evVec=temp_ev->getVector();
619 const ShapeType& evShape=temp_ev->getShape();
620 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
621 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
622 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
623 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
624 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
625 }
626 }
627 }
628
629 void
630 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
631 {
632 int sampleNo,dataPointNo;
633 int numSamples = getNumSamples();
634 int numDataPointsPerSample = getNumDPPSample();
635 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
636 if (temp_ev==0) {
637 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
638 }
639 ValueType& vec=getVector();
640 const ShapeType& shape=getShape();
641 ValueType& evVec=temp_ev->getVector();
642 const ShapeType& evShape=temp_ev->getShape();
643 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
644 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
647 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
648 }
649 }
650 }
651 void
652 DataExpanded::eigenvalues(DataAbstract* ev)
653 {
654 int sampleNo,dataPointNo;
655 int numSamples = getNumSamples();
656 int numDataPointsPerSample = getNumDPPSample();
657 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
658 if (temp_ev==0) {
659 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
660 }
661 ValueType& vec=getVector();
662 const ShapeType& shape=getShape();
663 ValueType& evVec=temp_ev->getVector();
664 const ShapeType& evShape=temp_ev->getShape();
665 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
666 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
667 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
668 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
669 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
670 }
671 }
672 }
673 void
674 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
675 {
676 int numSamples = getNumSamples();
677 int numDataPointsPerSample = getNumDPPSample();
678 int sampleNo,dataPointNo;
679 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
680 if (temp_ev==0) {
681 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
682 }
683 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
684 if (temp_V==0) {
685 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
686 }
687 ValueType& vec=getVector();
688 const ShapeType& shape=getShape();
689 ValueType& evVec=temp_ev->getVector();
690 const ShapeType& evShape=temp_ev->getShape();
691 ValueType& VVec=temp_V->getVector();
692 const ShapeType& VShape=temp_V->getShape();
693 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
694 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
695 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
696 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
697 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
698 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
699 }
700 }
701 }
702
703 void
704 DataExpanded::setToZero(){
705 // TODO: Surely there is a more efficient way to do this????
706 // Why is there no memset here? Parallel issues?
707 int numSamples = getNumSamples();
708 int numDataPointsPerSample = getNumDPPSample();
709 DataTypes::ValueType::size_type n = getNoValues();
710 double* p;
711 int sampleNo,dataPointNo, i;
712 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
713 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
714 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
715 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
716 for (i=0; i<n ;++i) p[i]=0.;
717 }
718 }
719 }
720
721 /* Append MPI rank to file name if multiple MPI processes */
722 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
723 /* Make plenty of room for the mpi_rank number and terminating '\0' */
724 char *newFileName = (char *)malloc(strlen(fileName)+20);
725 strncpy(newFileName, fileName, strlen(fileName)+1);
726 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
727 return(newFileName);
728 }
729
730 void
731 DataExpanded::dump(const std::string fileName) const
732 {
733 #ifdef USE_NETCDF
734 const int ldims=2+DataTypes::maxRank;
735 const NcDim* ncdims[ldims];
736 NcVar *var, *ids;
737 int rank = getRank();
738 int type= getFunctionSpace().getTypeCode();
739 int ndims =0;
740 long dims[ldims];
741 const double* d_ptr=&(m_data[0]);
742 const DataTypes::ShapeType& shape = getShape();
743 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
744 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
745 #ifdef PASO_MPI
746 MPI_Status status;
747 #endif
748
749 #ifdef PASO_MPI
750 /* Serialize NetCDF I/O */
751 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
752 #endif
753
754 // netCDF error handler
755 NcError err(NcError::verbose_nonfatal);
756 // Create the file.
757 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
758 NcFile dataFile(newFileName, NcFile::Replace);
759 // check if writing was successful
760 if (!dataFile.is_valid())
761 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
762 if (!dataFile.add_att("type_id",2) )
763 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
764 if (!dataFile.add_att("rank",rank) )
765 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
766 if (!dataFile.add_att("function_space_type",type))
767 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
768 ndims=rank+2;
769 if ( rank >0 ) {
770 dims[0]=shape[0];
771 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
772 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
773 }
774 if ( rank >1 ) {
775 dims[1]=shape[1];
776 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
777 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
778 }
779 if ( rank >2 ) {
780 dims[2]=shape[2];
781 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
782 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
783 }
784 if ( rank >3 ) {
785 dims[3]=shape[3];
786 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
787 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
788 }
789 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
790 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
791 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
792 dims[rank+1]=getFunctionSpace().getNumSamples();
793 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
794 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
795
796 if (getFunctionSpace().getNumSamples()>0)
797 {
798
799 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
800 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
801 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
802 if (! (ids->put(ids_p,dims[rank+1])) )
803 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
804 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
805 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
806 if (! (var->put(d_ptr,dims)) )
807 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
808 }
809 #ifdef PASO_MPI
810 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
811 #endif
812 #else
813 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
814 #endif
815 }
816
817 void
818 DataExpanded::setTaggedValue(int tagKey,
819 const DataTypes::ShapeType& pointshape,
820 const DataTypes::ValueType& value,
821 int dataOffset)
822 {
823 int numSamples = getNumSamples();
824 int numDataPointsPerSample = getNumDPPSample();
825 int sampleNo,dataPointNo, i;
826 DataTypes::ValueType::size_type n = getNoValues();
827 double* p;
828 const double* in=&value[0+dataOffset];
829
830 if (value.size() != n) {
831 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
832 }
833
834 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
835 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
836 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
837 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
838 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
839 for (i=0; i<n ;++i) p[i]=in[i];
840 }
841 }
842 }
843 }
844
845
846 void
847 DataExpanded::reorderByReferenceIDs(int *reference_ids)
848 {
849 int numSamples = getNumSamples();
850 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
851 int sampleNo, sampleNo2,i;
852 double* p,*p2;
853 register double rtmp;
854 FunctionSpace fs=getFunctionSpace();
855
856 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
857 const int id_in=reference_ids[sampleNo];
858 const int id=fs.getReferenceIDOfSample(sampleNo);
859 if (id!=id_in) {
860 bool matched=false;
861 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
862 if (id == reference_ids[sampleNo2]) {
863 p=&(m_data[getPointOffset(sampleNo,0)]);
864 p2=&(m_data[getPointOffset(sampleNo2,0)]);
865 for (i=0; i<n ;i++) {
866 rtmp=p[i];
867 p[i]=p2[i];
868 p2[i]=rtmp;
869 }
870 reference_ids[sampleNo]=id;
871 reference_ids[sampleNo2]=id_in;
872 matched=true;
873 break;
874 }
875 }
876 if (! matched) {
877 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
878 }
879 }
880 }
881 }
882
883 DataTypes::ValueType&
884 DataExpanded::getVector()
885 {
886 return m_data.getData();
887 }
888
889 const DataTypes::ValueType&
890 DataExpanded::getVector() const
891 {
892 return m_data.getData();
893 }
894
895 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26