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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (hide annotations)
Wed Sep 17 01:45:46 2008 UTC (11 years, 1 month ago) by jfenwick
File size: 33970 byte(s)
Merged noarrayview branch onto trunk.


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 1796 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 1796 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 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
216     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
217 jgs 108 //
218 jgs 117 // check shape
219 jfenwick 1796 if (getRank()!=region.size()) {
220 jgs 102 throw DataException("Error - Invalid slice region.");
221     }
222 jfenwick 1796 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
223     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 1796 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 1796 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 1796 /* getPointDataView().copySliceFrom(getPointOffset(i,j),
239 jgs 122 tempDataExp->getPointDataView(),
240     tempDataExp->getPointOffset(i,j),
241 jfenwick 1796 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 1796
248 jgs 82 }
249     }
250 jgs 102 }
251 jgs 82
252 jgs 102 void
253 jfenwick 1796 DataExpanded::copy(const DataConstant& value)
254 jgs 102 {
255 jfenwick 1796 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 1796 //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 1796
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 1796 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 1796 // 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 1796 if (!DataTypes::checkShape(getShape(),tempShape)) {
316     throw DataException(DataTypes::createShapeErrorMessage(
317 jgs 102 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
318 jfenwick 1796 tempShape,getShape()));
319 jgs 82 }
320 jgs 102 //
321 jgs 117 // now copy over the data
322 jfenwick 1796 //copy(temp_dataView);
323     getVector().copyFromNumArray(value);
324 jgs 102 }
325 jgs 82
326 jfenwick 1796
327 jgs 102 void
328 jfenwick 1796 DataExpanded::initialise(int noSamples,
329 jgs 102 int noDataPointsPerSample)
330     {
331     //
332 jgs 117 // resize data array to the required size
333 jfenwick 1796 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
334    
335 jgs 102 //
336 jfenwick 1796 // // create the data view of the data array
337     // DataArrayView temp(m_data.getData(),shape);
338     // setPointDataView(temp);
339    
340     // // 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 1796 // DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
352    
353     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 1796 offset=getPointOffset(i,j);
357 jgs 102 stringstream suffix;
358 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
359 jfenwick 1796 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 1796 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 1796 // 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 1796 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
390    
391 matt 1319 void
392 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
393     //
394     // Get the number of samples and data-points per sample.
395     int numSamples = getNumSamples();
396     int numDataPointsPerSample = getNumDPPSample();
397 jfenwick 1796 int dataPointRank = getRank();
398     ShapeType dataPointShape = getShape();
399 gross 922 if (numSamples*numDataPointsPerSample > 0) {
400     //TODO: global error handling
401 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
402 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
403     }
404 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
405 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
406     }
407 jfenwick 1796 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
408     ValueType& vec=getVector();
409 gross 922 if (dataPointRank==0) {
410 jfenwick 1796 vec[0]=value;
411 gross 922 } else if (dataPointRank==1) {
412     for (int i=0; i<dataPointShape[0]; i++) {
413 jfenwick 1796 vec[offset+i]=value;
414 gross 922 }
415     } else if (dataPointRank==2) {
416     for (int i=0; i<dataPointShape[0]; i++) {
417     for (int j=0; j<dataPointShape[1]; j++) {
418 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
419 gross 922 }
420     }
421     } else if (dataPointRank==3) {
422     for (int i=0; i<dataPointShape[0]; i++) {
423     for (int j=0; j<dataPointShape[1]; j++) {
424     for (int k=0; k<dataPointShape[2]; k++) {
425 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
426 gross 922 }
427     }
428     }
429     } else if (dataPointRank==4) {
430     for (int i=0; i<dataPointShape[0]; i++) {
431     for (int j=0; j<dataPointShape[1]; j++) {
432     for (int k=0; k<dataPointShape[2]; k++) {
433     for (int l=0; l<dataPointShape[3]; l++) {
434 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
435 gross 922 }
436     }
437     }
438     }
439 matt 1319 }
440 gross 922 }
441     }
442 matt 1319 void
443 gross 921 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
444     //
445     // Get the number of samples and data-points per sample.
446     int numSamples = getNumSamples();
447     int numDataPointsPerSample = getNumDPPSample();
448 jfenwick 1796 int dataPointRank = getRank();
449     const ShapeType& shape = getShape();
450 gross 921 //
451     // check rank:
452     if (value.getrank()!=dataPointRank)
453     throw DataException("Rank of numarray does not match Data object rank");
454     if (numSamples*numDataPointsPerSample > 0) {
455     //TODO: global error handling
456 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
457 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
458     }
459 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
460 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
461     }
462 jfenwick 1796 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
463     ValueType& vec=getVector();
464 gross 921 if (dataPointRank==0) {
465 jfenwick 1796 vec[0]=extract<double>(value[0]);
466 gross 921 } else if (dataPointRank==1) {
467 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
468     vec[i]=extract<double>(value[i]);
469 gross 921 }
470     } else if (dataPointRank==2) {
471 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
472     for (int j=0; j<shape[1]; j++) {
473     vec[getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
474 gross 921 }
475     }
476     } else if (dataPointRank==3) {
477 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
478     for (int j=0; j<shape[1]; j++) {
479     for (int k=0; k<shape[2]; k++) {
480     vec[getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
481 gross 921 }
482     }
483     }
484     } else if (dataPointRank==4) {
485 jfenwick 1796 for (int i=0; i<shape[0]; i++) {
486     for (int j=0; j<shape[1]; j++) {
487     for (int k=0; k<shape[2]; k++) {
488     for (int l=0; l<shape[3]; l++) {
489     vec[getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
490 gross 921 }
491     }
492     }
493     }
494 matt 1319 }
495 gross 921 }
496     }
497 jgs 126 void
498     DataExpanded::copyAll(const boost::python::numeric::array& value) {
499     //
500     // Get the number of samples and data-points per sample.
501     int numSamples = getNumSamples();
502     int numDataPointsPerSample = getNumDPPSample();
503 jfenwick 1796 int dataPointRank = getRank();
504     const ShapeType& dataPointShape = getShape();
505 jgs 126 //
506     // check rank:
507     if (value.getrank()!=dataPointRank+1)
508     throw DataException("Rank of numarray does not match Data object rank");
509     if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
510     throw DataException("leading dimension of numarray is too small");
511     //
512 jfenwick 1796 ValueType& vec=getVector();
513 jgs 126 int dataPoint = 0;
514     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
515     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
516 jfenwick 1796 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
517     ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
518 jgs 126 if (dataPointRank==0) {
519 jfenwick 1796 vec[offset]=extract<double>(value[dataPoint]);
520 jgs 126 } else if (dataPointRank==1) {
521     for (int i=0; i<dataPointShape[0]; i++) {
522 jfenwick 1796 vec[offset+i]=extract<double>(value[dataPoint][i]);
523 jgs 126 }
524     } else if (dataPointRank==2) {
525     for (int i=0; i<dataPointShape[0]; i++) {
526     for (int j=0; j<dataPointShape[1]; j++) {
527 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
528     // dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
529 jgs 126 }
530     }
531     } else if (dataPointRank==3) {
532     for (int i=0; i<dataPointShape[0]; i++) {
533     for (int j=0; j<dataPointShape[1]; j++) {
534     for (int k=0; k<dataPointShape[2]; k++) {
535 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
536 jgs 126 }
537     }
538     }
539     } else if (dataPointRank==4) {
540     for (int i=0; i<dataPointShape[0]; i++) {
541     for (int j=0; j<dataPointShape[1]; j++) {
542     for (int k=0; k<dataPointShape[2]; k++) {
543     for (int l=0; l<dataPointShape[3]; l++) {
544 jfenwick 1796 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
545 jgs 126 }
546     }
547     }
548     }
549     }
550     dataPoint++;
551     }
552     }
553     }
554 gross 580 void
555 ksteube 775 DataExpanded::symmetric(DataAbstract* ev)
556     {
557     int sampleNo,dataPointNo;
558     int numSamples = getNumSamples();
559     int numDataPointsPerSample = getNumDPPSample();
560     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
561     if (temp_ev==0) {
562     throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
563     }
564 jfenwick 1796 // DataArrayView& thisView=getPointDataView();
565     // DataArrayView& evView=ev->getPointDataView();
566     ValueType& vec=getVector();
567     const ShapeType& shape=getShape();
568     ValueType& evVec=temp_ev->getVector();
569     const ShapeType& evShape=temp_ev->getShape();
570 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
571     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
572     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
573 jfenwick 1796 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
574     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
575 ksteube 775 }
576     }
577     }
578     void
579     DataExpanded::nonsymmetric(DataAbstract* ev)
580     {
581     int sampleNo,dataPointNo;
582     int numSamples = getNumSamples();
583     int numDataPointsPerSample = getNumDPPSample();
584     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
585     if (temp_ev==0) {
586     throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
587     }
588 jfenwick 1796 // DataArrayView& thisView=getPointDataView();
589     // DataArrayView& evView=ev->getPointDataView();
590     ValueType& vec=getVector();
591     const ShapeType& shape=getShape();
592     ValueType& evVec=temp_ev->getVector();
593     const ShapeType& evShape=temp_ev->getShape();
594 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
595     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
596     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
597 jfenwick 1796 // DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
598     // evView,ev->getPointOffset(sampleNo,dataPointNo));
599     DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
600     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
601 ksteube 775 }
602     }
603     }
604     void
605 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
606 ksteube 775 {
607     int sampleNo,dataPointNo;
608     int numSamples = getNumSamples();
609     int numDataPointsPerSample = getNumDPPSample();
610     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
611     if (temp_ev==0) {
612 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
613 ksteube 775 }
614 jfenwick 1796 ValueType& vec=getVector();
615     const ShapeType& shape=getShape();
616     ValueType& evVec=temp_ev->getVector();
617     const ShapeType& evShape=temp_ev->getShape();
618 ksteube 775 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
619     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
620     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
621 jfenwick 1796 /* DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
622     evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);*/
623     DataMaths::trace(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     DataExpanded::transpose(DataAbstract* ev, int axis_offset)
631     {
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     throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
638     }
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 ksteube 775 #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::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
647     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
648 ksteube 775 }
649     }
650     }
651 gross 800
652 ksteube 775 void
653 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
654 gross 800 {
655     int sampleNo,dataPointNo;
656     int numSamples = getNumSamples();
657     int numDataPointsPerSample = getNumDPPSample();
658     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
659     if (temp_ev==0) {
660 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
661 gross 800 }
662 jfenwick 1796 ValueType& vec=getVector();
663     const ShapeType& shape=getShape();
664     ValueType& evVec=temp_ev->getVector();
665     const ShapeType& evShape=temp_ev->getShape();
666 gross 800 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
667     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
668     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
669 jfenwick 1796 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
670     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
671 gross 800 }
672     }
673     }
674     void
675 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
676     {
677 gross 584 int sampleNo,dataPointNo;
678 gross 580 int numSamples = getNumSamples();
679     int numDataPointsPerSample = getNumDPPSample();
680     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
681     if (temp_ev==0) {
682     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
683     }
684 jfenwick 1796 ValueType& vec=getVector();
685     const ShapeType& shape=getShape();
686     ValueType& evVec=temp_ev->getVector();
687     const ShapeType& evShape=temp_ev->getShape();
688 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
689 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
690     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
691 jfenwick 1796 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
692     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
693 gross 580 }
694     }
695     }
696     void
697     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
698     {
699     int numSamples = getNumSamples();
700     int numDataPointsPerSample = getNumDPPSample();
701 gross 584 int sampleNo,dataPointNo;
702 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
703     if (temp_ev==0) {
704     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
705     }
706     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
707     if (temp_V==0) {
708     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
709     }
710 jfenwick 1796 ValueType& vec=getVector();
711     const ShapeType& shape=getShape();
712     ValueType& evVec=temp_ev->getVector();
713     const ShapeType& evShape=temp_ev->getShape();
714     ValueType& VVec=temp_V->getVector();
715     const ShapeType& VShape=temp_V->getShape();
716 gross 580 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
717 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
718     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
719 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
720     evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
721     VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
722 gross 580 }
723     }
724     }
725    
726 gross 950 void
727 gross 1118 DataExpanded::setToZero(){
728 jfenwick 1796 // TODO: Surely there is a more efficient way to do this????
729     // Why is there no memset here? Parallel issues?
730 gross 1118 int numSamples = getNumSamples();
731     int numDataPointsPerSample = getNumDPPSample();
732 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
733 gross 1118 double* p;
734     int sampleNo,dataPointNo, i;
735     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
736     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
737     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
738     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
739 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
740 gross 1118 }
741     }
742     }
743    
744 ksteube 1748 /* Append MPI rank to file name if multiple MPI processes */
745     char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
746     /* Make plenty of room for the mpi_rank number and terminating '\0' */
747     char *newFileName = (char *)malloc(strlen(fileName)+20);
748     strncpy(newFileName, fileName, strlen(fileName)+1);
749     if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
750     return(newFileName);
751     }
752 gross 1118
753     void
754 gross 950 DataExpanded::dump(const std::string fileName) const
755     {
756 gross 1023 #ifdef USE_NETCDF
757 jfenwick 1796 const int ldims=2+DataTypes::maxRank;
758 gross 967 const NcDim* ncdims[ldims];
759     NcVar *var, *ids;
760 jfenwick 1796 int rank = getRank();
761 gross 967 int type= getFunctionSpace().getTypeCode();
762     int ndims =0;
763     long dims[ldims];
764 gross 1141 const double* d_ptr=&(m_data[0]);
765 jfenwick 1796 const DataTypes::ShapeType& shape = getShape();
766 ksteube 1748 int mpi_iam=0, mpi_num=1;
767 gross 965
768 gross 967 // netCDF error handler
769     NcError err(NcError::verbose_nonfatal);
770     // Create the file.
771 ksteube 1748 #ifdef PASO_MPI
772     MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
773     MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
774     #endif
775     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
776     NcFile dataFile(newFileName, NcFile::Replace);
777 gross 967 // check if writing was successful
778     if (!dataFile.is_valid())
779     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
780 gross 1141 if (!dataFile.add_att("type_id",2) )
781 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
782     if (!dataFile.add_att("rank",rank) )
783     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
784     if (!dataFile.add_att("function_space_type",type))
785     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
786     ndims=rank+2;
787     if ( rank >0 ) {
788     dims[0]=shape[0];
789     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
790     throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");
791     }
792     if ( rank >1 ) {
793     dims[1]=shape[1];
794     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
795     throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");
796     }
797     if ( rank >2 ) {
798     dims[2]=shape[2];
799     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
800     throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");
801     }
802     if ( rank >3 ) {
803     dims[3]=shape[3];
804     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
805     throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");
806     }
807     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
808     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
809     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
810     dims[rank+1]=getFunctionSpace().getNumSamples();
811     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
812     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
813    
814     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
815     throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
816     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
817     if (! (ids->put(ids_p,dims[rank+1])) )
818     throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
819    
820     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
821     throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
822 gross 1141 if (! (var->put(d_ptr,dims)) )
823 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
824 gross 1023 #else
825     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
826     #endif
827 gross 950 }
828 gross 1358
829 jfenwick 1796 // void
830     // DataExpanded::setTaggedValue(int tagKey,
831     // const DataArrayView& value)
832     // {
833     // int numSamples = getNumSamples();
834     // int numDataPointsPerSample = getNumDPPSample();
835     // int sampleNo,dataPointNo, i;
836     // DataTypes::ValueType::size_type n = getNoValues();
837     // double* p,*in=&(value.getData()[0]);
838     //
839     // if (value.noValues() != n) {
840     // throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
841     // }
842     //
843     // #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
844     // for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
845     // if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
846     // for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
847     // p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
848     // for (i=0; i<n ;++i) p[i]=in[i];
849     // }
850     // }
851     // }
852     // }
853    
854     void
855 gross 1358 DataExpanded::setTaggedValue(int tagKey,
856 jfenwick 1796 const DataTypes::ShapeType& pointshape,
857     const DataTypes::ValueType& value,
858     int dataOffset)
859 gross 1358 {
860     int numSamples = getNumSamples();
861     int numDataPointsPerSample = getNumDPPSample();
862     int sampleNo,dataPointNo, i;
863 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues();
864     double* p;
865     const double* in=&value[0+dataOffset];
866 gross 1358
867 jfenwick 1796 if (value.size() != n) {
868 gross 1358 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
869     }
870    
871     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
872     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
873     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
874     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
875     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
876 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
877 gross 1358 }
878     }
879     }
880     }
881    
882 jfenwick 1796
883 gross 1487 void
884     DataExpanded::reorderByReferenceIDs(int *reference_ids)
885     {
886     int numSamples = getNumSamples();
887 jfenwick 1796 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
888 gross 1487 int sampleNo, sampleNo2,i;
889     double* p,*p2;
890     register double rtmp;
891     FunctionSpace fs=getFunctionSpace();
892 gross 1358
893 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
894     const int id_in=reference_ids[sampleNo];
895     const int id=fs.getReferenceIDOfSample(sampleNo);
896     if (id!=id_in) {
897     bool matched=false;
898     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
899     if (id == reference_ids[sampleNo2]) {
900     p=&(m_data[getPointOffset(sampleNo,0)]);
901     p2=&(m_data[getPointOffset(sampleNo2,0)]);
902 gross 1513 for (i=0; i<n ;i++) {
903 gross 1487 rtmp=p[i];
904     p[i]=p2[i];
905     p2[i]=rtmp;
906     }
907     reference_ids[sampleNo]=id;
908     reference_ids[sampleNo2]=id_in;
909     matched=true;
910     break;
911     }
912     }
913 phornby 1628 if (! matched) {
914 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
915     }
916     }
917     }
918     }
919    
920 jfenwick 1796 DataTypes::ValueType&
921     DataExpanded::getVector()
922     {
923     return m_data.getData();
924     }
925 gross 1487
926 jfenwick 1796 const DataTypes::ValueType&
927     DataExpanded::getVector() const
928     {
929     return m_data.getData();
930     }
931    
932 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