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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (hide annotations)
Mon Nov 10 01:21:39 2008 UTC (10 years, 11 months ago) by jfenwick
File size: 31840 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 ksteube 1800 #include "Data.h"
16 jgs 480 #include "DataExpanded.h"
17 jgs 474 #include "DataException.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20 gross 1023 #ifdef USE_NETCDF
21 ksteube 1312 #include <netcdfcpp.h>
22 gross 1023 #endif
23 ksteube 1748 #ifdef PASO_MPI
24     #include <mpi.h>
25     #endif
26 jgs 82
27     #include <boost/python/extract.hpp>
28 jfenwick 1796 #include "DataMaths.h"
29 jgs 82
30     using namespace std;
31     using namespace boost::python;
32     using namespace boost;
33 jfenwick 1796 using namespace escript::DataTypes;
34 jgs 82
35     namespace escript {
36    
37 jgs 102 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
38     const FunctionSpace& what)
39 jfenwick 2005 : parent(what,DataTypes::shapeFromNumArray(value))
40 jgs 102 {
41     //
42 jgs 117 // initialise the data array for this object
43 jfenwick 1796 initialise(what.getNumSamples(),what.getNumDPPSample());
44 jgs 102 //
45     // copy the given value to every data point
46     copy(value);
47     }
48 jgs 82
49 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
50 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape()),
51 jgs 102 m_data(other.m_data)
52 jgs 117 {
53 jgs 102 }
54 jgs 82
55 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
56 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape())
57 jgs 117 {
58     //
59     // initialise the data array for this object
60 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
61 jgs 102 //
62 jgs 117 // DataConstant only has one value, copy this to every data point
63 jfenwick 1796 copy(other);
64 jgs 102 }
65 jgs 82
66 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
67 jfenwick 2005 : parent(other.getFunctionSpace(), other.getShape())
68 jgs 117 {
69     //
70 jgs 119 // initialise the data array for this object
71 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
72 jgs 117 //
73     // for each data point in this object, extract and copy the corresponding data
74     // value from the given DataTagged object
75 jgs 102 int i,j;
76 jfenwick 1796 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
77     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
78 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
79 jgs 117 for (i=0;i<numRows;i++) {
80     for (j=0;j<numCols;j++) {
81 jgs 102 try {
82 jfenwick 1796 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
83     other.getVector(),
84 jgs 122 other.getPointOffset(i,j));
85 jgs 82 }
86 jgs 102 catch (std::exception& e) {
87     cout << e.what() << endl;
88     }
89 jgs 82 }
90     }
91 jgs 102 }
92 jgs 82
93 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
94 jfenwick 1796 const DataTypes::RegionType& region)
95 jfenwick 2005 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
96 jgs 102 {
97     //
98     // get the shape of the slice
99 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
100 jgs 102 //
101     // initialise this Data object to the shape of the slice
102 jfenwick 1796 initialise(other.getNumSamples(),other.getNumDPPSample());
103 jgs 102 //
104 jgs 117 // copy the data
105 jfenwick 1796 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
106     DataTypes::ValueType::size_type numRows=m_data.getNumRows();
107     DataTypes::ValueType::size_type numCols=m_data.getNumCols();
108 jgs 102 int i,j;
109 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
110 jgs 117 for (i=0;i<numRows;i++) {
111     for (j=0;j<numCols;j++) {
112 jgs 102 try {
113 jfenwick 1796 // getPointDataView().copySlice(getPointOffset(i,j),
114     // other.getPointDataView(),
115     // other.getPointOffset(i,j),
116     // region_loop_range);
117     DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
118     other.getVector(),
119     other.getShape(),
120 jgs 122 other.getPointOffset(i,j),
121     region_loop_range);
122 jgs 82 }
123 jgs 102 catch (std::exception& e) {
124     cout << e.what() << endl;
125     }
126 jgs 82 }
127     }
128 jgs 102 }
129 jgs 82
130 jfenwick 1796 // DataExpanded::DataExpanded(const DataArrayView& value,
131     // const FunctionSpace& what)
132     // : DataAbstract(what)
133     // {
134     // //
135     // // get the shape of the given data value
136     // DataTypes::ShapeType tempShape=value.getShape();
137     // //
138     // // initialise this Data object to the shape of the given data value
139     // initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
140     // //
141     // // copy the given value to every data point
142     // copy(value);
143     // }
144 jgs 82
145 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
146 jfenwick 1796 const DataTypes::ShapeType &shape,
147     const DataTypes::ValueType &data)
148 jfenwick 2005 : parent(what,shape)
149 jgs 119 {
150 jfenwick 1796 EsysAssert(data.size()%getNoValues()==0,
151     "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
152    
153     if (data.size()==getNoValues())
154     {
155     ValueType& vec=m_data.getData();
156     //
157     // create the view of the data
158     initialise(what.getNumSamples(),what.getNumDPPSample());
159     // now we copy this value to all elements
160     for (int i=0;i<getLength();)
161     {
162     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 jfenwick 2005 DataExpanded::getPointOffset(int sampleNo,
365     int dataPointNo)
366     {
367     return m_data.index(sampleNo,dataPointNo);
368     }
369    
370     DataTypes::ValueType::size_type
371 jgs 102 DataExpanded::getLength() const
372     {
373 jgs 117 return m_data.size();
374 jgs 102 }
375 jgs 82
376 jgs 123
377    
378 matt 1319 void
379 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
380     //
381     // Get the number of samples and data-points per sample.
382     int numSamples = getNumSamples();
383     int numDataPointsPerSample = getNumDPPSample();
384 jfenwick 1796 int dataPointRank = getRank();
385     ShapeType dataPointShape = getShape();
386 gross 922 if (numSamples*numDataPointsPerSample > 0) {
387     //TODO: global error handling
388 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
389 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
390     }
391 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
392 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
393     }
394 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
395     ValueType& vec=getVector();
396 gross 922 if (dataPointRank==0) {
397 gross 1846 vec[offset]=value;
398 gross 922 } else if (dataPointRank==1) {
399     for (int i=0; i<dataPointShape[0]; i++) {
400 jfenwick 1796 vec[offset+i]=value;
401 gross 922 }
402     } else if (dataPointRank==2) {
403     for (int i=0; i<dataPointShape[0]; i++) {
404     for (int j=0; j<dataPointShape[1]; j++) {
405 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
406 gross 922 }
407     }
408     } else if (dataPointRank==3) {
409     for (int i=0; i<dataPointShape[0]; i++) {
410     for (int j=0; j<dataPointShape[1]; j++) {
411     for (int k=0; k<dataPointShape[2]; k++) {
412 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
413 gross 922 }
414     }
415     }
416     } else if (dataPointRank==4) {
417     for (int i=0; i<dataPointShape[0]; i++) {
418     for (int j=0; j<dataPointShape[1]; j++) {
419     for (int k=0; k<dataPointShape[2]; k++) {
420     for (int l=0; l<dataPointShape[3]; l++) {
421 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
422 gross 922 }
423     }
424     }
425     }
426 matt 1319 }
427 gross 922 }
428     }
429 matt 1319 void
430 gross 921 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
431     //
432     // Get the number of samples and data-points per sample.
433     int numSamples = getNumSamples();
434     int numDataPointsPerSample = getNumDPPSample();
435 jfenwick 1796 int dataPointRank = getRank();
436     const ShapeType& shape = getShape();
437 gross 921 //
438     // check rank:
439     if (value.getrank()!=dataPointRank)
440     throw DataException("Rank of numarray does not match Data object rank");
441     if (numSamples*numDataPointsPerSample > 0) {
442     //TODO: global error handling
443 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
444 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
445     }
446 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
447 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
448     }
449 gross 1962 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
450 jfenwick 1796 ValueType& vec=getVector();
451 gross 921 if (dataPointRank==0) {
452 gross 1962 vec[offset]=extract<double>(value[0]);
453 gross 921 } else if (dataPointRank==1) {
454 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
455 gross 1962 vec[offset+i]=extract<double>(value[i]);
456 gross 921 }
457     } else if (dataPointRank==2) {
458 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
459     for (int j=0; j<shape[1]; j++) {
460 gross 1962 vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
461 gross 921 }
462     }
463     } else if (dataPointRank==3) {
464 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
465     for (int j=0; j<shape[1]; j++) {
466     for (int k=0; k<shape[2]; k++) {
467 gross 1962 vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
468 gross 921 }
469     }
470     }
471     } else if (dataPointRank==4) {
472 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
473     for (int j=0; j<shape[1]; j++) {
474     for (int k=0; k<shape[2]; k++) {
475     for (int l=0; l<shape[3]; l++) {
476 gross 1962 vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
477 gross 921 }
478     }
479     }
480     }
481 matt 1319 }
482 gross 921 }
483     }
484 jgs 126 void
485     DataExpanded::copyAll(const boost::python::numeric::array& value) {
486     //
487     // Get the number of samples and data-points per sample.
488     int numSamples = getNumSamples();
489     int numDataPointsPerSample = getNumDPPSample();
490 jfenwick 1796 int dataPointRank = getRank();
491     const ShapeType& dataPointShape = getShape();
492 jgs 126 //
493     // check rank:
494     if (value.getrank()!=dataPointRank+1)
495     throw DataException("Rank of numarray does not match Data object rank");
496     if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
497     throw DataException("leading dimension of numarray is too small");
498     //
499 jfenwick 1796 ValueType& vec=getVector();
500 jgs 126 int dataPoint = 0;
501     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503 jfenwick 1796 ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
504 jgs 126 if (dataPointRank==0) {
505 jfenwick 1796 vec[offset]=extract<double>(value[dataPoint]);
506 jgs 126 } else if (dataPointRank==1) {
507     for (int i=0; i<dataPointShape[0]; i++) {
508 jfenwick 1796 vec[offset+i]=extract<double>(value[dataPoint][i]);
509 jgs 126 }
510     } else if (dataPointRank==2) {
511     for (int i=0; i<dataPointShape[0]; i++) {
512     for (int j=0; j<dataPointShape[1]; j++) {
513 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
514 jgs 126 }
515     }
516     } else if (dataPointRank==3) {
517     for (int i=0; i<dataPointShape[0]; i++) {
518     for (int j=0; j<dataPointShape[1]; j++) {
519     for (int k=0; k<dataPointShape[2]; k++) {
520 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
521 jgs 126 }
522     }
523     }
524     } else if (dataPointRank==4) {
525     for (int i=0; i<dataPointShape[0]; i++) {
526     for (int j=0; j<dataPointShape[1]; j++) {
527     for (int k=0; k<dataPointShape[2]; k++) {
528     for (int l=0; l<dataPointShape[3]; l++) {
529 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
530 jgs 126 }
531     }
532     }
533     }
534     }
535     dataPoint++;
536     }
537     }
538     }
539 gross 580 void
540 ksteube 775 DataExpanded::symmetric(DataAbstract* ev)
541     {
542     int sampleNo,dataPointNo;
543     int numSamples = getNumSamples();
544     int numDataPointsPerSample = getNumDPPSample();
545     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
546     if (temp_ev==0) {
547     throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
548     }
549 jfenwick 1796 ValueType& vec=getVector();
550     const ShapeType& shape=getShape();
551     ValueType& evVec=temp_ev->getVector();
552     const ShapeType& evShape=temp_ev->getShape();
553 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
554     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
555     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
556 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
557     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
558 ksteube 775 }
559     }
560     }
561     void
562     DataExpanded::nonsymmetric(DataAbstract* ev)
563     {
564     int sampleNo,dataPointNo;
565     int numSamples = getNumSamples();
566     int numDataPointsPerSample = getNumDPPSample();
567     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
568     if (temp_ev==0) {
569     throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
570     }
571 jfenwick 1796 ValueType& vec=getVector();
572     const ShapeType& shape=getShape();
573     ValueType& evVec=temp_ev->getVector();
574     const ShapeType& evShape=temp_ev->getShape();
575 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
576     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
577     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
578 jfenwick 1796 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
579     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
580 ksteube 775 }
581     }
582     }
583     void
584 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
585 ksteube 775 {
586     int sampleNo,dataPointNo;
587     int numSamples = getNumSamples();
588     int numDataPointsPerSample = getNumDPPSample();
589     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
590     if (temp_ev==0) {
591 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
592 ksteube 775 }
593 jfenwick 1796 ValueType& vec=getVector();
594     const ShapeType& shape=getShape();
595     ValueType& evVec=temp_ev->getVector();
596     const ShapeType& evShape=temp_ev->getShape();
597 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
598     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
599     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
600 jfenwick 1796 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
601     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
602 ksteube 775 }
603     }
604     }
605 gross 800
606 ksteube 775 void
607     DataExpanded::transpose(DataAbstract* ev, int axis_offset)
608     {
609     int sampleNo,dataPointNo;
610     int numSamples = getNumSamples();
611     int numDataPointsPerSample = getNumDPPSample();
612     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
613     if (temp_ev==0) {
614     throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
615     }
616 jfenwick 1796 ValueType& vec=getVector();
617     const ShapeType& shape=getShape();
618     ValueType& evVec=temp_ev->getVector();
619     const ShapeType& evShape=temp_ev->getShape();
620 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
621     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
622     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
623 jfenwick 1796 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
624     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
625 ksteube 775 }
626     }
627     }
628 gross 800
629 ksteube 775 void
630 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
631 gross 800 {
632     int sampleNo,dataPointNo;
633     int numSamples = getNumSamples();
634     int numDataPointsPerSample = getNumDPPSample();
635     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
636     if (temp_ev==0) {
637 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
638 gross 800 }
639 jfenwick 1796 ValueType& vec=getVector();
640     const ShapeType& shape=getShape();
641     ValueType& evVec=temp_ev->getVector();
642     const ShapeType& evShape=temp_ev->getShape();
643 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
644     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
647     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
648 gross 800 }
649     }
650     }
651     void
652 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
653     {
654 gross 584 int sampleNo,dataPointNo;
655 gross 580 int numSamples = getNumSamples();
656     int numDataPointsPerSample = getNumDPPSample();
657     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
658     if (temp_ev==0) {
659     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
660     }
661 jfenwick 1796 ValueType& vec=getVector();
662     const ShapeType& shape=getShape();
663     ValueType& evVec=temp_ev->getVector();
664     const ShapeType& evShape=temp_ev->getShape();
665 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
666 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
667     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
668 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
669     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
670 gross 580 }
671     }
672     }
673     void
674     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
675     {
676     int numSamples = getNumSamples();
677     int numDataPointsPerSample = getNumDPPSample();
678 gross 584 int sampleNo,dataPointNo;
679 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
680     if (temp_ev==0) {
681     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
682     }
683     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
684     if (temp_V==0) {
685     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
686     }
687 jfenwick 1796 ValueType& vec=getVector();
688     const ShapeType& shape=getShape();
689     ValueType& evVec=temp_ev->getVector();
690     const ShapeType& evShape=temp_ev->getShape();
691     ValueType& VVec=temp_V->getVector();
692     const ShapeType& VShape=temp_V->getShape();
693 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
694 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
695     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
696 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
697     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
698     VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
699 gross 580 }
700     }
701     }
702    
703 gross 950 void
704 gross 1118 DataExpanded::setToZero(){
705 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
706     // Why is there no memset here? Parallel issues?
707 gross 1118 int numSamples = getNumSamples();
708     int numDataPointsPerSample = getNumDPPSample();
709 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
710 gross 1118 double* p;
711     int sampleNo,dataPointNo, i;
712     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
713     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
714     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
715     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
716 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
717 gross 1118 }
718     }
719     }
720    
721 ksteube 1748 /* Append MPI rank to file name if multiple MPI processes */
722     char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
723     /* Make plenty of room for the mpi_rank number and terminating '\0' */
724     char *newFileName = (char *)malloc(strlen(fileName)+20);
725     strncpy(newFileName, fileName, strlen(fileName)+1);
726     if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
727     return(newFileName);
728     }
729 gross 1118
730     void
731 gross 950 DataExpanded::dump(const std::string fileName) const
732     {
733 gross 1023 #ifdef USE_NETCDF
734 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
735 gross 967 const NcDim* ncdims[ldims];
736     NcVar *var, *ids;
737 jfenwick 1796 int rank = getRank();
738 gross 967 int type= getFunctionSpace().getTypeCode();
739     int ndims =0;
740     long dims[ldims];
741 gross 1141 const double* d_ptr=&(m_data[0]);
742 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
743 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
744     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
745 ksteube 1800 #ifdef PASO_MPI
746     MPI_Status status;
747     #endif
748 gross 965
749 ksteube 1827 #ifdef PASO_MPI
750     /* Serialize NetCDF I/O */
751     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
752     #endif
753    
754 gross 967 // netCDF error handler
755     NcError err(NcError::verbose_nonfatal);
756     // Create the file.
757 ksteube 1748 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
758     NcFile dataFile(newFileName, NcFile::Replace);
759 gross 967 // check if writing was successful
760     if (!dataFile.is_valid())
761     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
762 gross 1141 if (!dataFile.add_att("type_id",2) )
763 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
764     if (!dataFile.add_att("rank",rank) )
765     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
766     if (!dataFile.add_att("function_space_type",type))
767     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
768     ndims=rank+2;
769     if ( rank >0 ) {
770     dims[0]=shape[0];
771     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
772 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
773 gross 967 }
774     if ( rank >1 ) {
775     dims[1]=shape[1];
776     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
777 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
778 gross 967 }
779     if ( rank >2 ) {
780     dims[2]=shape[2];
781     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
782 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
783 gross 967 }
784     if ( rank >3 ) {
785     dims[3]=shape[3];
786     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
787 ksteube 1800 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
788 gross 967 }
789     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
790     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
791     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
792     dims[rank+1]=getFunctionSpace().getNumSamples();
793     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
794     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
795    
796 jfenwick 1808 if (getFunctionSpace().getNumSamples()>0)
797     {
798    
799     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
800 gross 967 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
801 jfenwick 1808 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
802     if (! (ids->put(ids_p,dims[rank+1])) )
803 gross 967 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
804 jfenwick 1808 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
805 gross 967 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
806 jfenwick 1808 if (! (var->put(d_ptr,dims)) )
807 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
808 jfenwick 1808 }
809 ksteube 1800 #ifdef PASO_MPI
810 ksteube 1801 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
811 ksteube 1800 #endif
812 gross 1023 #else
813     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
814     #endif
815 gross 950 }
816 gross 1358
817 jfenwick 1796 void
818 gross 1358 DataExpanded::setTaggedValue(int tagKey,
819 jfenwick 1796 const DataTypes::ShapeType& pointshape,
820     const DataTypes::ValueType& value,
821     int dataOffset)
822 gross 1358 {
823     int numSamples = getNumSamples();
824     int numDataPointsPerSample = getNumDPPSample();
825     int sampleNo,dataPointNo, i;
826 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
827     double* p;
828     const double* in=&value[0+dataOffset];
829 gross 1358
830 jfenwick 1796 if (value.size() != n) {
831 gross 1358 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
832     }
833    
834     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
835     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
836     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
837     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
838     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
839 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
840 gross 1358 }
841     }
842     }
843     }
844    
845 jfenwick 1796
846 gross 1487 void
847     DataExpanded::reorderByReferenceIDs(int *reference_ids)
848     {
849     int numSamples = getNumSamples();
850 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
851 gross 1487 int sampleNo, sampleNo2,i;
852     double* p,*p2;
853     register double rtmp;
854     FunctionSpace fs=getFunctionSpace();
855 gross 1358
856 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
857     const int id_in=reference_ids[sampleNo];
858     const int id=fs.getReferenceIDOfSample(sampleNo);
859     if (id!=id_in) {
860     bool matched=false;
861     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
862     if (id == reference_ids[sampleNo2]) {
863     p=&(m_data[getPointOffset(sampleNo,0)]);
864     p2=&(m_data[getPointOffset(sampleNo2,0)]);
865 gross 1513 for (i=0; i<n ;i++) {
866 gross 1487 rtmp=p[i];
867     p[i]=p2[i];
868     p2[i]=rtmp;
869     }
870     reference_ids[sampleNo]=id;
871     reference_ids[sampleNo2]=id_in;
872     matched=true;
873     break;
874     }
875     }
876 phornby 1628 if (! matched) {
877 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
878     }
879     }
880     }
881     }
882    
883 jfenwick 1796 DataTypes::ValueType&
884     DataExpanded::getVector()
885     {
886     return m_data.getData();
887     }
888 gross 1487
889 jfenwick 1796 const DataTypes::ValueType&
890     DataExpanded::getVector() const
891     {
892     return m_data.getData();
893     }
894    
895 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26