/[escript]/branches/subworld2/escriptcore/src/DataExpanded.cpp
ViewVC logotype

Annotation of /branches/subworld2/escriptcore/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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