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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/DataArrayView.cpp
File size: 19310 byte(s)
Initial revision

1 // $Id$
2
3 /*
4 ******************************************************************************
5 * *
6 * COPYRIGHT ACcESS 2004 - All Rights Reserved *
7 * *
8 * This software is the property of ACcESS. No part of this code *
9 * may be copied in any form or by any means without the expressed written *
10 * consent of ACcESS. Copying, use or modification of this software *
11 * by any unauthorised person is illegal unless that person has a software *
12 * license agreement with ACcESS. *
13 * *
14 ******************************************************************************
15 */
16
17 #include "escript/Data/DataArrayView.h"
18 #include "escript/Data/DataException.h"
19
20 #include <sstream>
21 #include <vector>
22 #include <iostream>
23 #include <boost/python/extract.hpp>
24 #include <boost/python/object.hpp>
25
26 using namespace std;
27 using namespace boost::python;
28
29 namespace escript {
30
31 DataArrayView::DataArrayView():
32 m_shape(),
33 m_offset(),
34 m_noValues()
35 {
36 m_data=0;
37 }
38
39 DataArrayView::DataArrayView(ValueType& data,
40 const ShapeType& viewShape,
41 int offset):
42 m_data(&data),
43 m_shape(viewShape),
44 m_offset(offset),
45 m_noValues(noValues(viewShape))
46 {
47 //
48 // check the shape and throw an exception if an illegal shape is specified
49 if (viewShape.size() > m_maxRank) {
50 stringstream temp;
51 temp << "Error- Couldn't construct DataArrayView, invalid rank."
52 << "Input data rank: " << viewShape.size()
53 << " Maximum rank allowed for DataArrayView: "
54 << m_maxRank;
55 throw DataException(temp.str());
56 }
57 if (m_data->size()<(noValues()+m_offset)) {
58 stringstream temp;
59 temp << "Error- Couldn't construct DataArrayView, insufficient storage."
60 << " Shape requires: " << noValues()+m_offset << " values."
61 << " Data values available: " << m_data->size();
62 throw DataException(temp.str());
63 }
64 }
65
66 DataArrayView::DataArrayView(const DataArrayView& other):
67 m_data(other.m_data),
68 m_offset(other.m_offset),
69 m_shape(other.m_shape),
70 m_noValues(other.m_noValues)
71 {
72 }
73
74 bool
75 DataArrayView::isEmpty() const
76 {
77 return(m_data==0);
78 }
79
80 void
81 DataArrayView::copy(const boost::python::numeric::array& value)
82 {
83 DataArrayView::ShapeType tempShape;
84 for (int i=0; i<value.getrank(); ++i) {
85 tempShape.push_back(extract<int>(value.getshape()[i]));
86 }
87 EsysAssert((!isEmpty()&&checkShape(tempShape)),
88 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
89
90 if (value.getrank()==0) {
91 (*this)()=extract<double>(value[value.getshape()]);
92 } else if (value.getrank()==1) {
93 for (ValueType::size_type i=0;i<tempShape[0];++i) {
94 (*this)(i)=extract<double>(value[i]);
95 }
96 } else if (value.getrank()==2) {
97 for (ValueType::size_type i=0;i<tempShape[0];++i) {
98 for (ValueType::size_type j=0;j<tempShape[1];++j) {
99 (*this)(i,j)=extract<double>(value[i][j]);
100 }
101 }
102 } else if (value.getrank()==3) {
103 for (ValueType::size_type i=0;i<tempShape[0];++i) {
104 for (ValueType::size_type j=0;j<tempShape[1];++j) {
105 for (ValueType::size_type k=0;k<tempShape[2];++k) {
106 (*this)(i,j,k)=extract<double>(value[i][j][k]);
107 }
108 }
109 }
110 } else if (value.getrank()==4) {
111 for (ValueType::size_type i=0;i<tempShape[0];++i) {
112 for (ValueType::size_type j=0;j<tempShape[1];++j) {
113 for (ValueType::size_type k=0;k<tempShape[2];++k) {
114 for (ValueType::size_type m=0;m<tempShape[3];++m) {
115 (*this)(i,j,k,m)=extract<double>(value[i][j][k][m]);
116 }
117 }
118 }
119 }
120 }
121 }
122
123 void
124 DataArrayView::copy(const DataArrayView& other)
125 {
126 copy(m_offset,other);
127 }
128
129 void
130 DataArrayView::copy(ValueType::size_type offset,
131 const DataArrayView& other)
132 {
133 EsysAssert((!isEmpty()&&checkShape(other.getShape())),
134 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
135 memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
136 }
137
138 void
139 DataArrayView::copy(ValueType::size_type thisOffset,
140 const DataArrayView& other,
141 ValueType::size_type otherOffset)
142 {
143 EsysAssert((!isEmpty()&&checkShape(other.getShape())),
144 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
145 memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
146 }
147
148 void
149 DataArrayView::copy(ValueType::size_type offset,
150 ValueType::value_type value)
151 {
152 //
153 // fill the entire view with the single value
154 EsysAssert(!isEmpty(),"Error - View is empty");
155 ValueType temp(noValues(),value);
156 memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
157 }
158
159 int
160 DataArrayView::getRank() const
161 {
162 return m_shape.size();
163 }
164
165 const DataArrayView::ShapeType&
166 DataArrayView::getShape() const
167 {
168 return m_shape;
169 }
170
171 int
172 DataArrayView::noValues(const ShapeType& shape)
173 {
174 ShapeType::const_iterator i;
175 //
176 // An empty shape vector means rank 0 which contains 1 value
177 int noValues=1;
178 for (i=shape.begin();i!=shape.end();++i) {
179 noValues*=(*i);
180 }
181 return noValues;
182 }
183
184 int
185 DataArrayView::noValues() const
186 {
187 return m_noValues;
188 }
189
190 bool
191 DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
192 {
193 return (m_shape==other);
194 }
195
196 string
197 DataArrayView::createShapeErrorMessage(const string& messagePrefix,
198 const DataArrayView::ShapeType& other) const
199 {
200 stringstream temp;
201 temp << messagePrefix
202 << " This shape: " << shapeToString(m_shape)
203 << " Other shape: " << shapeToString(other);
204 return temp.str();
205 }
206
207 void
208 DataArrayView::setOffset(ValueType::size_type offset)
209 {
210 EsysAssert((m_data->size()>=(noValues()+offset)), "Invalid offset");
211 m_offset=offset;
212 }
213
214 DataArrayView::ShapeType
215 DataArrayView::getResultSliceShape(const RegionType& region)
216 {
217 RegionType::const_iterator i;
218 DataArrayView::ShapeType result;
219 for (i=region.begin();i!=region.end();++i) {
220 int dimSize=(i->second-i->first);
221 if (dimSize!=0) {
222 result.push_back(dimSize);
223 }
224 }
225 return result;
226 }
227
228 void
229 DataArrayView::copySlice(const DataArrayView& other,
230 const RegionType& region)
231 {
232 //
233 // version of slice that uses the default offset
234 copySlice(m_offset,other,other.m_offset,region);
235 }
236
237 void
238 DataArrayView::copySlice(ValueType::size_type thisOffset,
239 const DataArrayView& other,
240 ValueType::size_type otherOffset,
241 const RegionType& region)
242 {
243 //
244 // Assume region is okay, Other asserts will detect range errors
245 EsysAssert((other.getRank()==region.size()),"Error - Invalid slice region.");
246 EsysAssert((!isEmpty()&&checkShape(getResultSliceShape(region))),
247 createShapeErrorMessage("Error - Couldn't copy slice due to shape mismatch.",getResultSliceShape(region)));
248 //
249 // take a local copy of the region and modify for the case
250 // of 2:2 being equivalent to 2:3. The only difference in meaning being in
251 // determining the output shape.
252 RegionType localRegion=region;
253 for (int i=0;i<localRegion.size();++i) {
254 if (localRegion[i].first==localRegion[i].second) {
255 ++(localRegion[i].second);
256 }
257 }
258 int numCopy=0;
259 switch (localRegion.size()) {
260 case 0:
261 (*m_data)[thisOffset]=(*other.m_data)[otherOffset];
262 break;
263 case 1:
264 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
265 (*m_data)[numCopy+thisOffset]=(*other.m_data)[i+otherOffset];
266 ++numCopy;
267 }
268 break;
269 case 2:
270 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
271 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
272 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j)+otherOffset];
273 ++numCopy;
274 }
275 }
276 break;
277 case 3:
278 for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
279 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
280 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
281 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k)+otherOffset];
282 ++numCopy;
283 }
284 }
285 }
286 break;
287 case 4:
288 for (int m=localRegion[2].first;m<localRegion[2].second;++m) {
289 for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
290 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
291 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
292 (*m_data)[numCopy+thisOffset]=(*other.m_data)[other.relIndex(i,j,k,m)+otherOffset];
293 ++numCopy;
294 }
295 }
296 }
297 }
298 break;
299 default:
300 stringstream mess;
301 mess << "Error - (copySlice) Invalid rank.";
302 throw DataException(mess.str());
303 }
304 }
305
306 void
307 DataArrayView::copySliceFrom(const DataArrayView& other,
308 const RegionType& region)
309 {
310 //
311 // version of slice that uses the default offset
312 copySliceFrom(m_offset,other,other.m_offset,region);
313 }
314
315 void
316 DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
317 const DataArrayView& other,
318 ValueType::size_type otherOffset,
319 const RegionType& region)
320 {
321 //
322 // Assume region is okay, Other asserts will detect range errors
323 EsysAssert((getRank()==region.size()),"Error - Invalid slice region.");
324 EsysAssert((!isEmpty()&&other.checkShape(getResultSliceShape(region))),
325 createShapeErrorMessage("Error - Couldn't copy slice due to shape mismatch.",
326 getResultSliceShape(region)));
327 //
328 // take a local copy of the region and modify for the case
329 // of 2:2 being equivalent to 2:3. The only difference in meaning being in
330 // determining the output shape.
331 RegionType localRegion=region;
332 for (int i=0;i<localRegion.size();++i) {
333 if (localRegion[i].first==localRegion[i].second) {
334 ++(localRegion[i].second);
335 }
336 }
337 int numCopy=0;
338 switch (localRegion.size()) {
339 case 0:
340 (*m_data)[thisOffset]=(*other.m_data)[otherOffset];
341 break;
342 case 1:
343 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
344 (*m_data)[i+thisOffset]=(*other.m_data)[numCopy+otherOffset];
345 ++numCopy;
346 }
347 break;
348 case 2:
349 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
350 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
351 (*m_data)[relIndex(i,j)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
352 ++numCopy;
353 }
354 }
355 break;
356 case 3:
357 for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
358 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
359 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
360 (*m_data)[relIndex(i,j,k)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
361 ++numCopy;
362 }
363 }
364 }
365 break;
366 case 4:
367 for (int m=localRegion[2].first;m<localRegion[2].second;++m) {
368 for (int k=localRegion[2].first;k<localRegion[2].second;++k) {
369 for (int j=localRegion[1].first;j<localRegion[1].second;++j) {
370 for (int i=localRegion[0].first;i<localRegion[0].second;++i) {
371 (*m_data)[relIndex(i,j,k,m)+thisOffset]=(*other.m_data)[numCopy+otherOffset];
372 ++numCopy;
373 }
374 }
375 }
376 }
377 break;
378 default:
379 stringstream mess;
380 mess << "Error - (copySliceFrom) Invalid rank: " << localRegion.size() << " for region.";
381 throw DataException(mess.str());
382 }
383 }
384
385 DataArrayView::ShapeType
386 DataArrayView::determineResultShape(const DataArrayView& left,
387 const DataArrayView& right)
388 {
389 DataArrayView::ShapeType temp;
390 for (int i=0; i<(left.getRank()-1); ++i) {
391 temp.push_back(left.getShape()[i]);
392 }
393 for (int i=1; i<right.getRank(); ++i) {
394 temp.push_back(right.getShape()[i]);
395 }
396 return temp;
397 }
398
399 string
400 DataArrayView::toString(const string& suffix) const
401 {
402 EsysAssert(!isEmpty(),"Error - View is empty");
403 stringstream temp;
404 string finalSuffix=suffix;
405 if (suffix.length() > 0) {
406 finalSuffix+=" ";
407 }
408 switch (getRank()) {
409 case 0:
410 temp << finalSuffix << (*this)();
411 break;
412 case 1:
413 for (int i=0;i<getShape()[0];++i) {
414 temp << "(" << i << ") " << finalSuffix << (*this)(i);
415 if (i!=(getShape()[0]-1)) {
416 temp << endl;
417 }
418 }
419 break;
420 case 2:
421 for (int i=0;i<getShape()[0];++i) {
422 for (int j=0;j<getShape()[1];++j) {
423 temp << "(" << i << "," << j << ") " << finalSuffix << (*this)(i,j);
424 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
425 temp << endl;
426 }
427 }
428 }
429 break;
430 case 3:
431 for (int i=0;i<getShape()[0];++i) {
432 for (int j=0;j<getShape()[1];++j) {
433 for (int k=0;k<getShape()[2];++k) {
434 temp << "("
435 << i << ","
436 << j << ","
437 << k << ") "
438 << finalSuffix << (*this)(i,j,k);
439 //
440 // don't put an endl after the last value
441 if (!(i==(getShape()[0]-1) &&
442 j==(getShape()[1]-1) &&
443 k==(getShape()[2]-1))) {
444 temp << endl;
445 }
446 }
447 }
448 }
449 break;
450 case 4:
451 // break;
452 default:
453 stringstream mess;
454 mess << "Error - (toString) Invalid rank: " << getRank();
455 throw DataException(mess.str());
456 }
457 return temp.str();
458 }
459
460 string
461 DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
462 {
463 stringstream temp;
464 temp << "(";
465 for (int i=0;i<shape.size();++i) {
466 temp << shape[i];
467 if (shape.size()-i != 1) {
468 temp << ",";
469 }
470 }
471 temp << ")";
472 return temp.str();
473 }
474
475 void
476 DataArrayView::matMult(const DataArrayView& left,
477 const DataArrayView& right,
478 DataArrayView& result)
479 {
480 if (left.getRank()==0 || right.getRank()==0) {
481 stringstream temp;
482 temp << "Error - (matMult) Invalid for rank 0 objects.";
483 throw DataException(temp.str());
484 }
485
486 if (left.getShape()[left.getRank()-1]!=right.getShape()[0]) {
487 stringstream temp;
488 temp << "Error - (matMult) Dimension: " << left.getRank()
489 << ", size: " << left.getShape()[left.getRank()-1]
490 << " of LHS and dimension: 1, size: " << right.getShape()[0]
491 << " of RHS don't match.";
492 throw DataException(temp.str());
493 }
494 int outputRank=left.getRank()+right.getRank()-2;
495
496 if (outputRank < 0) {
497 stringstream temp;
498 temp << "Error - (matMult) Left and right operands cannot be multiplied "
499 << "as they have incompatable rank.";
500 throw DataException(temp.str());
501 }
502
503 if (outputRank != result.getRank()) {
504 stringstream temp;
505 temp << "Error - (matMult) Rank of result array is: "
506 << result.getRank()
507 << " it must be: " << outputRank;
508 throw DataException(temp.str());
509 }
510
511 for (int i=0; i<(left.getRank()-1); ++i) {
512 if (left.getShape()[i]!=result.getShape()[i]) {
513 stringstream temp;
514 temp << "Error - (matMult) Dimension: " << i
515 << " of LHS and output array don't match.";
516 throw DataException(temp.str());
517 }
518 }
519
520 for (int i=1; i<right.getRank(); ++i) {
521 if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
522 stringstream temp;
523 temp << "Error - (matMult) Dimension: " << i
524 << ", size: " << right.getShape()[i]
525 << " of RHS and dimension: " << i+left.getRank()-1
526 << ", size: " << result.getShape()[i+left.getRank()-1]
527 << " of output array don't match.";
528 throw DataException(temp.str());
529 }
530 }
531
532 switch (left.getRank()) {
533 case 1:
534 switch (right.getRank()) {
535 case 1:
536 result()=0;
537 for (int i=0;i<left.getShape()[0];++i) {
538 result()+=left(i)*right(i);
539 }
540 break;
541 case 2:
542 for (int i=0;i<result.getShape()[0];++i) {
543 result(i)=0;
544 for (int j=0;j<right.getShape()[0];++j) {
545 result(i)+=left(j)*right(j,i);
546 }
547 }
548 break;
549 default:
550 stringstream temp;
551 temp << "Error - (matMult) Invalid rank. Programming error.";
552 throw DataException(temp.str());
553 break;
554 }
555 break;
556 case 2:
557 switch (right.getRank()) {
558 case 1:
559 result()=0;
560 for (int i=0;i<left.getShape()[0];++i) {
561 result(i)=0;
562 for (int j=0;j<left.getShape()[1];++j) {
563 result(i)+=left(i,j)*right(i);
564 }
565 }
566 break;
567 case 2:
568 for (int i=0;i<result.getShape()[0];++i) {
569 for (int j=0;j<result.getShape()[1];++j) {
570 result(i,j)=0;
571 for (int jR=0;jR<right.getShape()[0];++jR) {
572 result(i,j)+=left(i,jR)*right(jR,j);
573 }
574 }
575 }
576 break;
577 default:
578 stringstream temp;
579 temp << "Error - (matMult) Invalid rank. Programming error.";
580 throw DataException(temp.str());
581 break;
582 }
583 break;
584 default:
585 stringstream temp;
586 temp << "Error - (matMult) Not supported for rank: " << left.getRank();
587 throw DataException(temp.str());
588 break;
589 }
590 }
591
592 DataArrayView::RegionType
593 DataArrayView::getSliceRegion(const boost::python::object& key) const
594 {
595 int slice_len=0,i=0,rank=getRank(),out_len=0;
596 extract<tuple> key_tuple(key);
597 if (key_tuple.check()) {
598 slice_len=extract<int>(key.attr("__len__")());
599 } else {
600 slice_len=1;
601 }
602 if (slice_len>=rank) {
603 out_len=slice_len;
604 } else {
605 out_len=rank;
606 }
607 DataArrayView::RegionType out(out_len);
608 if (key_tuple.check()) {
609 for (i=0;i<slice_len;i++) out[i]=getSliceRange(getShape()[i],key[i]);
610 } else {
611 out[0]=getSliceRange(getShape()[0],key);
612 }
613 for (i=slice_len;i<rank;i++) out[i]=std::pair<int,int>(0,getShape()[i]);
614 // throw DataException("Error - number of slices is bigger then the rank.");
615 return out;
616 }
617
618 std::pair<int,int>
619 getSliceRange(const int s,
620 const boost::python::object& key)
621 {
622 int s1=s;
623 int s0=0;
624 /* if the key is an int we set start=end=key */
625 extract<int> slice_int(key);
626 if (slice_int.check()) {
627 s0=slice_int();
628 s1=s0;
629 } else {
630 extract<int> step(key.attr("step"));
631 if (step.check() && step()!=1) {
632 throw DataException("Error - Data does not support increments in slicing ");
633 } else {
634 extract<int> start(key.attr("start"));
635 if (start.check()) s0=start();
636 extract<int> stop(key.attr("stop"));
637 if (stop.check()) s1=stop();
638 }
639 }
640 return std::pair<int,int>(s0,s1);
641 }
642
643 bool operator==(const DataArrayView& left, const DataArrayView& right)
644 {
645 //
646 // The views are considered equal if they have the same shape and each value
647 // of the view is the same.
648 bool result=(left.m_shape==right.m_shape);
649 if (result) {
650 for (int i=0;i<left.noValues();++i) {
651 result=result && (left.getData(i)==right.getData(i));
652 }
653 }
654 return result;
655 }
656
657 bool operator!=(const DataArrayView& left, const DataArrayView& right)
658 {
659 return (!operator==(left,right));
660 }
661
662 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26