/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Annotation of /branches/arrayview_from_1695_trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1781 - (hide annotations)
Thu Sep 11 05:03:14 2008 UTC (10 years, 8 months ago) by jfenwick
File size: 34344 byte(s)
Branch commit

Merged changes from trunk version 1695 up to and including version 1779.


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