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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 1 month ago) by ksteube
File size: 26448 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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