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

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

Parent Directory Parent Directory | Revision Log Revision Log


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