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

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

Parent Directory Parent Directory | Revision Log Revision Log


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