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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 856 - (show annotations)
Tue Sep 26 01:00:36 2006 UTC (16 years, 6 months ago) by gross
File size: 20168 byte(s)
data print shows now the element reference number
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 FunctionSpace fs=getFunctionSpace();
257 //
258 // create a temporary view as the offset will be changed
259 DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
260 for (int i=0;i<m_data.getNumRows();i++) {
261 for (int j=0;j<m_data.getNumCols();j++) {
262 tempView.setOffset(m_data.index(i,j));
263 stringstream suffix;
264 suffix << "( id: " << i << ", ref: " << fs.getReferenceNoFromSampleNo(i) << ", pnt: " << j << ")";
265 temp << tempView.toString(suffix.str());
266 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
267 temp << endl;
268 }
269 }
270 }
271 return temp.str();
272 }
273
274 DataArrayView::ValueType::size_type
275 DataExpanded::getPointOffset(int sampleNo,
276 int dataPointNo) const
277 {
278 return m_data.index(sampleNo,dataPointNo);
279 }
280
281 DataArrayView
282 DataExpanded::getDataPoint(int sampleNo,
283 int dataPointNo)
284 {
285 DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
286 return temp;
287 }
288
289 DataArrayView::ValueType::size_type
290 DataExpanded::getLength() const
291 {
292 return m_data.size();
293 }
294
295 void
296 DataExpanded::setRefValue(int ref,
297 const DataArray& value)
298 {
299 //
300 // Get the number of samples and data-points per sample.
301 int numSamples = getNumSamples();
302 int numDPPSample = getNumDPPSample();
303
304 //
305 // Determine the sample number which corresponds to this reference number.
306 int sampleNo = -1;
307 int tempRef = -1;
308 for (int n=0; n<numSamples; n++) {
309 tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
310 if (tempRef == ref) {
311 sampleNo = n;
312 break;
313 }
314 }
315 if (sampleNo == -1) {
316 throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");
317 }
318
319 for (int n=0; n<numDPPSample; n++) {
320 //
321 // Get *each* data-point in the sample in turn.
322 DataArrayView pointView = getDataPoint(sampleNo, n);
323 //
324 // Assign the values in the DataArray to this data-point.
325 pointView.copy(value.getView());
326 }
327 }
328
329 void
330 DataExpanded::getRefValue(int ref,
331 DataArray& value)
332 {
333 //
334 // Get the number of samples and data-points per sample.
335 int numSamples = getNumSamples();
336 int numDPPSample = getNumDPPSample();
337
338 //
339 // Determine the sample number which corresponds to this reference number
340 int sampleNo = -1;
341 int tempRef = -1;
342 for (int n=0; n<numSamples; n++) {
343 tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
344 if (tempRef == ref) {
345 sampleNo = n;
346 break;
347 }
348 }
349 if (sampleNo == -1) {
350 throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");
351 }
352
353 //
354 // Get the *first* data-point associated with this sample number.
355 DataArrayView pointView = getDataPoint(sampleNo, 0);
356
357 //
358 // Load the values from this data-point into the DataArray
359 value.getView().copy(pointView);
360 }
361
362 int
363 DataExpanded::archiveData(ofstream& archiveFile,
364 const DataArrayView::ValueType::size_type noValues) const
365 {
366 return(m_data.archiveData(archiveFile, noValues));
367 }
368
369 int
370 DataExpanded::extractData(ifstream& archiveFile,
371 const DataArrayView::ValueType::size_type noValues)
372 {
373 return(m_data.extractData(archiveFile, noValues));
374 }
375
376 void
377 DataExpanded::copyAll(const boost::python::numeric::array& value) {
378 //
379 // Get the number of samples and data-points per sample.
380 int numSamples = getNumSamples();
381 int numDataPointsPerSample = getNumDPPSample();
382 int dataPointRank = getPointDataView().getRank();
383 ShapeType dataPointShape = getPointDataView().getShape();
384 //
385 // check rank:
386 if (value.getrank()!=dataPointRank+1)
387 throw DataException("Rank of numarray does not match Data object rank");
388 if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
389 throw DataException("leading dimension of numarray is too small");
390 //
391 int dataPoint = 0;
392 for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
393 for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
394 DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
395 if (dataPointRank==0) {
396 dataPointView()=extract<double>(value[dataPoint]);
397 } else if (dataPointRank==1) {
398 for (int i=0; i<dataPointShape[0]; i++) {
399 dataPointView(i)=extract<double>(value[dataPoint][i]);
400 }
401 } else if (dataPointRank==2) {
402 for (int i=0; i<dataPointShape[0]; i++) {
403 for (int j=0; j<dataPointShape[1]; j++) {
404 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
405 }
406 }
407 } else if (dataPointRank==3) {
408 for (int i=0; i<dataPointShape[0]; i++) {
409 for (int j=0; j<dataPointShape[1]; j++) {
410 for (int k=0; k<dataPointShape[2]; k++) {
411 dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
412 }
413 }
414 }
415 } else if (dataPointRank==4) {
416 for (int i=0; i<dataPointShape[0]; i++) {
417 for (int j=0; j<dataPointShape[1]; j++) {
418 for (int k=0; k<dataPointShape[2]; k++) {
419 for (int l=0; l<dataPointShape[3]; l++) {
420 dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
421 }
422 }
423 }
424 }
425 }
426 dataPoint++;
427 }
428 }
429 }
430 void
431 DataExpanded::symmetric(DataAbstract* ev)
432 {
433 int sampleNo,dataPointNo;
434 int numSamples = getNumSamples();
435 int numDataPointsPerSample = getNumDPPSample();
436 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
437 if (temp_ev==0) {
438 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
439 }
440 DataArrayView& thisView=getPointDataView();
441 DataArrayView& evView=ev->getPointDataView();
442 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
443 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
444 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
445 DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
446 evView,ev->getPointOffset(sampleNo,dataPointNo));
447 }
448 }
449 }
450 void
451 DataExpanded::nonsymmetric(DataAbstract* ev)
452 {
453 int sampleNo,dataPointNo;
454 int numSamples = getNumSamples();
455 int numDataPointsPerSample = getNumDPPSample();
456 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
457 if (temp_ev==0) {
458 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
459 }
460 DataArrayView& thisView=getPointDataView();
461 DataArrayView& evView=ev->getPointDataView();
462 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
463 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
464 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
465 DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
466 evView,ev->getPointOffset(sampleNo,dataPointNo));
467 }
468 }
469 }
470 void
471 DataExpanded::trace(DataAbstract* ev, int axis_offset)
472 {
473 int sampleNo,dataPointNo;
474 int numSamples = getNumSamples();
475 int numDataPointsPerSample = getNumDPPSample();
476 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
477 if (temp_ev==0) {
478 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
479 }
480 DataArrayView& thisView=getPointDataView();
481 DataArrayView& evView=ev->getPointDataView();
482 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
483 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
484 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
485 DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
486 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
487 }
488 }
489 }
490
491 void
492 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
493 {
494 int sampleNo,dataPointNo;
495 int numSamples = getNumSamples();
496 int numDataPointsPerSample = getNumDPPSample();
497 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
498 if (temp_ev==0) {
499 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
500 }
501 DataArrayView& thisView=getPointDataView();
502 DataArrayView& evView=ev->getPointDataView();
503 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
504 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
505 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
506 DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
507 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
508 }
509 }
510 }
511
512 void
513 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
514 {
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 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
521 }
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 DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
528 evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
529 }
530 }
531 }
532 void
533 DataExpanded::eigenvalues(DataAbstract* ev)
534 {
535 int sampleNo,dataPointNo;
536 int numSamples = getNumSamples();
537 int numDataPointsPerSample = getNumDPPSample();
538 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
539 if (temp_ev==0) {
540 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
541 }
542 DataArrayView& thisView=getPointDataView();
543 DataArrayView& evView=ev->getPointDataView();
544 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
545 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
546 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
547 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
548 evView,ev->getPointOffset(sampleNo,dataPointNo));
549 }
550 }
551 }
552 void
553 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
554 {
555 int numSamples = getNumSamples();
556 int numDataPointsPerSample = getNumDPPSample();
557 int sampleNo,dataPointNo;
558 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
559 if (temp_ev==0) {
560 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
561 }
562 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
563 if (temp_V==0) {
564 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
565 }
566 DataArrayView& thisView=getPointDataView();
567 DataArrayView& evView=ev->getPointDataView();
568 DataArrayView& VView=V->getPointDataView();
569 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
570 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
571 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
572 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
573 evView,ev->getPointOffset(sampleNo,dataPointNo),
574 VView,V->getPointOffset(sampleNo,dataPointNo),
575 tol);
576 }
577 }
578 }
579
580 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26