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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1846 - (hide annotations)
Fri Oct 3 06:42:35 2008 UTC (10 years, 10 months ago) by gross
File size: 31596 byte(s)
bug fixed.
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 jfenwick 1796 ValueType& vec=getVector();
443 gross 921 if (dataPointRank==0) {
444 jfenwick 1796 vec[0]=extract<double>(value[0]);
445 gross 921 } else if (dataPointRank==1) {
446 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
447     vec[i]=extract<double>(value[i]);
448 gross 921 }
449     } else if (dataPointRank==2) {
450 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
451     for (int j=0; j<shape[1]; j++) {
452     vec[getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
453 gross 921 }
454     }
455     } else if (dataPointRank==3) {
456 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
457     for (int j=0; j<shape[1]; j++) {
458     for (int k=0; k<shape[2]; k++) {
459     vec[getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
460 gross 921 }
461     }
462     }
463     } else if (dataPointRank==4) {
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     for (int l=0; l<shape[3]; l++) {
468     vec[getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
469 gross 921 }
470     }
471     }
472     }
473 matt 1319 }
474 gross 921 }
475     }
476 jgs 126 void
477     DataExpanded::copyAll(const boost::python::numeric::array& value) {
478     //
479     // Get the number of samples and data-points per sample.
480     int numSamples = getNumSamples();
481     int numDataPointsPerSample = getNumDPPSample();
482 jfenwick 1796 int dataPointRank = getRank();
483     const ShapeType& dataPointShape = getShape();
484 jgs 126 //
485     // check rank:
486     if (value.getrank()!=dataPointRank+1)
487     throw DataException("Rank of numarray does not match Data object rank");
488     if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
489     throw DataException("leading dimension of numarray is too small");
490     //
491 jfenwick 1796 ValueType& vec=getVector();
492 jgs 126 int dataPoint = 0;
493     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
494     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
495 jfenwick 1796 ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
496 jgs 126 if (dataPointRank==0) {
497 jfenwick 1796 vec[offset]=extract<double>(value[dataPoint]);
498 jgs 126 } else if (dataPointRank==1) {
499     for (int i=0; i<dataPointShape[0]; i++) {
500 jfenwick 1796 vec[offset+i]=extract<double>(value[dataPoint][i]);
501 jgs 126 }
502     } else if (dataPointRank==2) {
503     for (int i=0; i<dataPointShape[0]; i++) {
504     for (int j=0; j<dataPointShape[1]; j++) {
505 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
506 jgs 126 }
507     }
508     } else if (dataPointRank==3) {
509     for (int i=0; i<dataPointShape[0]; i++) {
510     for (int j=0; j<dataPointShape[1]; j++) {
511     for (int k=0; k<dataPointShape[2]; k++) {
512 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
513 jgs 126 }
514     }
515     }
516     } else if (dataPointRank==4) {
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     for (int l=0; l<dataPointShape[3]; l++) {
521 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
522 jgs 126 }
523     }
524     }
525     }
526     }
527     dataPoint++;
528     }
529     }
530     }
531 gross 580 void
532 ksteube 775 DataExpanded::symmetric(DataAbstract* ev)
533     {
534     int sampleNo,dataPointNo;
535     int numSamples = getNumSamples();
536     int numDataPointsPerSample = getNumDPPSample();
537     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
538     if (temp_ev==0) {
539     throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
540     }
541 jfenwick 1796 ValueType& vec=getVector();
542     const ShapeType& shape=getShape();
543     ValueType& evVec=temp_ev->getVector();
544     const ShapeType& evShape=temp_ev->getShape();
545 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
546     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
547     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
548 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
549     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
550 ksteube 775 }
551     }
552     }
553     void
554     DataExpanded::nonsymmetric(DataAbstract* ev)
555     {
556     int sampleNo,dataPointNo;
557     int numSamples = getNumSamples();
558     int numDataPointsPerSample = getNumDPPSample();
559     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
560     if (temp_ev==0) {
561     throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
562     }
563 jfenwick 1796 ValueType& vec=getVector();
564     const ShapeType& shape=getShape();
565     ValueType& evVec=temp_ev->getVector();
566     const ShapeType& evShape=temp_ev->getShape();
567 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
568     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
569     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
570 jfenwick 1796 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
571     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
572 ksteube 775 }
573     }
574     }
575     void
576 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
577 ksteube 775 {
578     int sampleNo,dataPointNo;
579     int numSamples = getNumSamples();
580     int numDataPointsPerSample = getNumDPPSample();
581     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
582     if (temp_ev==0) {
583 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
584 ksteube 775 }
585 jfenwick 1796 ValueType& vec=getVector();
586     const ShapeType& shape=getShape();
587     ValueType& evVec=temp_ev->getVector();
588     const ShapeType& evShape=temp_ev->getShape();
589 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
590     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
591     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
592 jfenwick 1796 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
593     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
594 ksteube 775 }
595     }
596     }
597 gross 800
598 ksteube 775 void
599     DataExpanded::transpose(DataAbstract* ev, int axis_offset)
600     {
601     int sampleNo,dataPointNo;
602     int numSamples = getNumSamples();
603     int numDataPointsPerSample = getNumDPPSample();
604     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
605     if (temp_ev==0) {
606     throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
607     }
608 jfenwick 1796 ValueType& vec=getVector();
609     const ShapeType& shape=getShape();
610     ValueType& evVec=temp_ev->getVector();
611     const ShapeType& evShape=temp_ev->getShape();
612 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
613     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
614     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
615 jfenwick 1796 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
616     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
617 ksteube 775 }
618     }
619     }
620 gross 800
621 ksteube 775 void
622 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
623 gross 800 {
624     int sampleNo,dataPointNo;
625     int numSamples = getNumSamples();
626     int numDataPointsPerSample = getNumDPPSample();
627     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
628     if (temp_ev==0) {
629 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
630 gross 800 }
631 jfenwick 1796 ValueType& vec=getVector();
632     const ShapeType& shape=getShape();
633     ValueType& evVec=temp_ev->getVector();
634     const ShapeType& evShape=temp_ev->getShape();
635 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
636     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
637     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
638 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
639     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
640 gross 800 }
641     }
642     }
643     void
644 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
645     {
646 gross 584 int sampleNo,dataPointNo;
647 gross 580 int numSamples = getNumSamples();
648     int numDataPointsPerSample = getNumDPPSample();
649     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
650     if (temp_ev==0) {
651     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
652     }
653 jfenwick 1796 ValueType& vec=getVector();
654     const ShapeType& shape=getShape();
655     ValueType& evVec=temp_ev->getVector();
656     const ShapeType& evShape=temp_ev->getShape();
657 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
658 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
659     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
660 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
661     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
662 gross 580 }
663     }
664     }
665     void
666     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
667     {
668     int numSamples = getNumSamples();
669     int numDataPointsPerSample = getNumDPPSample();
670 gross 584 int sampleNo,dataPointNo;
671 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
672     if (temp_ev==0) {
673     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
674     }
675     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
676     if (temp_V==0) {
677     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
678     }
679 jfenwick 1796 ValueType& vec=getVector();
680     const ShapeType& shape=getShape();
681     ValueType& evVec=temp_ev->getVector();
682     const ShapeType& evShape=temp_ev->getShape();
683     ValueType& VVec=temp_V->getVector();
684     const ShapeType& VShape=temp_V->getShape();
685 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
686 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
687     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
688 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
689     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
690     VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
691 gross 580 }
692     }
693     }
694    
695 gross 950 void
696 gross 1118 DataExpanded::setToZero(){
697 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
698     // Why is there no memset here? Parallel issues?
699 gross 1118 int numSamples = getNumSamples();
700     int numDataPointsPerSample = getNumDPPSample();
701 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
702 gross 1118 double* p;
703     int sampleNo,dataPointNo, i;
704     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
705     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
706     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
707     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
708 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
709 gross 1118 }
710     }
711     }
712    
713 ksteube 1748 /* Append MPI rank to file name if multiple MPI processes */
714     char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
715     /* Make plenty of room for the mpi_rank number and terminating '\0' */
716     char *newFileName = (char *)malloc(strlen(fileName)+20);
717     strncpy(newFileName, fileName, strlen(fileName)+1);
718     if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
719     return(newFileName);
720     }
721 gross 1118
722     void
723 gross 950 DataExpanded::dump(const std::string fileName) const
724     {
725 gross 1023 #ifdef USE_NETCDF
726 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
727 gross 967 const NcDim* ncdims[ldims];
728     NcVar *var, *ids;
729 jfenwick 1796 int rank = getRank();
730 gross 967 int type= getFunctionSpace().getTypeCode();
731     int ndims =0;
732     long dims[ldims];
733 gross 1141 const double* d_ptr=&(m_data[0]);
734 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
735 ksteube 1800 int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
736     int mpi_num=getFunctionSpace().getDomain().getMPISize();
737     #ifdef PASO_MPI
738     MPI_Status status;
739     #endif
740 gross 965
741 ksteube 1827 #ifdef PASO_MPI
742     /* Serialize NetCDF I/O */
743     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
744     #endif
745    
746 gross 967 // netCDF error handler
747     NcError err(NcError::verbose_nonfatal);
748     // Create the file.
749 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
750     NcFile dataFile(newFileName, NcFile::Replace);
751 gross 967 // check if writing was successful
752     if (!dataFile.is_valid())
753     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
754 gross 1141 if (!dataFile.add_att("type_id",2) )
755 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
756     if (!dataFile.add_att("rank",rank) )
757     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
758     if (!dataFile.add_att("function_space_type",type))
759     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
760     ndims=rank+2;
761     if ( rank >0 ) {
762     dims[0]=shape[0];
763     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
764 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
765 gross 967 }
766     if ( rank >1 ) {
767     dims[1]=shape[1];
768     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
769 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
770 gross 967 }
771     if ( rank >2 ) {
772     dims[2]=shape[2];
773     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
774 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
775 gross 967 }
776     if ( rank >3 ) {
777     dims[3]=shape[3];
778     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
779 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
780 gross 967 }
781     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
782     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
783     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
784     dims[rank+1]=getFunctionSpace().getNumSamples();
785     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
786     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
787    
788 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
789     {
790    
791     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
792 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
793 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
794     if (! (ids->put(ids_p,dims[rank+1])) )
795 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
796 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
797 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
798 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
799 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
800 jfenwick 1808 }
801 ksteube 1800 #ifdef PASO_MPI
802 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
803 ksteube 1800 #endif
804 gross 1023 #else
805     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
806     #endif
807 gross 950 }
808 gross 1358
809 jfenwick 1796 void
810 gross 1358 DataExpanded::setTaggedValue(int tagKey,
811 jfenwick 1796 const DataTypes::ShapeType& pointshape,
812     const DataTypes::ValueType& value,
813     int dataOffset)
814 gross 1358 {
815     int numSamples = getNumSamples();
816     int numDataPointsPerSample = getNumDPPSample();
817     int sampleNo,dataPointNo, i;
818 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
819     double* p;
820     const double* in=&value[0+dataOffset];
821 gross 1358
822 jfenwick 1796 if (value.size() != n) {
823 gross 1358 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
824     }
825    
826     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
827     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
828     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
829     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
830     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
831 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
832 gross 1358 }
833     }
834     }
835     }
836    
837 jfenwick 1796
838 gross 1487 void
839     DataExpanded::reorderByReferenceIDs(int *reference_ids)
840     {
841     int numSamples = getNumSamples();
842 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
843 gross 1487 int sampleNo, sampleNo2,i;
844     double* p,*p2;
845     register double rtmp;
846     FunctionSpace fs=getFunctionSpace();
847 gross 1358
848 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
849     const int id_in=reference_ids[sampleNo];
850     const int id=fs.getReferenceIDOfSample(sampleNo);
851     if (id!=id_in) {
852     bool matched=false;
853     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
854     if (id == reference_ids[sampleNo2]) {
855     p=&(m_data[getPointOffset(sampleNo,0)]);
856     p2=&(m_data[getPointOffset(sampleNo2,0)]);
857 gross 1513 for (i=0; i<n ;i++) {
858 gross 1487 rtmp=p[i];
859     p[i]=p2[i];
860     p2[i]=rtmp;
861     }
862     reference_ids[sampleNo]=id;
863     reference_ids[sampleNo2]=id_in;
864     matched=true;
865     break;
866     }
867     }
868 phornby 1628 if (! matched) {
869 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
870     }
871     }
872     }
873     }
874    
875 jfenwick 1796 DataTypes::ValueType&
876     DataExpanded::getVector()
877     {
878     return m_data.getData();
879     }
880 gross 1487
881 jfenwick 1796 const DataTypes::ValueType&
882     DataExpanded::getVector() const
883     {
884     return m_data.getData();
885     }
886    
887 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