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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 757 - (hide annotations)
Mon Jun 26 13:12:56 2006 UTC (13 years, 3 months ago) by woo409
File size: 16984 byte(s)
+ Merge of intelc_win32 branch (revision 741:755) with trunk. Tested on iVEC altix (run_tests and py_tests all pass)

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 jgs 82
19     #include <boost/python/extract.hpp>
20    
21     using namespace std;
22     using namespace boost::python;
23     using namespace boost;
24    
25     namespace escript {
26    
27 jgs 102 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
28     const FunctionSpace& what)
29     : DataAbstract(what)
30     {
31     DataArrayView::ShapeType tempShape;
32     //
33 jgs 117 // extract the shape of the python numarray
34     for (int i=0; i<value.getrank(); i++) {
35 jgs 102 tempShape.push_back(extract<int>(value.getshape()[i]));
36 jgs 82 }
37 jgs 117 //
38     // initialise the data array for this object
39 jgs 102 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
40     //
41     // copy the given value to every data point
42     copy(value);
43     }
44 jgs 82
45 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
46     : DataAbstract(other.getFunctionSpace()),
47     m_data(other.m_data)
48 jgs 117 {
49 jgs 102 //
50     // create the view for the data
51     DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
52     setPointDataView(temp);
53     }
54 jgs 82
55 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
56     : DataAbstract(other.getFunctionSpace())
57 jgs 117 {
58     //
59     // initialise the data array for this object
60 jgs 102 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
61     //
62 jgs 117 // DataConstant only has one value, copy this to every data point
63 jgs 102 copy(other.getPointDataView());
64     }
65 jgs 82
66 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
67     : DataAbstract(other.getFunctionSpace())
68 jgs 117 {
69     //
70 jgs 119 // initialise the data array for this object
71 jgs 102 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
72 jgs 117 //
73     // for each data point in this object, extract and copy the corresponding data
74     // value from the given DataTagged object
75 jgs 102 int i,j;
76     DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
77     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
78 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
79 jgs 117 for (i=0;i<numRows;i++) {
80     for (j=0;j<numCols;j++) {
81 jgs 102 try {
82 jgs 122 getPointDataView().copy(getPointOffset(i,j),
83     other.getPointDataView(),
84     other.getPointOffset(i,j));
85 jgs 82 }
86 jgs 102 catch (std::exception& e) {
87     cout << e.what() << endl;
88     }
89 jgs 82 }
90     }
91 jgs 102 }
92 jgs 82
93 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
94     const DataArrayView::RegionType& region)
95     : DataAbstract(other.getFunctionSpace())
96     {
97     //
98     // get the shape of the slice
99     DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
100     //
101     // initialise this Data object to the shape of the slice
102     initialise(shape,other.getNumSamples(),other.getNumDPPSample());
103     //
104 jgs 117 // copy the data
105 jgs 108 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
106 jgs 102 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
107     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
108     int i,j;
109 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
110 jgs 117 for (i=0;i<numRows;i++) {
111     for (j=0;j<numCols;j++) {
112 jgs 102 try {
113 jgs 122 getPointDataView().copySlice(getPointOffset(i,j),
114     other.getPointDataView(),
115     other.getPointOffset(i,j),
116     region_loop_range);
117 jgs 82 }
118 jgs 102 catch (std::exception& e) {
119     cout << e.what() << endl;
120     }
121 jgs 82 }
122     }
123 jgs 102 }
124 jgs 82
125 jgs 102 DataExpanded::DataExpanded(const DataArrayView& value,
126     const FunctionSpace& what)
127     : DataAbstract(what)
128     {
129 jgs 117 //
130     // get the shape of the given data value
131 jgs 102 DataArrayView::ShapeType tempShape=value.getShape();
132 jgs 117 //
133     // initialise this Data object to the shape of the given data value
134 jgs 102 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
135 jgs 117 //
136     // copy the given value to every data point
137 jgs 102 copy(value);
138     }
139 jgs 82
140 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
141     const DataArrayView::ShapeType &shape,
142     const DataArrayView::ValueType &data)
143     : DataAbstract(what)
144     {
145     //
146 jgs 160 // create the view of the data
147     initialise(shape,what.getNumSamples(),what.getNumDPPSample());
148     //
149 jgs 119 // copy the data in the correct format
150     m_data.getData()=data;
151     }
152    
153 jgs 102 DataExpanded::~DataExpanded()
154     {
155     }
156 jgs 82
157 jgs 102 void
158     DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)
159     {
160     if (getPointDataView().getRank()!=0) {
161     stringstream temp;
162     temp << "Error - Can only reshape Data with data points of rank 0. "
163     << "This Data has data points with rank: "
164     << getPointDataView().getRank();
165     throw DataException(temp.str());
166 jgs 82 }
167 jgs 117 //
168     // create the new DataBlocks2D data array, and a corresponding DataArrayView
169     DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape));
170 jgs 102 DataArrayView newView(newData.getData(),shape);
171     //
172     // Copy the original data to every value for the new shape
173     int i,j;
174     int nRows=m_data.getNumRows();
175     int nCols=m_data.getNumCols();
176 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
177 jgs 117 for (i=0;i<nRows;i++) {
178     for (j=0;j<nCols;j++) {
179 jgs 102 // NOTE: An exception may be thown from this call if
180     // DOASSERT is on which of course will play
181     // havoc with the omp threads. Run single threaded
182     // if using DOASSERT.
183 jgs 117 newView.copy(newData.index(i,j),m_data(i,j));
184 jgs 82 }
185     }
186 jgs 117 // swap the new data array for the original
187 jgs 102 m_data.Swap(newData);
188 jgs 117 // set the corresponding DataArrayView
189 jgs 102 DataArrayView temp(m_data.getData(),shape);
190     setPointDataView(temp);
191     }
192 jgs 82
193 jgs 102 DataAbstract*
194     DataExpanded::getSlice(const DataArrayView::RegionType& region) const
195     {
196     return new DataExpanded(*this,region);
197     }
198    
199     void
200     DataExpanded::setSlice(const DataAbstract* value,
201     const DataArrayView::RegionType& region)
202     {
203     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
204     if (tempDataExp==0) {
205     throw DataException("Programming error - casting to DataExpanded.");
206     }
207 jgs 108 //
208 jgs 117 // get shape of slice
209 jgs 108 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
210     DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
211     //
212 jgs 117 // check shape
213 jgs 102 if (getPointDataView().getRank()!=region.size()) {
214     throw DataException("Error - Invalid slice region.");
215     }
216 woo409 757 if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
217 jgs 102 throw DataException (value->getPointDataView().createShapeErrorMessage(
218 jgs 108 "Error - Couldn't copy slice due to shape mismatch.",shape));
219 jgs 102 }
220     //
221 jgs 117 // copy the data from the slice into this object
222 jgs 102 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
223     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
224     int i, j;
225 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
226 jgs 102 for (i=0;i<numRows;i++) {
227     for (j=0;j<numCols;j++) {
228 jgs 122 getPointDataView().copySliceFrom(getPointOffset(i,j),
229     tempDataExp->getPointDataView(),
230     tempDataExp->getPointOffset(i,j),
231     region_loop_range);
232 jgs 82 }
233     }
234 jgs 102 }
235 jgs 82
236 jgs 102 void
237     DataExpanded::copy(const DataArrayView& value)
238     {
239     //
240 jgs 117 // copy a single value to every data point in this object
241 jgs 102 int nRows=m_data.getNumRows();
242     int nCols=m_data.getNumCols();
243 jgs 117 int i,j;
244 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
245 jgs 117 for (i=0;i<nRows;i++) {
246     for (j=0;j<nCols;j++) {
247 jgs 102 // NOTE: An exception may be thown from this call if
248     // DOASSERT is on which of course will play
249     // havoc with the omp threads. Run single threaded
250     // if using DOASSERT.
251     getPointDataView().copy(m_data.index(i,j),value);
252 jgs 82 }
253     }
254 jgs 102 }
255 jgs 82
256 jgs 102 void
257     DataExpanded::copy(const boost::python::numeric::array& value)
258     {
259     //
260 jgs 117 // first convert the numarray into a DataArray object
261 jgs 102 DataArray temp(value);
262     //
263 jgs 117 // check the input shape matches this shape
264 jgs 102 if (!getPointDataView().checkShape(temp.getView().getShape())) {
265     throw DataException(getPointDataView().createShapeErrorMessage(
266     "Error - (DataExpanded) Cannot copy due to shape mismatch.",
267     temp.getView().getShape()));
268 jgs 82 }
269 jgs 102 //
270 jgs 117 // now copy over the data
271 jgs 102 copy(temp.getView());
272     }
273 jgs 82
274 jgs 102 void
275     DataExpanded::initialise(const DataArrayView::ShapeType& shape,
276     int noSamples,
277     int noDataPointsPerSample)
278     {
279     //
280 jgs 117 // resize data array to the required size
281 jgs 102 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
282     //
283 jgs 117 // create the data view of the data array
284 jgs 102 DataArrayView temp(m_data.getData(),shape);
285     setPointDataView(temp);
286     }
287 jgs 82
288 jgs 102 string
289     DataExpanded::toString() const
290     {
291     stringstream temp;
292     //
293     // create a temporary view as the offset will be changed
294     DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
295 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
296     for (int j=0;j<m_data.getNumCols();j++) {
297 jgs 102 tempView.setOffset(m_data.index(i,j));
298     stringstream suffix;
299     suffix << "(" << i << "," << j << ")";
300     temp << tempView.toString(suffix.str());
301     if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
302     temp << endl;
303 jgs 82 }
304     }
305     }
306 jgs 102 return temp.str();
307     }
308 jgs 82
309 jgs 102 DataArrayView::ValueType::size_type
310     DataExpanded::getPointOffset(int sampleNo,
311     int dataPointNo) const
312     {
313     return m_data.index(sampleNo,dataPointNo);
314     }
315 jgs 82
316 jgs 102 DataArrayView
317     DataExpanded::getDataPoint(int sampleNo,
318     int dataPointNo)
319     {
320     DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
321     return temp;
322     }
323 jgs 82
324 jgs 102 DataArrayView::ValueType::size_type
325     DataExpanded::getLength() const
326     {
327 jgs 117 return m_data.size();
328 jgs 102 }
329 jgs 82
330 jgs 110 void
331     DataExpanded::setRefValue(int ref,
332     const DataArray& value)
333     {
334     //
335 jgs 113 // Get the number of samples and data-points per sample.
336 jgs 110 int numSamples = getNumSamples();
337     int numDPPSample = getNumDPPSample();
338    
339     //
340 jgs 113 // Determine the sample number which corresponds to this reference number.
341 jgs 110 int sampleNo = -1;
342     int tempRef = -1;
343     for (int n=0; n<numSamples; n++) {
344     tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
345     if (tempRef == ref) {
346     sampleNo = n;
347     break;
348     }
349     }
350     if (sampleNo == -1) {
351     throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");
352     }
353    
354 jgs 113 for (int n=0; n<numDPPSample; n++) {
355     //
356 jgs 126 // Get *each* data-point in the sample in turn.
357 jgs 113 DataArrayView pointView = getDataPoint(sampleNo, n);
358     //
359     // Assign the values in the DataArray to this data-point.
360     pointView.copy(value.getView());
361     }
362 jgs 110 }
363    
364     void
365     DataExpanded::getRefValue(int ref,
366     DataArray& value)
367     {
368     //
369 jgs 113 // Get the number of samples and data-points per sample.
370 jgs 110 int numSamples = getNumSamples();
371     int numDPPSample = getNumDPPSample();
372    
373     //
374 jgs 113 // Determine the sample number which corresponds to this reference number
375 jgs 110 int sampleNo = -1;
376     int tempRef = -1;
377     for (int n=0; n<numSamples; n++) {
378     tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
379     if (tempRef == ref) {
380     sampleNo = n;
381     break;
382     }
383     }
384     if (sampleNo == -1) {
385     throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");
386     }
387    
388 jgs 113 //
389 jgs 126 // Get the *first* data-point associated with this sample number.
390 jgs 110 DataArrayView pointView = getDataPoint(sampleNo, 0);
391    
392     //
393     // Load the values from this data-point into the DataArray
394     value.getView().copy(pointView);
395     }
396    
397 jgs 123 int
398     DataExpanded::archiveData(ofstream& archiveFile,
399     const DataArrayView::ValueType::size_type noValues) const
400     {
401     return(m_data.archiveData(archiveFile, noValues));
402     }
403    
404     int
405     DataExpanded::extractData(ifstream& archiveFile,
406     const DataArrayView::ValueType::size_type noValues)
407     {
408     return(m_data.extractData(archiveFile, noValues));
409     }
410    
411 jgs 126 void
412     DataExpanded::copyAll(const boost::python::numeric::array& value) {
413     //
414     // Get the number of samples and data-points per sample.
415     int numSamples = getNumSamples();
416     int numDataPointsPerSample = getNumDPPSample();
417     int dataPointRank = getPointDataView().getRank();
418     ShapeType dataPointShape = getPointDataView().getShape();
419     //
420     // check rank:
421     if (value.getrank()!=dataPointRank+1)
422     throw DataException("Rank of numarray does not match Data object rank");
423     if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
424     throw DataException("leading dimension of numarray is too small");
425     //
426     int dataPoint = 0;
427     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
428     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
429     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
430     if (dataPointRank==0) {
431     dataPointView()=extract<double>(value[dataPoint]);
432     } else if (dataPointRank==1) {
433     for (int i=0; i<dataPointShape[0]; i++) {
434     dataPointView(i)=extract<double>(value[dataPoint][i]);
435     }
436     } else if (dataPointRank==2) {
437     for (int i=0; i<dataPointShape[0]; i++) {
438     for (int j=0; j<dataPointShape[1]; j++) {
439     dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
440     }
441     }
442     } else if (dataPointRank==3) {
443     for (int i=0; i<dataPointShape[0]; i++) {
444     for (int j=0; j<dataPointShape[1]; j++) {
445     for (int k=0; k<dataPointShape[2]; k++) {
446     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
447     }
448     }
449     }
450     } else if (dataPointRank==4) {
451     for (int i=0; i<dataPointShape[0]; i++) {
452     for (int j=0; j<dataPointShape[1]; j++) {
453     for (int k=0; k<dataPointShape[2]; k++) {
454     for (int l=0; l<dataPointShape[3]; l++) {
455     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
456     }
457     }
458     }
459     }
460     }
461     dataPoint++;
462     }
463     }
464     }
465 gross 580 void
466     DataExpanded::eigenvalues(DataAbstract* ev)
467     {
468 gross 584 int sampleNo,dataPointNo;
469 gross 580 int numSamples = getNumSamples();
470     int numDataPointsPerSample = getNumDPPSample();
471     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
472     if (temp_ev==0) {
473     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
474     }
475     DataArrayView& thisView=getPointDataView();
476     DataArrayView& evView=ev->getPointDataView();
477     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
478 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
479     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
480 gross 580 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
481     evView,ev->getPointOffset(sampleNo,dataPointNo));
482     }
483     }
484     }
485     void
486     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
487     {
488     int numSamples = getNumSamples();
489     int numDataPointsPerSample = getNumDPPSample();
490 gross 584 int sampleNo,dataPointNo;
491 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
492     if (temp_ev==0) {
493     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
494     }
495     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
496     if (temp_V==0) {
497     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
498     }
499     DataArrayView& thisView=getPointDataView();
500     DataArrayView& evView=ev->getPointDataView();
501     DataArrayView& VView=V->getPointDataView();
502     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
503 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
504     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
505 gross 580 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
506     evView,ev->getPointOffset(sampleNo,dataPointNo),
507     VView,V->getPointOffset(sampleNo,dataPointNo),
508     tol);
509     }
510     }
511     }
512    
513 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