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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1800 - (hide annotations)
Thu Sep 18 05:28:20 2008 UTC (10 years, 11 months ago) by ksteube
File size: 34385 byte(s)
Serialized parallel I/O when writing mesh or data to NetCDF file on multiple MPI processors.
Added domain method getMPIComm() to complement getMPISize() and getMPIRank().

1 jgs 82
2 ksteube 1312 /* $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 ksteube 1800 #include "Data.h"
17 jgs 480 #include "DataExpanded.h"
18 jgs 474 #include "DataException.h"
19     #include "DataConstant.h"
20     #include "DataTagged.h"
21 gross 1023 #ifdef USE_NETCDF
22 ksteube 1312 #include <netcdfcpp.h>
23 gross 1023 #endif
24 ksteube 1748 #ifdef PASO_MPI
25     #include <mpi.h>
26     #endif
27 jgs 82
28     #include <boost/python/extract.hpp>
29 jfenwick 1796 #include "DataMaths.h"
30 jgs 82
31     using namespace std;
32     using namespace boost::python;
33     using namespace boost;
34 jfenwick 1796 using namespace escript::DataTypes;
35 jgs 82
36     namespace escript {
37    
38 jgs 102 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
39     const FunctionSpace& what)
40 jfenwick 1796 : DataAbstract(what,DataTypes::shapeFromNumArray(value))
41 jgs 102 {
42 jfenwick 1796 /* DataTypes::ShapeType tempShape;
43 jgs 102 //
44 jgs 117 // extract the shape of the python numarray
45     for (int i=0; i<value.getrank(); i++) {
46 jgs 102 tempShape.push_back(extract<int>(value.getshape()[i]));
47 jfenwick 1796 }*/
48 jgs 117 //
49     // initialise the data array for this object
50 jfenwick 1796 initialise(what.getNumSamples(),what.getNumDPPSample());
51 jgs 102 //
52     // copy the given value to every data point
53     copy(value);
54     }
55 jgs 82
56 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
57 jfenwick 1796 : DataAbstract(other.getFunctionSpace(), other.getShape()),
58 jgs 102 m_data(other.m_data)
59 jgs 117 {
60 jfenwick 1796 /*
61 jgs 102 //
62     // create the view for the data
63     DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
64     setPointDataView(temp);
65 jfenwick 1796
66     // keep local shape in sync
67     setShape(other.getPointDataView().getShape());*/
68 jgs 102 }
69 jgs 82
70 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
71 jfenwick 1796 : DataAbstract(other.getFunctionSpace(), other.getShape())
72 jgs 117 {
73     //
74     // initialise the data array for this object
75 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
76 jgs 102 //
77 jgs 117 // DataConstant only has one value, copy this to every data point
78 jfenwick 1796 copy(other);
79 jgs 102 }
80 jgs 82
81 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
82 jfenwick 1796 : DataAbstract(other.getFunctionSpace(), other.getShape())
83 jgs 117 {
84     //
85 jgs 119 // initialise the data array for this object
86 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
87 jgs 117 //
88     // for each data point in this object, extract and copy the corresponding data
89     // value from the given DataTagged object
90 jgs 102 int i,j;
91 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
92     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
93 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
94 jgs 117 for (i=0;i<numRows;i++) {
95     for (j=0;j<numCols;j++) {
96 jgs 102 try {
97 jfenwick 1796 // 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 jgs 122 other.getPointOffset(i,j));
103 jgs 82 }
104 jgs 102 catch (std::exception& e) {
105     cout << e.what() << endl;
106     }
107 jgs 82 }
108     }
109 jgs 102 }
110 jgs 82
111 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
112 jfenwick 1796 const DataTypes::RegionType& region)
113     : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
114 jgs 102 {
115     //
116     // get the shape of the slice
117 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
118 jgs 102 //
119     // initialise this Data object to the shape of the slice
120 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
121 jgs 102 //
122 jgs 117 // copy the data
123 jfenwick 1796 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 jgs 102 int i,j;
127 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
128 jgs 117 for (i=0;i<numRows;i++) {
129     for (j=0;j<numCols;j++) {
130 jgs 102 try {
131 jfenwick 1796 // 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 jgs 122 other.getPointOffset(i,j),
139     region_loop_range);
140 jgs 82 }
141 jgs 102 catch (std::exception& e) {
142     cout << e.what() << endl;
143     }
144 jgs 82 }
145     }
146 jgs 102 }
147 jgs 82
148 jfenwick 1796 // 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 jgs 82
163 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
164 jfenwick 1796 const DataTypes::ShapeType &shape,
165     const DataTypes::ValueType &data)
166     : DataAbstract(what,shape)
167 jgs 119 {
168 jfenwick 1796 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 jgs 119 }
195    
196 jgs 102 DataExpanded::~DataExpanded()
197     {
198     }
199 jgs 82
200 jgs 102 DataAbstract*
201 jfenwick 1799 DataExpanded::deepCopy()
202     {
203     return new DataExpanded(*this);
204     }
205    
206    
207     DataAbstract*
208 jfenwick 1796 DataExpanded::getSlice(const DataTypes::RegionType& region) const
209 jgs 102 {
210     return new DataExpanded(*this,region);
211     }
212    
213     void
214     DataExpanded::setSlice(const DataAbstract* value,
215 jfenwick 1796 const DataTypes::RegionType& region)
216 jgs 102 {
217     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
218     if (tempDataExp==0) {
219     throw DataException("Programming error - casting to DataExpanded.");
220     }
221 jgs 108 //
222 jgs 117 // get shape of slice
223 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
224     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
225 jgs 108 //
226 jgs 117 // check shape
227 jfenwick 1796 if (getRank()!=region.size()) {
228 jgs 102 throw DataException("Error - Invalid slice region.");
229     }
230 jfenwick 1796 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 jgs 102 }
234     //
235 jgs 117 // copy the data from the slice into this object
236 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
237     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
238 jgs 102 int i, j;
239 jfenwick 1796 ValueType& vec=getVector();
240     const ShapeType& mshape=getShape();
241     const ValueType& tVec=tempDataExp->getVector();
242     const ShapeType& tShape=tempDataExp->getShape();
243 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
244 jgs 102 for (i=0;i<numRows;i++) {
245     for (j=0;j<numCols;j++) {
246 jfenwick 1796 /* getPointDataView().copySliceFrom(getPointOffset(i,j),
247 jgs 122 tempDataExp->getPointDataView(),
248     tempDataExp->getPointOffset(i,j),
249 jfenwick 1796 region_loop_range);*/
250     DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
251     tVec,
252     tShape,
253     tempDataExp->getPointOffset(i,j),
254 jgs 122 region_loop_range);
255 jfenwick 1796
256 jgs 82 }
257     }
258 jgs 102 }
259 jgs 82
260 jgs 102 void
261 jfenwick 1796 DataExpanded::copy(const DataConstant& value)
262 jgs 102 {
263 jfenwick 1796 EsysAssert((checkShape(getShape(), value.getShape())),
264     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
265    
266 jgs 102 //
267 jgs 117 // copy a single value to every data point in this object
268 jgs 102 int nRows=m_data.getNumRows();
269     int nCols=m_data.getNumCols();
270 jgs 117 int i,j;
271 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
272 jgs 117 for (i=0;i<nRows;i++) {
273     for (j=0;j<nCols;j++) {
274 matt 1319 // NOTE: An exception may be thown from this call if
275 jgs 102 // DOASSERT is on which of course will play
276     // havoc with the omp threads. Run single threaded
277 matt 1319 // if using DOASSERT.
278 jfenwick 1796 //getPointDataView().copy(getPointOffset(i,j),value);
279     DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
280 jgs 82 }
281     }
282 jgs 102 }
283 jgs 82
284 jfenwick 1796
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 jgs 102 void
306 matt 1319 DataExpanded::copy(const boost::python::numeric::array& value)
307 jgs 102 {
308 matt 1319
309     // extract the shape of the numarray
310 jfenwick 1796 DataTypes::ShapeType tempShape;
311 matt 1319 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 jfenwick 1796 // 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 matt 1319
321 jgs 102 //
322 jgs 117 // check the input shape matches this shape
323 jfenwick 1796 if (!DataTypes::checkShape(getShape(),tempShape)) {
324     throw DataException(DataTypes::createShapeErrorMessage(
325 jgs 102 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
326 jfenwick 1796 tempShape,getShape()));
327 jgs 82 }
328 jgs 102 //
329 jgs 117 // now copy over the data
330 jfenwick 1796 //copy(temp_dataView);
331     getVector().copyFromNumArray(value);
332 jgs 102 }
333 jgs 82
334 jfenwick 1796
335 jgs 102 void
336 jfenwick 1796 DataExpanded::initialise(int noSamples,
337 jgs 102 int noDataPointsPerSample)
338     {
339     //
340 jgs 117 // resize data array to the required size
341 jfenwick 1796 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
342    
343 jgs 102 //
344 jfenwick 1796 // // 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 jgs 102 }
351 jgs 82
352 jgs 102 string
353     DataExpanded::toString() const
354     {
355     stringstream temp;
356 gross 856 FunctionSpace fs=getFunctionSpace();
357 jgs 102 //
358     // create a temporary view as the offset will be changed
359 jfenwick 1796 // DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
360    
361     int offset=0;
362 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
363     for (int j=0;j<m_data.getNumCols();j++) {
364 jfenwick 1796 offset=getPointOffset(i,j);
365 jgs 102 stringstream suffix;
366 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
367 jfenwick 1796 temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
368 jgs 102 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
369     temp << endl;
370 jgs 82 }
371     }
372     }
373 ksteube 1312 return temp.str();
374 jgs 102 }
375 jgs 82
376 jfenwick 1796 DataTypes::ValueType::size_type
377 jgs 102 DataExpanded::getPointOffset(int sampleNo,
378     int dataPointNo) const
379     {
380     return m_data.index(sampleNo,dataPointNo);
381     }
382 jgs 82
383 jfenwick 1796 // 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 jgs 82
391 jfenwick 1796 DataTypes::ValueType::size_type
392 jgs 102 DataExpanded::getLength() const
393     {
394 jgs 117 return m_data.size();
395 jgs 102 }
396 jgs 82
397 jgs 123
398    
399 matt 1319 void
400 gross 922 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 jfenwick 1796 int dataPointRank = getRank();
406     ShapeType dataPointShape = getShape();
407 gross 922 if (numSamples*numDataPointsPerSample > 0) {
408     //TODO: global error handling
409 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
410 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
411     }
412 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
413 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
414     }
415 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
416     ValueType& vec=getVector();
417 gross 922 if (dataPointRank==0) {
418 jfenwick 1796 vec[0]=value;
419 gross 922 } else if (dataPointRank==1) {
420     for (int i=0; i<dataPointShape[0]; i++) {
421 jfenwick 1796 vec[offset+i]=value;
422 gross 922 }
423     } else if (dataPointRank==2) {
424     for (int i=0; i<dataPointShape[0]; i++) {
425     for (int j=0; j<dataPointShape[1]; j++) {
426 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
427 gross 922 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
434 gross 922 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
443 gross 922 }
444     }
445     }
446     }
447 matt 1319 }
448 gross 922 }
449     }
450 matt 1319 void
451 gross 921 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 jfenwick 1796 int dataPointRank = getRank();
457     const ShapeType& shape = getShape();
458 gross 921 //
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 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
465 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
466     }
467 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
468 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
469     }
470 jfenwick 1796 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
471     ValueType& vec=getVector();
472 gross 921 if (dataPointRank==0) {
473 jfenwick 1796 vec[0]=extract<double>(value[0]);
474 gross 921 } else if (dataPointRank==1) {
475 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
476     vec[i]=extract<double>(value[i]);
477 gross 921 }
478     } else if (dataPointRank==2) {
479 jfenwick 1796 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 gross 921 }
483     }
484     } else if (dataPointRank==3) {
485 jfenwick 1796 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 gross 921 }
490     }
491     }
492     } else if (dataPointRank==4) {
493 jfenwick 1796 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 gross 921 }
499     }
500     }
501     }
502 matt 1319 }
503 gross 921 }
504     }
505 jgs 126 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 jfenwick 1796 int dataPointRank = getRank();
512     const ShapeType& dataPointShape = getShape();
513 jgs 126 //
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 jfenwick 1796 ValueType& vec=getVector();
521 jgs 126 int dataPoint = 0;
522     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
523     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
524 jfenwick 1796 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
525     ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
526 jgs 126 if (dataPointRank==0) {
527 jfenwick 1796 vec[offset]=extract<double>(value[dataPoint]);
528 jgs 126 } else if (dataPointRank==1) {
529     for (int i=0; i<dataPointShape[0]; i++) {
530 jfenwick 1796 vec[offset+i]=extract<double>(value[dataPoint][i]);
531 jgs 126 }
532     } else if (dataPointRank==2) {
533     for (int i=0; i<dataPointShape[0]; i++) {
534     for (int j=0; j<dataPointShape[1]; j++) {
535 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
536     // dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
537 jgs 126 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
544 jgs 126 }
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 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
553 jgs 126 }
554     }
555     }
556     }
557     }
558     dataPoint++;
559     }
560     }
561     }
562 gross 580 void
563 ksteube 775 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 jfenwick 1796 // 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 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
579     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
580     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
581 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
582     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
583 ksteube 775 }
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 jfenwick 1796 // 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 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
603     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
604     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
605 jfenwick 1796 // 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 ksteube 775 }
610     }
611     }
612     void
613 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
614 ksteube 775 {
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 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
621 ksteube 775 }
622 jfenwick 1796 ValueType& vec=getVector();
623     const ShapeType& shape=getShape();
624     ValueType& evVec=temp_ev->getVector();
625     const ShapeType& evShape=temp_ev->getShape();
626 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
627     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
628     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
629 jfenwick 1796 /* 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 ksteube 775 }
634     }
635     }
636 gross 800
637 ksteube 775 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 jfenwick 1796 ValueType& vec=getVector();
648     const ShapeType& shape=getShape();
649     ValueType& evVec=temp_ev->getVector();
650     const ShapeType& evShape=temp_ev->getShape();
651 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
652     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
653     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
654 jfenwick 1796 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
655     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
656 ksteube 775 }
657     }
658     }
659 gross 800
660 ksteube 775 void
661 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
662 gross 800 {
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 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
669 gross 800 }
670 jfenwick 1796 ValueType& vec=getVector();
671     const ShapeType& shape=getShape();
672     ValueType& evVec=temp_ev->getVector();
673     const ShapeType& evShape=temp_ev->getShape();
674 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
675     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
676     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
677 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
678     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
679 gross 800 }
680     }
681     }
682     void
683 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
684     {
685 gross 584 int sampleNo,dataPointNo;
686 gross 580 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 jfenwick 1796 ValueType& vec=getVector();
693     const ShapeType& shape=getShape();
694     ValueType& evVec=temp_ev->getVector();
695     const ShapeType& evShape=temp_ev->getShape();
696 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
697 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
698     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
699 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
700     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
701 gross 580 }
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 gross 584 int sampleNo,dataPointNo;
710 gross 580 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 jfenwick 1796 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 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
725 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
726     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
727 jfenwick 1796 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 gross 580 }
731     }
732     }
733    
734 gross 950 void
735 gross 1118 DataExpanded::setToZero(){
736 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
737     // Why is there no memset here? Parallel issues?
738 gross 1118 int numSamples = getNumSamples();
739     int numDataPointsPerSample = getNumDPPSample();
740 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
741 gross 1118 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 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
748 gross 1118 }
749     }
750     }
751    
752 ksteube 1748 /* 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 gross 1118
761     void
762 gross 950 DataExpanded::dump(const std::string fileName) const
763     {
764 gross 1023 #ifdef USE_NETCDF
765 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
766 gross 967 const NcDim* ncdims[ldims];
767     NcVar *var, *ids;
768 jfenwick 1796 int rank = getRank();
769 gross 967 int type= getFunctionSpace().getTypeCode();
770     int ndims =0;
771     long dims[ldims];
772 gross 1141 const double* d_ptr=&(m_data[0]);
773 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
774 ksteube 1800 int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
775     int mpi_num=getFunctionSpace().getDomain().getMPISize();
776     #ifdef PASO_MPI
777     MPI_Status status;
778     #endif
779 gross 965
780 gross 967 // netCDF error handler
781     NcError err(NcError::verbose_nonfatal);
782     // Create the file.
783 ksteube 1748 #ifdef PASO_MPI
784 ksteube 1800 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, getFunctionSpace().getDomain().getMPIComm(), &status);
785 ksteube 1748 #endif
786     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
787     NcFile dataFile(newFileName, NcFile::Replace);
788 gross 967 // check if writing was successful
789     if (!dataFile.is_valid())
790     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
791 gross 1141 if (!dataFile.add_att("type_id",2) )
792 gross 967 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 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
802 gross 967 }
803     if ( rank >1 ) {
804     dims[1]=shape[1];
805     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
806 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
807 gross 967 }
808     if ( rank >2 ) {
809     dims[2]=shape[2];
810     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
811 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
812 gross 967 }
813     if ( rank >3 ) {
814     dims[3]=shape[3];
815     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
816 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
817 gross 967 }
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 gross 1141 if (! (var->put(d_ptr,dims)) )
834 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
835 ksteube 1800 #ifdef PASO_MPI
836     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, getFunctionSpace().getDomain().getMPIComm());
837     #endif
838 gross 1023 #else
839     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
840     #endif
841 gross 950 }
842 gross 1358
843 jfenwick 1796 // 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 gross 1358 DataExpanded::setTaggedValue(int tagKey,
870 jfenwick 1796 const DataTypes::ShapeType& pointshape,
871     const DataTypes::ValueType& value,
872     int dataOffset)
873 gross 1358 {
874     int numSamples = getNumSamples();
875     int numDataPointsPerSample = getNumDPPSample();
876     int sampleNo,dataPointNo, i;
877 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
878     double* p;
879     const double* in=&value[0+dataOffset];
880 gross 1358
881 jfenwick 1796 if (value.size() != n) {
882 gross 1358 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 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
891 gross 1358 }
892     }
893     }
894     }
895    
896 jfenwick 1796
897 gross 1487 void
898     DataExpanded::reorderByReferenceIDs(int *reference_ids)
899     {
900     int numSamples = getNumSamples();
901 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
902 gross 1487 int sampleNo, sampleNo2,i;
903     double* p,*p2;
904     register double rtmp;
905     FunctionSpace fs=getFunctionSpace();
906 gross 1358
907 gross 1487 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 gross 1513 for (i=0; i<n ;i++) {
917 gross 1487 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 phornby 1628 if (! matched) {
928 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
929     }
930     }
931     }
932     }
933    
934 jfenwick 1796 DataTypes::ValueType&
935     DataExpanded::getVector()
936     {
937     return m_data.getData();
938     }
939 gross 1487
940 jfenwick 1796 const DataTypes::ValueType&
941     DataExpanded::getVector() const
942     {
943     return m_data.getData();
944     }
945    
946 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26