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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 775 - (show annotations)
Mon Jul 10 04:00:08 2006 UTC (13 years, 3 months ago) by ksteube
File size: 20473 byte(s)
Modified the following python methods in escript/py_src/util.py to
call faster C++ methods:
	escript_trace
	escript_transpose
	escript_symmetric
	escript_nonsymmetric

1 // $Id$
2 /*
3 ************************************************************
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 */
13
14 #include "DataExpanded.h"
15 #include "DataException.h"
16 #include "DataConstant.h"
17 #include "DataTagged.h"
18
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 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
28 const FunctionSpace& what)
29 : DataAbstract(what)
30 {
31 DataArrayView::ShapeType tempShape;
32 //
33 // extract the shape of the python numarray
34 for (int i=0; i<value.getrank(); i++) {
35 tempShape.push_back(extract<int>(value.getshape()[i]));
36 }
37 //
38 // initialise the data array for this object
39 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
40 //
41 // copy the given value to every data point
42 copy(value);
43 }
44
45 DataExpanded::DataExpanded(const DataExpanded& other)
46 : DataAbstract(other.getFunctionSpace()),
47 m_data(other.m_data)
48 {
49 //
50 // create the view for the data
51 DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
52 setPointDataView(temp);
53 }
54
55 DataExpanded::DataExpanded(const DataConstant& other)
56 : DataAbstract(other.getFunctionSpace())
57 {
58 //
59 // initialise the data array for this object
60 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
61 //
62 // DataConstant only has one value, copy this to every data point
63 copy(other.getPointDataView());
64 }
65
66 DataExpanded::DataExpanded(const DataTagged& other)
67 : DataAbstract(other.getFunctionSpace())
68 {
69 //
70 // initialise the data array for this object
71 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
72 //
73 // for each data point in this object, extract and copy the corresponding data
74 // value from the given DataTagged object
75 int i,j;
76 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
77 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
78 #pragma omp parallel for private(i,j) schedule(static)
79 for (i=0;i<numRows;i++) {
80 for (j=0;j<numCols;j++) {
81 try {
82 getPointDataView().copy(getPointOffset(i,j),
83 other.getPointDataView(),
84 other.getPointOffset(i,j));
85 }
86 catch (std::exception& e) {
87 cout << e.what() << endl;
88 }
89 }
90 }
91 }
92
93 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 // copy the data
105 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
106 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
107 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
108 int i,j;
109 #pragma omp parallel for private(i,j) schedule(static)
110 for (i=0;i<numRows;i++) {
111 for (j=0;j<numCols;j++) {
112 try {
113 getPointDataView().copySlice(getPointOffset(i,j),
114 other.getPointDataView(),
115 other.getPointOffset(i,j),
116 region_loop_range);
117 }
118 catch (std::exception& e) {
119 cout << e.what() << endl;
120 }
121 }
122 }
123 }
124
125 DataExpanded::DataExpanded(const DataArrayView& value,
126 const FunctionSpace& what)
127 : DataAbstract(what)
128 {
129 //
130 // get the shape of the given data value
131 DataArrayView::ShapeType tempShape=value.getShape();
132 //
133 // initialise this Data object to the shape of the given data value
134 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
135 //
136 // copy the given value to every data point
137 copy(value);
138 }
139
140 DataExpanded::DataExpanded(const FunctionSpace& what,
141 const DataArrayView::ShapeType &shape,
142 const DataArrayView::ValueType &data)
143 : DataAbstract(what)
144 {
145 //
146 // create the view of the data
147 initialise(shape,what.getNumSamples(),what.getNumDPPSample());
148 //
149 // copy the data in the correct format
150 m_data.getData()=data;
151 }
152
153 DataExpanded::~DataExpanded()
154 {
155 }
156
157 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 }
167 //
168 // create the new DataBlocks2D data array, and a corresponding DataArrayView
169 DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape));
170 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 #pragma omp parallel for private(i,j) schedule(static)
177 for (i=0;i<nRows;i++) {
178 for (j=0;j<nCols;j++) {
179 // 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 newView.copy(newData.index(i,j),m_data(i,j));
184 }
185 }
186 // swap the new data array for the original
187 m_data.Swap(newData);
188 // set the corresponding DataArrayView
189 DataArrayView temp(m_data.getData(),shape);
190 setPointDataView(temp);
191 }
192
193 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 //
208 // get shape of slice
209 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
210 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
211 //
212 // check shape
213 if (getPointDataView().getRank()!=region.size()) {
214 throw DataException("Error - Invalid slice region.");
215 }
216 if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
217 throw DataException (value->getPointDataView().createShapeErrorMessage(
218 "Error - Couldn't copy slice due to shape mismatch.",shape));
219 }
220 //
221 // copy the data from the slice into this object
222 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
223 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
224 int i, j;
225 #pragma omp parallel for private(i,j) schedule(static)
226 for (i=0;i<numRows;i++) {
227 for (j=0;j<numCols;j++) {
228 getPointDataView().copySliceFrom(getPointOffset(i,j),
229 tempDataExp->getPointDataView(),
230 tempDataExp->getPointOffset(i,j),
231 region_loop_range);
232 }
233 }
234 }
235
236 void
237 DataExpanded::copy(const DataArrayView& value)
238 {
239 //
240 // copy a single value to every data point in this object
241 int nRows=m_data.getNumRows();
242 int nCols=m_data.getNumCols();
243 int i,j;
244 #pragma omp parallel for private(i,j) schedule(static)
245 for (i=0;i<nRows;i++) {
246 for (j=0;j<nCols;j++) {
247 // 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 }
253 }
254 }
255
256 void
257 DataExpanded::copy(const boost::python::numeric::array& value)
258 {
259 //
260 // first convert the numarray into a DataArray object
261 DataArray temp(value);
262 //
263 // check the input shape matches this shape
264 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 }
269 //
270 // now copy over the data
271 copy(temp.getView());
272 }
273
274 void
275 DataExpanded::initialise(const DataArrayView::ShapeType& shape,
276 int noSamples,
277 int noDataPointsPerSample)
278 {
279 //
280 // resize data array to the required size
281 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
282 //
283 // create the data view of the data array
284 DataArrayView temp(m_data.getData(),shape);
285 setPointDataView(temp);
286 }
287
288 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 for (int i=0;i<m_data.getNumRows();i++) {
296 for (int j=0;j<m_data.getNumCols();j++) {
297 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 }
304 }
305 }
306 return temp.str();
307 }
308
309 DataArrayView::ValueType::size_type
310 DataExpanded::getPointOffset(int sampleNo,
311 int dataPointNo) const
312 {
313 return m_data.index(sampleNo,dataPointNo);
314 }
315
316 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
324 DataArrayView::ValueType::size_type
325 DataExpanded::getLength() const
326 {
327 return m_data.size();
328 }
329
330 void
331 DataExpanded::setRefValue(int ref,
332 const DataArray& value)
333 {
334 //
335 // Get the number of samples and data-points per sample.
336 int numSamples = getNumSamples();
337 int numDPPSample = getNumDPPSample();
338
339 //
340 // Determine the sample number which corresponds to this reference number.
341 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 for (int n=0; n<numDPPSample; n++) {
355 //
356 // Get *each* data-point in the sample in turn.
357 DataArrayView pointView = getDataPoint(sampleNo, n);
358 //
359 // Assign the values in the DataArray to this data-point.
360 pointView.copy(value.getView());
361 }
362 }
363
364 void
365 DataExpanded::getRefValue(int ref,
366 DataArray& value)
367 {
368 //
369 // Get the number of samples and data-points per sample.
370 int numSamples = getNumSamples();
371 int numDPPSample = getNumDPPSample();
372
373 //
374 // Determine the sample number which corresponds to this reference number
375 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 //
389 // Get the *first* data-point associated with this sample number.
390 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 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 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 void
466 DataExpanded::symmetric(DataAbstract* ev)
467 {
468 int sampleNo,dataPointNo;
469 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::symmetric: 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 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
479 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
480 DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
481 evView,ev->getPointOffset(sampleNo,dataPointNo));
482 }
483 }
484 }
485 void
486 DataExpanded::nonsymmetric(DataAbstract* ev)
487 {
488 int sampleNo,dataPointNo;
489 int numSamples = getNumSamples();
490 int numDataPointsPerSample = getNumDPPSample();
491 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
492 if (temp_ev==0) {
493 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
494 }
495 DataArrayView& thisView=getPointDataView();
496 DataArrayView& evView=ev->getPointDataView();
497 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
498 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
499 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
500 DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
501 evView,ev->getPointOffset(sampleNo,dataPointNo));
502 }
503 }
504 }
505 void
506 DataExpanded::matrixtrace(DataAbstract* ev, int axis_offset)
507 {
508 int sampleNo,dataPointNo;
509 int numSamples = getNumSamples();
510 int numDataPointsPerSample = getNumDPPSample();
511 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
512 if (temp_ev==0) {
513 throw DataException("Error - DataExpanded::matrixtrace: casting to DataExpanded failed (propably a programming error).");
514 }
515 DataArrayView& thisView=getPointDataView();
516 DataArrayView& evView=ev->getPointDataView();
517 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
518 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
519 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
520 DataArrayView::matrixtrace(thisView,getPointOffset(sampleNo,dataPointNo),
521 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
522 }
523 }
524 }
525 void
526 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
527 {
528 int sampleNo,dataPointNo;
529 int numSamples = getNumSamples();
530 int numDataPointsPerSample = getNumDPPSample();
531 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
532 if (temp_ev==0) {
533 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
534 }
535 DataArrayView& thisView=getPointDataView();
536 DataArrayView& evView=ev->getPointDataView();
537 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
538 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
539 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
540 DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
541 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
542 }
543 }
544 }
545 void
546 DataExpanded::eigenvalues(DataAbstract* ev)
547 {
548 int sampleNo,dataPointNo;
549 int numSamples = getNumSamples();
550 int numDataPointsPerSample = getNumDPPSample();
551 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
552 if (temp_ev==0) {
553 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
554 }
555 DataArrayView& thisView=getPointDataView();
556 DataArrayView& evView=ev->getPointDataView();
557 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
558 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
559 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
560 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
561 evView,ev->getPointOffset(sampleNo,dataPointNo));
562 }
563 }
564 }
565 void
566 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
567 {
568 int numSamples = getNumSamples();
569 int numDataPointsPerSample = getNumDPPSample();
570 int sampleNo,dataPointNo;
571 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
572 if (temp_ev==0) {
573 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
574 }
575 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
576 if (temp_V==0) {
577 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
578 }
579 DataArrayView& thisView=getPointDataView();
580 DataArrayView& evView=ev->getPointDataView();
581 DataArrayView& VView=V->getPointDataView();
582 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
583 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
584 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
585 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
586 evView,ev->getPointOffset(sampleNo,dataPointNo),
587 VView,V->getPointOffset(sampleNo,dataPointNo),
588 tol);
589 }
590 }
591 }
592
593 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26