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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 967 - (hide annotations)
Tue Feb 13 09:40:12 2007 UTC (12 years, 6 months ago) by gross
File size: 25761 byte(s)
dump and load of expanded data via netCDF added. some test are still missing.
1 jgs 102 // $Id$
2 jgs 82 /*
3 elspeth 615 ************************************************************
4     * Copyright 2006 by ACcESS MNRF *
5     * *
6     * http://www.access.edu.au *
7     * Primary Business: Queensland, Australia *
8     * Licensed under the Open Software License version 3.0 *
9     * http://www.opensource.org/licenses/osl-3.0.php *
10     * *
11     ************************************************************
12 jgs 82 */
13    
14 jgs 480 #include "DataExpanded.h"
15 jgs 474 #include "DataException.h"
16     #include "DataConstant.h"
17     #include "DataTagged.h"
18 gross 967 #include <netcdfcpp.h>
19 jgs 82
20     #include <boost/python/extract.hpp>
21    
22     using namespace std;
23     using namespace boost::python;
24     using namespace boost;
25    
26     namespace escript {
27    
28 jgs 102 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
29     const FunctionSpace& what)
30     : DataAbstract(what)
31     {
32     DataArrayView::ShapeType tempShape;
33     //
34 jgs 117 // extract the shape of the python numarray
35     for (int i=0; i<value.getrank(); i++) {
36 jgs 102 tempShape.push_back(extract<int>(value.getshape()[i]));
37 jgs 82 }
38 jgs 117 //
39     // initialise the data array for this object
40 jgs 102 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
41     //
42     // copy the given value to every data point
43     copy(value);
44     }
45 jgs 82
46 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
47     : DataAbstract(other.getFunctionSpace()),
48     m_data(other.m_data)
49 jgs 117 {
50 jgs 102 //
51     // create the view for the data
52     DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
53     setPointDataView(temp);
54     }
55 jgs 82
56 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
57     : DataAbstract(other.getFunctionSpace())
58 jgs 117 {
59     //
60     // initialise the data array for this object
61 jgs 102 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
62     //
63 jgs 117 // DataConstant only has one value, copy this to every data point
64 jgs 102 copy(other.getPointDataView());
65     }
66 jgs 82
67 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
68     : DataAbstract(other.getFunctionSpace())
69 jgs 117 {
70     //
71 jgs 119 // initialise the data array for this object
72 jgs 102 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
73 jgs 117 //
74     // for each data point in this object, extract and copy the corresponding data
75     // value from the given DataTagged object
76 jgs 102 int i,j;
77     DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
78     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
79 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
80 jgs 117 for (i=0;i<numRows;i++) {
81     for (j=0;j<numCols;j++) {
82 jgs 102 try {
83 jgs 122 getPointDataView().copy(getPointOffset(i,j),
84     other.getPointDataView(),
85     other.getPointOffset(i,j));
86 jgs 82 }
87 jgs 102 catch (std::exception& e) {
88     cout << e.what() << endl;
89     }
90 jgs 82 }
91     }
92 jgs 102 }
93 jgs 82
94 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
95     const DataArrayView::RegionType& region)
96     : DataAbstract(other.getFunctionSpace())
97     {
98     //
99     // get the shape of the slice
100     DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
101     //
102     // initialise this Data object to the shape of the slice
103     initialise(shape,other.getNumSamples(),other.getNumDPPSample());
104     //
105 jgs 117 // copy the data
106 jgs 108 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
107 jgs 102 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
108     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
109     int i,j;
110 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
111 jgs 117 for (i=0;i<numRows;i++) {
112     for (j=0;j<numCols;j++) {
113 jgs 102 try {
114 jgs 122 getPointDataView().copySlice(getPointOffset(i,j),
115     other.getPointDataView(),
116     other.getPointOffset(i,j),
117     region_loop_range);
118 jgs 82 }
119 jgs 102 catch (std::exception& e) {
120     cout << e.what() << endl;
121     }
122 jgs 82 }
123     }
124 jgs 102 }
125 jgs 82
126 jgs 102 DataExpanded::DataExpanded(const DataArrayView& value,
127     const FunctionSpace& what)
128     : DataAbstract(what)
129     {
130 jgs 117 //
131     // get the shape of the given data value
132 jgs 102 DataArrayView::ShapeType tempShape=value.getShape();
133 jgs 117 //
134     // initialise this Data object to the shape of the given data value
135 jgs 102 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
136 jgs 117 //
137     // copy the given value to every data point
138 jgs 102 copy(value);
139     }
140 jgs 82
141 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
142     const DataArrayView::ShapeType &shape,
143     const DataArrayView::ValueType &data)
144     : DataAbstract(what)
145     {
146     //
147 jgs 160 // create the view of the data
148     initialise(shape,what.getNumSamples(),what.getNumDPPSample());
149     //
150 jgs 119 // copy the data in the correct format
151     m_data.getData()=data;
152     }
153    
154 jgs 102 DataExpanded::~DataExpanded()
155     {
156     }
157 jgs 82
158 jgs 102 DataAbstract*
159     DataExpanded::getSlice(const DataArrayView::RegionType& region) const
160     {
161     return new DataExpanded(*this,region);
162     }
163    
164     void
165     DataExpanded::setSlice(const DataAbstract* value,
166     const DataArrayView::RegionType& region)
167     {
168     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
169     if (tempDataExp==0) {
170     throw DataException("Programming error - casting to DataExpanded.");
171     }
172 jgs 108 //
173 jgs 117 // get shape of slice
174 jgs 108 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
175     DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
176     //
177 jgs 117 // check shape
178 jgs 102 if (getPointDataView().getRank()!=region.size()) {
179     throw DataException("Error - Invalid slice region.");
180     }
181 woo409 757 if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
182 jgs 102 throw DataException (value->getPointDataView().createShapeErrorMessage(
183 jgs 108 "Error - Couldn't copy slice due to shape mismatch.",shape));
184 jgs 102 }
185     //
186 jgs 117 // copy the data from the slice into this object
187 jgs 102 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
188     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
189     int i, j;
190 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
191 jgs 102 for (i=0;i<numRows;i++) {
192     for (j=0;j<numCols;j++) {
193 jgs 122 getPointDataView().copySliceFrom(getPointOffset(i,j),
194     tempDataExp->getPointDataView(),
195     tempDataExp->getPointOffset(i,j),
196     region_loop_range);
197 jgs 82 }
198     }
199 jgs 102 }
200 jgs 82
201 jgs 102 void
202     DataExpanded::copy(const DataArrayView& value)
203     {
204     //
205 jgs 117 // copy a single value to every data point in this object
206 jgs 102 int nRows=m_data.getNumRows();
207     int nCols=m_data.getNumCols();
208 jgs 117 int i,j;
209 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
210 jgs 117 for (i=0;i<nRows;i++) {
211     for (j=0;j<nCols;j++) {
212 jgs 102 // NOTE: An exception may be thown from this call if
213     // DOASSERT is on which of course will play
214     // havoc with the omp threads. Run single threaded
215     // if using DOASSERT.
216     getPointDataView().copy(m_data.index(i,j),value);
217 jgs 82 }
218     }
219 jgs 102 }
220 jgs 82
221 jgs 102 void
222     DataExpanded::copy(const boost::python::numeric::array& value)
223     {
224     //
225 jgs 117 // first convert the numarray into a DataArray object
226 jgs 102 DataArray temp(value);
227     //
228 jgs 117 // check the input shape matches this shape
229 jgs 102 if (!getPointDataView().checkShape(temp.getView().getShape())) {
230     throw DataException(getPointDataView().createShapeErrorMessage(
231     "Error - (DataExpanded) Cannot copy due to shape mismatch.",
232     temp.getView().getShape()));
233 jgs 82 }
234 jgs 102 //
235 jgs 117 // now copy over the data
236 jgs 102 copy(temp.getView());
237     }
238 jgs 82
239 jgs 102 void
240     DataExpanded::initialise(const DataArrayView::ShapeType& shape,
241     int noSamples,
242     int noDataPointsPerSample)
243     {
244     //
245 jgs 117 // resize data array to the required size
246 jgs 102 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
247     //
248 jgs 117 // create the data view of the data array
249 jgs 102 DataArrayView temp(m_data.getData(),shape);
250     setPointDataView(temp);
251     }
252 jgs 82
253 jgs 102 string
254     DataExpanded::toString() const
255     {
256     stringstream temp;
257 gross 856 FunctionSpace fs=getFunctionSpace();
258 jgs 102 //
259     // create a temporary view as the offset will be changed
260     DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
261 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
262     for (int j=0;j<m_data.getNumCols();j++) {
263 jgs 102 tempView.setOffset(m_data.index(i,j));
264     stringstream suffix;
265 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
266 jgs 102 temp << tempView.toString(suffix.str());
267     if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
268     temp << endl;
269 jgs 82 }
270     }
271     }
272 jgs 102 return temp.str();
273     }
274 jgs 82
275 jgs 102 DataArrayView::ValueType::size_type
276     DataExpanded::getPointOffset(int sampleNo,
277     int dataPointNo) const
278     {
279     return m_data.index(sampleNo,dataPointNo);
280     }
281 jgs 82
282 jgs 102 DataArrayView
283     DataExpanded::getDataPoint(int sampleNo,
284     int dataPointNo)
285     {
286     DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
287     return temp;
288     }
289 jgs 82
290 jgs 102 DataArrayView::ValueType::size_type
291     DataExpanded::getLength() const
292     {
293 jgs 117 return m_data.size();
294 jgs 102 }
295 jgs 82
296 jgs 123 int
297     DataExpanded::archiveData(ofstream& archiveFile,
298     const DataArrayView::ValueType::size_type noValues) const
299     {
300     return(m_data.archiveData(archiveFile, noValues));
301     }
302    
303     int
304     DataExpanded::extractData(ifstream& archiveFile,
305     const DataArrayView::ValueType::size_type noValues)
306     {
307     return(m_data.extractData(archiveFile, noValues));
308     }
309    
310 gross 921 void
311 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
312     //
313     // Get the number of samples and data-points per sample.
314     int numSamples = getNumSamples();
315     int numDataPointsPerSample = getNumDPPSample();
316     int dataPointRank = getPointDataView().getRank();
317     ShapeType dataPointShape = getPointDataView().getShape();
318     if (numSamples*numDataPointsPerSample > 0) {
319     //TODO: global error handling
320 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
321 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
322     }
323 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
324 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
325     }
326     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
327     if (dataPointRank==0) {
328     dataPointView()=value;
329     } else if (dataPointRank==1) {
330     for (int i=0; i<dataPointShape[0]; i++) {
331     dataPointView(i)=value;
332     }
333     } else if (dataPointRank==2) {
334     for (int i=0; i<dataPointShape[0]; i++) {
335     for (int j=0; j<dataPointShape[1]; j++) {
336     dataPointView(i,j)=value;
337     }
338     }
339     } else if (dataPointRank==3) {
340     for (int i=0; i<dataPointShape[0]; i++) {
341     for (int j=0; j<dataPointShape[1]; j++) {
342     for (int k=0; k<dataPointShape[2]; k++) {
343     dataPointView(i,j,k)=value;
344     }
345     }
346     }
347     } else if (dataPointRank==4) {
348     for (int i=0; i<dataPointShape[0]; i++) {
349     for (int j=0; j<dataPointShape[1]; j++) {
350     for (int k=0; k<dataPointShape[2]; k++) {
351     for (int l=0; l<dataPointShape[3]; l++) {
352     dataPointView(i,j,k,l)=value;
353     }
354     }
355     }
356     }
357     }
358     }
359     }
360     void
361 gross 921 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
362     //
363     // Get the number of samples and data-points per sample.
364     int numSamples = getNumSamples();
365     int numDataPointsPerSample = getNumDPPSample();
366     int dataPointRank = getPointDataView().getRank();
367     ShapeType dataPointShape = getPointDataView().getShape();
368     //
369     // check rank:
370     if (value.getrank()!=dataPointRank)
371     throw DataException("Rank of numarray does not match Data object rank");
372     if (numSamples*numDataPointsPerSample > 0) {
373     //TODO: global error handling
374 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
375 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
376     }
377 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
378 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
379     }
380     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
381     if (dataPointRank==0) {
382     dataPointView()=extract<double>(value[0]);
383     } else if (dataPointRank==1) {
384     for (int i=0; i<dataPointShape[0]; i++) {
385     dataPointView(i)=extract<double>(value[i]);
386     }
387     } else if (dataPointRank==2) {
388     for (int i=0; i<dataPointShape[0]; i++) {
389     for (int j=0; j<dataPointShape[1]; j++) {
390     dataPointView(i,j)=extract<double>(value[i][j]);
391     }
392     }
393     } else if (dataPointRank==3) {
394     for (int i=0; i<dataPointShape[0]; i++) {
395     for (int j=0; j<dataPointShape[1]; j++) {
396     for (int k=0; k<dataPointShape[2]; k++) {
397     dataPointView(i,j,k)=extract<double>(value[i][j][k]);
398     }
399     }
400     }
401     } else if (dataPointRank==4) {
402     for (int i=0; i<dataPointShape[0]; i++) {
403     for (int j=0; j<dataPointShape[1]; j++) {
404     for (int k=0; k<dataPointShape[2]; k++) {
405     for (int l=0; l<dataPointShape[3]; l++) {
406     dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);
407     }
408     }
409     }
410     }
411     }
412     }
413     }
414 jgs 126 void
415     DataExpanded::copyAll(const boost::python::numeric::array& value) {
416     //
417     // Get the number of samples and data-points per sample.
418     int numSamples = getNumSamples();
419     int numDataPointsPerSample = getNumDPPSample();
420     int dataPointRank = getPointDataView().getRank();
421     ShapeType dataPointShape = getPointDataView().getShape();
422     //
423     // check rank:
424     if (value.getrank()!=dataPointRank+1)
425     throw DataException("Rank of numarray does not match Data object rank");
426     if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
427     throw DataException("leading dimension of numarray is too small");
428     //
429     int dataPoint = 0;
430     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
431     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
432     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
433     if (dataPointRank==0) {
434     dataPointView()=extract<double>(value[dataPoint]);
435     } else if (dataPointRank==1) {
436     for (int i=0; i<dataPointShape[0]; i++) {
437     dataPointView(i)=extract<double>(value[dataPoint][i]);
438     }
439     } else if (dataPointRank==2) {
440     for (int i=0; i<dataPointShape[0]; i++) {
441     for (int j=0; j<dataPointShape[1]; j++) {
442     dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
443     }
444     }
445     } else if (dataPointRank==3) {
446     for (int i=0; i<dataPointShape[0]; i++) {
447     for (int j=0; j<dataPointShape[1]; j++) {
448     for (int k=0; k<dataPointShape[2]; k++) {
449     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
450     }
451     }
452     }
453     } else if (dataPointRank==4) {
454     for (int i=0; i<dataPointShape[0]; i++) {
455     for (int j=0; j<dataPointShape[1]; j++) {
456     for (int k=0; k<dataPointShape[2]; k++) {
457     for (int l=0; l<dataPointShape[3]; l++) {
458     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
459     }
460     }
461     }
462     }
463     }
464     dataPoint++;
465     }
466     }
467     }
468 gross 580 void
469 ksteube 775 DataExpanded::symmetric(DataAbstract* ev)
470     {
471     int sampleNo,dataPointNo;
472     int numSamples = getNumSamples();
473     int numDataPointsPerSample = getNumDPPSample();
474     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
475     if (temp_ev==0) {
476     throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
477     }
478     DataArrayView& thisView=getPointDataView();
479     DataArrayView& evView=ev->getPointDataView();
480     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
481     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
482     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
483     DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
484     evView,ev->getPointOffset(sampleNo,dataPointNo));
485     }
486     }
487     }
488     void
489     DataExpanded::nonsymmetric(DataAbstract* ev)
490     {
491     int sampleNo,dataPointNo;
492     int numSamples = getNumSamples();
493     int numDataPointsPerSample = getNumDPPSample();
494     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
495     if (temp_ev==0) {
496     throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
497     }
498     DataArrayView& thisView=getPointDataView();
499     DataArrayView& evView=ev->getPointDataView();
500     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
501     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503     DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
504     evView,ev->getPointOffset(sampleNo,dataPointNo));
505     }
506     }
507     }
508     void
509 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
510 ksteube 775 {
511     int sampleNo,dataPointNo;
512     int numSamples = getNumSamples();
513     int numDataPointsPerSample = getNumDPPSample();
514     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
515     if (temp_ev==0) {
516 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
517 ksteube 775 }
518     DataArrayView& thisView=getPointDataView();
519     DataArrayView& evView=ev->getPointDataView();
520     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
521     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
522     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
523 gross 800 DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
524 ksteube 775 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
525     }
526     }
527     }
528 gross 800
529 ksteube 775 void
530     DataExpanded::transpose(DataAbstract* ev, int axis_offset)
531     {
532     int sampleNo,dataPointNo;
533     int numSamples = getNumSamples();
534     int numDataPointsPerSample = getNumDPPSample();
535     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
536     if (temp_ev==0) {
537     throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
538     }
539     DataArrayView& thisView=getPointDataView();
540     DataArrayView& evView=ev->getPointDataView();
541     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
542     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
543     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
544     DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
545     evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
546     }
547     }
548     }
549 gross 800
550 ksteube 775 void
551 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
552 gross 800 {
553     int sampleNo,dataPointNo;
554     int numSamples = getNumSamples();
555     int numDataPointsPerSample = getNumDPPSample();
556     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
557     if (temp_ev==0) {
558 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
559 gross 800 }
560     DataArrayView& thisView=getPointDataView();
561     DataArrayView& evView=ev->getPointDataView();
562     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
563     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
564     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
565 gross 804 DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
566     evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
567 gross 800 }
568     }
569     }
570     void
571 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
572     {
573 gross 584 int sampleNo,dataPointNo;
574 gross 580 int numSamples = getNumSamples();
575     int numDataPointsPerSample = getNumDPPSample();
576     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
577     if (temp_ev==0) {
578     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
579     }
580     DataArrayView& thisView=getPointDataView();
581     DataArrayView& evView=ev->getPointDataView();
582     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
583 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
584     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
585 gross 580 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
586     evView,ev->getPointOffset(sampleNo,dataPointNo));
587     }
588     }
589     }
590     void
591     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
592     {
593     int numSamples = getNumSamples();
594     int numDataPointsPerSample = getNumDPPSample();
595 gross 584 int sampleNo,dataPointNo;
596 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
597     if (temp_ev==0) {
598     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
599     }
600     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
601     if (temp_V==0) {
602     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
603     }
604     DataArrayView& thisView=getPointDataView();
605     DataArrayView& evView=ev->getPointDataView();
606     DataArrayView& VView=V->getPointDataView();
607     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
608 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
609     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
610 gross 580 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
611     evView,ev->getPointOffset(sampleNo,dataPointNo),
612     VView,V->getPointOffset(sampleNo,dataPointNo),
613     tol);
614     }
615     }
616     }
617    
618 gross 950 void
619     DataExpanded::dump(const std::string fileName) const
620     {
621 gross 967 #ifdef PASO_MPI
622     throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");
623     #endif
624     const int ldims=2+DataArrayView::maxRank;
625     const NcDim* ncdims[ldims];
626     NcVar *var, *ids;
627     int rank = getPointDataView().getRank();
628     int type= getFunctionSpace().getTypeCode();
629     int ndims =0;
630     long dims[ldims];
631     DataArrayView::ShapeType shape = getPointDataView().getShape();
632 gross 965
633 gross 967 // netCDF error handler
634     NcError err(NcError::verbose_nonfatal);
635     // Create the file.
636     NcFile dataFile(fileName.c_str(), NcFile::Replace);
637     // check if writing was successful
638     if (!dataFile.is_valid())
639     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
640     if (!dataFile.add_att("type","expanded") )
641     throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
642     if (!dataFile.add_att("rank",rank) )
643     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
644     if (!dataFile.add_att("function_space_type",type))
645     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
646     ndims=rank+2;
647     if ( rank >0 ) {
648     dims[0]=shape[0];
649     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
650     throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");
651     }
652     if ( rank >1 ) {
653     dims[1]=shape[1];
654     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
655     throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");
656     }
657     if ( rank >2 ) {
658     dims[2]=shape[2];
659     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
660     throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");
661     }
662     if ( rank >3 ) {
663     dims[3]=shape[3];
664     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
665     throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");
666     }
667     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
668     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
669     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
670     dims[rank+1]=getFunctionSpace().getNumSamples();
671     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
672     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
673    
674     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
675     throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
676     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
677     if (! (ids->put(ids_p,dims[rank+1])) )
678     throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
679    
680     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
681     throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
682     if (! (var->put(&m_data[0],dims)) )
683     throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
684 gross 950 }
685    
686    
687 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