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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1962 - (hide annotations)
Tue Nov 4 07:17:31 2008 UTC (11 years ago) by gross
File size: 31705 byte(s)
there is an offet missing when copy data to a data point.


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 1796 : DataAbstract(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 1796 : DataAbstract(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 1796 : DataAbstract(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 1796 : DataAbstract(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     : DataAbstract(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     : DataAbstract(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     for (int j=0;j<getNoValues();++j,++i)
163     {
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     getVector().copyFromNumArray(value);
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 jgs 102 DataExpanded::getLength() const
365     {
366 jgs 117 return m_data.size();
367 jgs 102 }
368 jgs 82
369 jgs 123
370    
371 matt 1319 void
372 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
373     //
374     // Get the number of samples and data-points per sample.
375     int numSamples = getNumSamples();
376     int numDataPointsPerSample = getNumDPPSample();
377 jfenwick 1796 int dataPointRank = getRank();
378     ShapeType dataPointShape = getShape();
379 gross 922 if (numSamples*numDataPointsPerSample > 0) {
380     //TODO: global error handling
381 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
382 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
383     }
384 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
385 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
386     }
387 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
388     ValueType& vec=getVector();
389 gross 922 if (dataPointRank==0) {
390 gross 1846 vec[offset]=value;
391 gross 922 } else if (dataPointRank==1) {
392     for (int i=0; i<dataPointShape[0]; i++) {
393 jfenwick 1796 vec[offset+i]=value;
394 gross 922 }
395     } else if (dataPointRank==2) {
396     for (int i=0; i<dataPointShape[0]; i++) {
397     for (int j=0; j<dataPointShape[1]; j++) {
398 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
399 gross 922 }
400     }
401     } else if (dataPointRank==3) {
402     for (int i=0; i<dataPointShape[0]; i++) {
403     for (int j=0; j<dataPointShape[1]; j++) {
404     for (int k=0; k<dataPointShape[2]; k++) {
405 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
406 gross 922 }
407     }
408     }
409     } else if (dataPointRank==4) {
410     for (int i=0; i<dataPointShape[0]; i++) {
411     for (int j=0; j<dataPointShape[1]; j++) {
412     for (int k=0; k<dataPointShape[2]; k++) {
413     for (int l=0; l<dataPointShape[3]; l++) {
414 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
415 gross 922 }
416     }
417     }
418     }
419 matt 1319 }
420 gross 922 }
421     }
422 matt 1319 void
423 gross 921 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
424     //
425     // Get the number of samples and data-points per sample.
426     int numSamples = getNumSamples();
427     int numDataPointsPerSample = getNumDPPSample();
428 jfenwick 1796 int dataPointRank = getRank();
429     const ShapeType& shape = getShape();
430 gross 921 //
431     // check rank:
432     if (value.getrank()!=dataPointRank)
433     throw DataException("Rank of numarray does not match Data object rank");
434     if (numSamples*numDataPointsPerSample > 0) {
435     //TODO: global error handling
436 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
437 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
438     }
439 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
440 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
441     }
442 gross 1962 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
443 jfenwick 1796 ValueType& vec=getVector();
444 gross 921 if (dataPointRank==0) {
445 gross 1962 vec[offset]=extract<double>(value[0]);
446 gross 921 } else if (dataPointRank==1) {
447 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
448 gross 1962 vec[offset+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 gross 1962 vec[offset+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 gross 1962 vec[offset+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 gross 1962 vec[offset+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 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
737     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
738 ksteube 1800 #ifdef PASO_MPI
739     MPI_Status status;
740     #endif
741 gross 965
742 ksteube 1827 #ifdef PASO_MPI
743     /* Serialize NetCDF I/O */
744     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
745     #endif
746    
747 gross 967 // netCDF error handler
748     NcError err(NcError::verbose_nonfatal);
749     // Create the file.
750 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
751     NcFile dataFile(newFileName, NcFile::Replace);
752 gross 967 // check if writing was successful
753     if (!dataFile.is_valid())
754     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
755 gross 1141 if (!dataFile.add_att("type_id",2) )
756 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
757     if (!dataFile.add_att("rank",rank) )
758     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
759     if (!dataFile.add_att("function_space_type",type))
760     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
761     ndims=rank+2;
762     if ( rank >0 ) {
763     dims[0]=shape[0];
764     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
765 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
766 gross 967 }
767     if ( rank >1 ) {
768     dims[1]=shape[1];
769     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
770 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
771 gross 967 }
772     if ( rank >2 ) {
773     dims[2]=shape[2];
774     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
775 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
776 gross 967 }
777     if ( rank >3 ) {
778     dims[3]=shape[3];
779     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
780 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
781 gross 967 }
782     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
783     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
784     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
785     dims[rank+1]=getFunctionSpace().getNumSamples();
786     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
787     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
788    
789 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
790     {
791    
792     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
793 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
794 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
795     if (! (ids->put(ids_p,dims[rank+1])) )
796 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
797 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
798 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
799 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
800 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
801 jfenwick 1808 }
802 ksteube 1800 #ifdef PASO_MPI
803 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
804 ksteube 1800 #endif
805 gross 1023 #else
806     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
807     #endif
808 gross 950 }
809 gross 1358
810 jfenwick 1796 void
811 gross 1358 DataExpanded::setTaggedValue(int tagKey,
812 jfenwick 1796 const DataTypes::ShapeType& pointshape,
813     const DataTypes::ValueType& value,
814     int dataOffset)
815 gross 1358 {
816     int numSamples = getNumSamples();
817     int numDataPointsPerSample = getNumDPPSample();
818     int sampleNo,dataPointNo, i;
819 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
820     double* p;
821     const double* in=&value[0+dataOffset];
822 gross 1358
823 jfenwick 1796 if (value.size() != n) {
824 gross 1358 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
825     }
826    
827     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
828     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
829     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
830     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
831     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
832 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
833 gross 1358 }
834     }
835     }
836     }
837    
838 jfenwick 1796
839 gross 1487 void
840     DataExpanded::reorderByReferenceIDs(int *reference_ids)
841     {
842     int numSamples = getNumSamples();
843 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
844 gross 1487 int sampleNo, sampleNo2,i;
845     double* p,*p2;
846     register double rtmp;
847     FunctionSpace fs=getFunctionSpace();
848 gross 1358
849 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
850     const int id_in=reference_ids[sampleNo];
851     const int id=fs.getReferenceIDOfSample(sampleNo);
852     if (id!=id_in) {
853     bool matched=false;
854     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
855     if (id == reference_ids[sampleNo2]) {
856     p=&(m_data[getPointOffset(sampleNo,0)]);
857     p2=&(m_data[getPointOffset(sampleNo2,0)]);
858 gross 1513 for (i=0; i<n ;i++) {
859 gross 1487 rtmp=p[i];
860     p[i]=p2[i];
861     p2[i]=rtmp;
862     }
863     reference_ids[sampleNo]=id;
864     reference_ids[sampleNo2]=id_in;
865     matched=true;
866     break;
867     }
868     }
869 phornby 1628 if (! matched) {
870 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
871     }
872     }
873     }
874     }
875    
876 jfenwick 1796 DataTypes::ValueType&
877     DataExpanded::getVector()
878     {
879     return m_data.getData();
880     }
881 gross 1487
882 jfenwick 1796 const DataTypes::ValueType&
883     DataExpanded::getVector() const
884     {
885     return m_data.getData();
886     }
887    
888 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