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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1801 - (show annotations)
Fri Sep 19 01:37:09 2008 UTC (11 years, 1 month ago) by ksteube
File size: 34327 byte(s)
Fixed serialization of I/O for MPI...code didn't compile without MPI

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26