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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 854 - (show annotations)
Thu Sep 21 05:29:42 2006 UTC (13 years ago) by gross
File size: 20069 byte(s)
Some modifications to the binary operations +,-,*/, pow. 
The code is a bit simpler now and more efficient has there is
no reseising required now. the resizing method has been removed as
it is very, very inefficient. Even serial code should be faster now.
It is now forbidden to do an inplace update of scalar data object with an object 
of rank >0 as this is very slow (and does not make much sense). 


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 DataAbstract*
158 DataExpanded::getSlice(const DataArrayView::RegionType& region) const
159 {
160 return new DataExpanded(*this,region);
161 }
162
163 void
164 DataExpanded::setSlice(const DataAbstract* value,
165 const DataArrayView::RegionType& region)
166 {
167 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
168 if (tempDataExp==0) {
169 throw DataException("Programming error - casting to DataExpanded.");
170 }
171 //
172 // get shape of slice
173 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
174 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
175 //
176 // check shape
177 if (getPointDataView().getRank()!=region.size()) {
178 throw DataException("Error - Invalid slice region.");
179 }
180 if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
181 throw DataException (value->getPointDataView().createShapeErrorMessage(
182 "Error - Couldn't copy slice due to shape mismatch.",shape));
183 }
184 //
185 // copy the data from the slice into this object
186 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
187 DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
188 int i, j;
189 #pragma omp parallel for private(i,j) schedule(static)
190 for (i=0;i<numRows;i++) {
191 for (j=0;j<numCols;j++) {
192 getPointDataView().copySliceFrom(getPointOffset(i,j),
193 tempDataExp->getPointDataView(),
194 tempDataExp->getPointOffset(i,j),
195 region_loop_range);
196 }
197 }
198 }
199
200 void
201 DataExpanded::copy(const DataArrayView& value)
202 {
203 //
204 // copy a single value to every data point in this object
205 int nRows=m_data.getNumRows();
206 int nCols=m_data.getNumCols();
207 int i,j;
208 #pragma omp parallel for private(i,j) schedule(static)
209 for (i=0;i<nRows;i++) {
210 for (j=0;j<nCols;j++) {
211 // NOTE: An exception may be thown from this call if
212 // DOASSERT is on which of course will play
213 // havoc with the omp threads. Run single threaded
214 // if using DOASSERT.
215 getPointDataView().copy(m_data.index(i,j),value);
216 }
217 }
218 }
219
220 void
221 DataExpanded::copy(const boost::python::numeric::array& value)
222 {
223 //
224 // first convert the numarray into a DataArray object
225 DataArray temp(value);
226 //
227 // check the input shape matches this shape
228 if (!getPointDataView().checkShape(temp.getView().getShape())) {
229 throw DataException(getPointDataView().createShapeErrorMessage(
230 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
231 temp.getView().getShape()));
232 }
233 //
234 // now copy over the data
235 copy(temp.getView());
236 }
237
238 void
239 DataExpanded::initialise(const DataArrayView::ShapeType& shape,
240 int noSamples,
241 int noDataPointsPerSample)
242 {
243 //
244 // resize data array to the required size
245 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
246 //
247 // create the data view of the data array
248 DataArrayView temp(m_data.getData(),shape);
249 setPointDataView(temp);
250 }
251
252 string
253 DataExpanded::toString() const
254 {
255 stringstream temp;
256 //
257 // create a temporary view as the offset will be changed
258 DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
259 for (int i=0;i<m_data.getNumRows();i++) {
260 for (int j=0;j<m_data.getNumCols();j++) {
261 tempView.setOffset(m_data.index(i,j));
262 stringstream suffix;
263 suffix << "(" << i << "," << j << ")";
264 temp << tempView.toString(suffix.str());
265 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
266 temp << endl;
267 }
268 }
269 }
270 return temp.str();
271 }
272
273 DataArrayView::ValueType::size_type
274 DataExpanded::getPointOffset(int sampleNo,
275 int dataPointNo) const
276 {
277 return m_data.index(sampleNo,dataPointNo);
278 }
279
280 DataArrayView
281 DataExpanded::getDataPoint(int sampleNo,
282 int dataPointNo)
283 {
284 DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
285 return temp;
286 }
287
288 DataArrayView::ValueType::size_type
289 DataExpanded::getLength() const
290 {
291 return m_data.size();
292 }
293
294 void
295 DataExpanded::setRefValue(int ref,
296 const DataArray& value)
297 {
298 //
299 // Get the number of samples and data-points per sample.
300 int numSamples = getNumSamples();
301 int numDPPSample = getNumDPPSample();
302
303 //
304 // Determine the sample number which corresponds to this reference number.
305 int sampleNo = -1;
306 int tempRef = -1;
307 for (int n=0; n<numSamples; n++) {
308 tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
309 if (tempRef == ref) {
310 sampleNo = n;
311 break;
312 }
313 }
314 if (sampleNo == -1) {
315 throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");
316 }
317
318 for (int n=0; n<numDPPSample; n++) {
319 //
320 // Get *each* data-point in the sample in turn.
321 DataArrayView pointView = getDataPoint(sampleNo, n);
322 //
323 // Assign the values in the DataArray to this data-point.
324 pointView.copy(value.getView());
325 }
326 }
327
328 void
329 DataExpanded::getRefValue(int ref,
330 DataArray& value)
331 {
332 //
333 // Get the number of samples and data-points per sample.
334 int numSamples = getNumSamples();
335 int numDPPSample = getNumDPPSample();
336
337 //
338 // Determine the sample number which corresponds to this reference number
339 int sampleNo = -1;
340 int tempRef = -1;
341 for (int n=0; n<numSamples; n++) {
342 tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
343 if (tempRef == ref) {
344 sampleNo = n;
345 break;
346 }
347 }
348 if (sampleNo == -1) {
349 throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");
350 }
351
352 //
353 // Get the *first* data-point associated with this sample number.
354 DataArrayView pointView = getDataPoint(sampleNo, 0);
355
356 //
357 // Load the values from this data-point into the DataArray
358 value.getView().copy(pointView);
359 }
360
361 int
362 DataExpanded::archiveData(ofstream& archiveFile,
363 const DataArrayView::ValueType::size_type noValues) const
364 {
365 return(m_data.archiveData(archiveFile, noValues));
366 }
367
368 int
369 DataExpanded::extractData(ifstream& archiveFile,
370 const DataArrayView::ValueType::size_type noValues)
371 {
372 return(m_data.extractData(archiveFile, noValues));
373 }
374
375 void
376 DataExpanded::copyAll(const boost::python::numeric::array& value) {
377 //
378 // Get the number of samples and data-points per sample.
379 int numSamples = getNumSamples();
380 int numDataPointsPerSample = getNumDPPSample();
381 int dataPointRank = getPointDataView().getRank();
382 ShapeType dataPointShape = getPointDataView().getShape();
383 //
384 // check rank:
385 if (value.getrank()!=dataPointRank+1)
386 throw DataException("Rank of numarray does not match Data object rank");
387 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
388 throw DataException("leading dimension of numarray is too small");
389 //
390 int dataPoint = 0;
391 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
392 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
393 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
394 if (dataPointRank==0) {
395 dataPointView()=extract<double>(value[dataPoint]);
396 } else if (dataPointRank==1) {
397 for (int i=0; i<dataPointShape[0]; i++) {
398 dataPointView(i)=extract<double>(value[dataPoint][i]);
399 }
400 } else if (dataPointRank==2) {
401 for (int i=0; i<dataPointShape[0]; i++) {
402 for (int j=0; j<dataPointShape[1]; j++) {
403 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
404 }
405 }
406 } else if (dataPointRank==3) {
407 for (int i=0; i<dataPointShape[0]; i++) {
408 for (int j=0; j<dataPointShape[1]; j++) {
409 for (int k=0; k<dataPointShape[2]; k++) {
410 dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
411 }
412 }
413 }
414 } else if (dataPointRank==4) {
415 for (int i=0; i<dataPointShape[0]; i++) {
416 for (int j=0; j<dataPointShape[1]; j++) {
417 for (int k=0; k<dataPointShape[2]; k++) {
418 for (int l=0; l<dataPointShape[3]; l++) {
419 dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
420 }
421 }
422 }
423 }
424 }
425 dataPoint++;
426 }
427 }
428 }
429 void
430 DataExpanded::symmetric(DataAbstract* ev)
431 {
432 int sampleNo,dataPointNo;
433 int numSamples = getNumSamples();
434 int numDataPointsPerSample = getNumDPPSample();
435 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
436 if (temp_ev==0) {
437 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
438 }
439 DataArrayView& thisView=getPointDataView();
440 DataArrayView& evView=ev->getPointDataView();
441 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
442 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
443 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
444 DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
445 evView,ev->getPointOffset(sampleNo,dataPointNo));
446 }
447 }
448 }
449 void
450 DataExpanded::nonsymmetric(DataAbstract* ev)
451 {
452 int sampleNo,dataPointNo;
453 int numSamples = getNumSamples();
454 int numDataPointsPerSample = getNumDPPSample();
455 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
456 if (temp_ev==0) {
457 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
458 }
459 DataArrayView& thisView=getPointDataView();
460 DataArrayView& evView=ev->getPointDataView();
461 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
462 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
463 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
464 DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
465 evView,ev->getPointOffset(sampleNo,dataPointNo));
466 }
467 }
468 }
469 void
470 DataExpanded::trace(DataAbstract* ev, int axis_offset)
471 {
472 int sampleNo,dataPointNo;
473 int numSamples = getNumSamples();
474 int numDataPointsPerSample = getNumDPPSample();
475 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
476 if (temp_ev==0) {
477 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
478 }
479 DataArrayView& thisView=getPointDataView();
480 DataArrayView& evView=ev->getPointDataView();
481 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
482 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
483 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
484 DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
485 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
486 }
487 }
488 }
489
490 void
491 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
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::transpose: 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::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
506 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
507 }
508 }
509 }
510
511 void
512 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
513 {
514 int sampleNo,dataPointNo;
515 int numSamples = getNumSamples();
516 int numDataPointsPerSample = getNumDPPSample();
517 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
518 if (temp_ev==0) {
519 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
520 }
521 DataArrayView& thisView=getPointDataView();
522 DataArrayView& evView=ev->getPointDataView();
523 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
524 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
525 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
526 DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
527 evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
528 }
529 }
530 }
531 void
532 DataExpanded::eigenvalues(DataAbstract* ev)
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::eigenvalues: 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::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
547 evView,ev->getPointOffset(sampleNo,dataPointNo));
548 }
549 }
550 }
551 void
552 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
553 {
554 int numSamples = getNumSamples();
555 int numDataPointsPerSample = getNumDPPSample();
556 int sampleNo,dataPointNo;
557 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
558 if (temp_ev==0) {
559 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
560 }
561 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
562 if (temp_V==0) {
563 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
564 }
565 DataArrayView& thisView=getPointDataView();
566 DataArrayView& evView=ev->getPointDataView();
567 DataArrayView& VView=V->getPointDataView();
568 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
569 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
570 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
571 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
572 evView,ev->getPointOffset(sampleNo,dataPointNo),
573 VView,V->getPointOffset(sampleNo,dataPointNo),
574 tol);
575 }
576 }
577 }
578
579 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26