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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1808 - (hide annotations)
Thu Sep 25 03:14:56 2008 UTC (10 years, 10 months ago) by jfenwick
File size: 31595 byte(s)
Removed some commented out lines.
Modified DataExpanded to not throw when creating objects with zero 
samples.
Modified toString() to report "(data contains no samples)" rather than 
printing a blank line.
Modified DataExpanded::dump() and load so that they do not attempt so 
save/load the ids and data fields if the data object contains no 
samples.



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