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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1799 - (show annotations)
Wed Sep 17 06:33:18 2008 UTC (11 years ago) by jfenwick
File size: 34049 byte(s)
Added Data::copySelf() [Note: this is exposed as copy() in python].
This method returns a pointer to a deep copy of the target.
There are c++ tests but no python tests for this yet.

All DataAbstracts now have a deepCopy() which simplifies the 
implementation of the compy methods.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26