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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1748 - (hide annotations)
Wed Sep 3 06:10:39 2008 UTC (10 years, 11 months ago) by ksteube
File size: 29630 byte(s)
MPI parallelism for Data().dump and load.  Use multiple NetCDF
files, one file per MPI process

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