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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1799 - (hide annotations)
Wed Sep 17 06:33:18 2008 UTC (11 years, 2 months ago) by jfenwick
File size: 34049 byte(s)
Added Data::copySelf() [Note: this is exposed as copy() in python].
This method returns a pointer to a deep copy of the target.
There are c++ tests but no python tests for this yet.

All DataAbstracts now have a deepCopy() which simplifies the 
implementation of the compy methods.


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