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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 478 - (show annotations)
Tue Jan 31 02:21:49 2006 UTC (13 years, 8 months ago) by jgs
File size: 26645 byte(s)
rationalise #includes

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 "DataArrayView.h"
18
19 #include "DataException.h"
20
21 #include <iostream>
22 #include <sstream>
23 #include <iostream>
24
25 #include <boost/python/extract.hpp>
26 #include <boost/shared_ptr.hpp>
27
28 using namespace std;
29 using namespace boost::python;
30
31 namespace escript {
32
33 DataArrayView::DataArrayView():
34 m_shape(),
35 m_offset(0),
36 m_noValues(0)
37 {
38 m_data=0;
39 }
40
41 DataArrayView::DataArrayView(ValueType& data,
42 const ShapeType& viewShape,
43 int offset):
44 m_data(&data),
45 m_shape(viewShape),
46 m_offset(offset),
47 m_noValues(noValues(viewShape))
48 {
49 //
50 // check the shape rank and size and throw an exception if a
51 // shape incompatible with the given data array is specified
52 if (viewShape.size() > m_maxRank) {
53 stringstream temp;
54 temp << "Error- Couldn't construct DataArrayView, invalid rank."
55 << "Input data rank: " << viewShape.size() << " "
56 << "Maximum rank allowed for DataArrayView: " << m_maxRank;
57 throw DataException(temp.str());
58 }
59 if (m_data->size() < (m_offset+noValues())) {
60 stringstream temp;
61 temp << "Error- Couldn't construct DataArrayView, insufficient data values. "
62 << "Shape requires: " << noValues()+m_offset << " values."
63 << "Data values available: " << m_data->size();
64 throw DataException(temp.str());
65 }
66 }
67
68 DataArrayView::DataArrayView(const DataArrayView& other):
69 m_data(other.m_data),
70 m_offset(other.m_offset),
71 m_shape(other.m_shape),
72 m_noValues(other.m_noValues)
73 {
74 }
75
76 bool
77 DataArrayView::isEmpty() const
78 {
79 return (m_data==0);
80 }
81
82 void
83 DataArrayView::copy(const boost::python::numeric::array& value)
84 {
85 ShapeType tempShape;
86 for (int i=0; i<value.getrank(); i++) {
87 tempShape.push_back(extract<int>(value.getshape()[i]));
88 }
89
90 EsysAssert((!isEmpty()&&checkShape(tempShape)),
91 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
92
93 if (value.getrank()==0) {
94 (*this)()=extract<double>(value[value.getshape()]);
95 } else if (value.getrank()==1) {
96 for (ValueType::size_type i=0;i<tempShape[0];i++) {
97 (*this)(i)=extract<double>(value[i]);
98 }
99 } else if (value.getrank()==2) {
100 for (ValueType::size_type i=0;i<tempShape[0];i++) {
101 for (ValueType::size_type j=0;j<tempShape[1];j++) {
102 (*this)(i,j)=extract<double>(value[i][j]);
103 }
104 }
105 } else if (value.getrank()==3) {
106 for (ValueType::size_type i=0;i<tempShape[0];i++) {
107 for (ValueType::size_type j=0;j<tempShape[1];j++) {
108 for (ValueType::size_type k=0;k<tempShape[2];k++) {
109 (*this)(i,j,k)=extract<double>(value[i][j][k]);
110 }
111 }
112 }
113 } else if (value.getrank()==4) {
114 for (ValueType::size_type i=0;i<tempShape[0];i++) {
115 for (ValueType::size_type j=0;j<tempShape[1];j++) {
116 for (ValueType::size_type k=0;k<tempShape[2];k++) {
117 for (ValueType::size_type l=0;l<tempShape[3];l++) {
118 (*this)(i,j,k,l)=extract<double>(value[i][j][k][l]);
119 }
120 }
121 }
122 }
123 }
124 }
125
126 void
127 DataArrayView::copy(const DataArrayView& other)
128 {
129 copy(m_offset,other);
130 }
131
132 void
133 DataArrayView::copy(ValueType::size_type offset,
134 const DataArrayView& other)
135 {
136 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
137 "Error - Couldn't copy due to insufficient storage.");
138 EsysAssert((checkShape(other.getShape())),
139 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
140 if (checkOffset(offset)) {
141 memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
142 } else {
143 throw DataException("Error - invalid offset specified.");
144 }
145 }
146
147 void
148 DataArrayView::copy(ValueType::size_type thisOffset,
149 const DataArrayView& other,
150 ValueType::size_type otherOffset)
151 {
152 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
153 "Error - Couldn't copy due to insufficient storage.");
154 EsysAssert((checkShape(other.getShape())),
155 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
156 if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
157 memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
158 } else {
159 throw DataException("Error - invalid offset specified.");
160 }
161 }
162
163 void
164 DataArrayView::copy(ValueType::size_type offset,
165 ValueType::value_type value)
166 {
167 EsysAssert((!isEmpty()&&checkOffset(offset)),
168 "Error - Couldn't copy due to insufficient storage.");
169 if (checkOffset(offset)) {
170 vector<double> temp(noValues(),value);
171 memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
172 } else {
173 throw DataException("Error - invalid offset specified.");
174 }
175 }
176
177 int
178 DataArrayView::getRank() const
179 {
180 return m_shape.size();
181 }
182
183 const
184 DataArrayView::ShapeType&
185 DataArrayView::getShape() const
186 {
187 return m_shape;
188 }
189
190 int
191 DataArrayView::noValues(const ShapeType& shape)
192 {
193 ShapeType::const_iterator i;
194 //
195 // An empty shape vector means rank 0 which contains 1 value
196 int noValues=1;
197 for (i=shape.begin();i!=shape.end();i++) {
198 noValues*=(*i);
199 }
200 return noValues;
201 }
202
203 int
204 DataArrayView::noValues(const RegionLoopRangeType& region)
205 {
206 //
207 // An empty region vector means rank 0 which contains 1 value
208 int noValues=1;
209 for (int i=0;i<region.size();i++) {
210 noValues*=region[i].second-region[i].first;
211 }
212 return noValues;
213 }
214
215 int
216 DataArrayView::noValues() const
217 {
218 return m_noValues;
219 }
220
221 bool
222 DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
223 {
224 return (m_shape==other);
225 }
226
227 string
228 DataArrayView::createShapeErrorMessage(const string& messagePrefix,
229 const DataArrayView::ShapeType& other) const
230 {
231 stringstream temp;
232 temp << messagePrefix
233 << " This shape: " << shapeToString(m_shape)
234 << " Other shape: " << shapeToString(other);
235 return temp.str();
236 }
237
238 DataArrayView::ValueType::size_type
239 DataArrayView::getOffset() const
240 {
241 return m_offset;
242 }
243
244 void
245 DataArrayView::setOffset(ValueType::size_type offset)
246 {
247 EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
248 if (checkOffset(offset)) {
249 m_offset=offset;
250 } else {
251 throw DataException("Error - invalid offset specified.");
252 }
253 }
254
255 void
256 DataArrayView::incrOffset()
257 {
258 EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
259 if (checkOffset(m_offset+noValues())) {
260 m_offset=m_offset+noValues();
261 } else {
262 throw DataException("Error - Cannot increment offset.");
263 }
264 }
265
266 bool
267 DataArrayView::checkOffset() const
268 {
269 return checkOffset(m_offset);
270 }
271
272 bool
273 DataArrayView::checkOffset(ValueType::size_type offset) const
274 {
275 return (m_data->size() >= (offset+noValues()));
276 }
277
278 DataArrayView::ValueType&
279 DataArrayView::getData() const
280 {
281 EsysAssert(!isEmpty(),"Error - View is empty");
282 return *m_data;
283 }
284
285 DataArrayView::ValueType::reference
286 DataArrayView::getData(ValueType::size_type i) const
287 {
288 //
289 // don't do any index checking so allow one past the end of the
290 // vector providing the equivalent of end()
291 return (*m_data)[m_offset+i];
292 }
293
294 DataArrayView::ShapeType
295 DataArrayView::getResultSliceShape(const RegionType& region)
296 {
297 int dimSize;
298 RegionType::const_iterator i;
299 ShapeType result;
300 for (i=region.begin();i!=region.end();i++) {
301 dimSize=((i->second)-(i->first));
302 if (dimSize!=0) {
303 result.push_back(dimSize);
304 }
305 }
306 return result;
307 }
308
309 DataArrayView::RegionType
310 DataArrayView::getSliceRegion(const boost::python::object& key) const
311 {
312 int slice_rank, i;
313 int this_rank=getRank();
314 RegionType out(this_rank);
315 /* allow for case where key is singular eg: [1], this implies we
316 want to generate a rank-1 dimension object, as opposed to eg: [1,2]
317 which implies we want to take a rank dimensional object with one
318 dimension of size 1 */
319 extract<tuple> key_tuple(key);
320 if (key_tuple.check()) {
321 slice_rank=extract<int> (key.attr("__len__")());
322 /* ensure slice is correctly dimensioned */
323 if (slice_rank>this_rank) {
324 throw DataException("Error - rank of slices does not match rank of slicee");
325 } else {
326 /* calculate values for slice range */
327 for (i=0;i<slice_rank;i++) {
328 out[i]=getSliceRange(key[i],getShape()[i]);
329 }
330 }
331 } else {
332 slice_rank=1;
333 if (slice_rank>this_rank) {
334 throw DataException("Error - rank of slices does not match rank of slicee");
335 } else {
336 out[0]=getSliceRange(key,getShape()[0]);
337 }
338 }
339 for (i=slice_rank;i<this_rank;i++) {
340 out[i]=std::pair<int,int>(0,getShape()[i]);
341 }
342 return out;
343 }
344
345 std::pair<int,int>
346 getSliceRange(const boost::python::object& key,
347 const int shape)
348 {
349 /* default slice range is range of entire shape dimension */
350 int s0=0, s1=shape;;
351 extract<int> slice_int(key);
352 if (slice_int.check()) {
353 /* if the key is a single int set start=key and end=key */
354 /* in this case, we want to return a rank-1 dimension object from
355 this object, taken from a single index value for one of this
356 object's dimensions */
357 s0=slice_int();
358 s1=s0;
359 } else {
360 /* if key is a pair extract begin and end values */
361 extract<int> step(key.attr("step"));
362 if (step.check() && step()!=1) {
363 throw DataException("Error - Data does not support increments in slicing ");
364 } else {
365 extract<int> start(key.attr("start"));
366 if (start.check()) {
367 s0=start();
368 }
369 extract<int> stop(key.attr("stop"));
370 if (stop.check()) {
371 s1=stop();
372 }
373 }
374 }
375 if (s0 < 0)
376 throw DataException("Error - slice index out of range.");
377 if (s0 == s1 and s1 >= shape)
378 throw DataException("Error - slice index out of range.");
379 if (s0 != s1 and s1>shape)
380 throw DataException("Error - slice index out of range.");
381 if (s0 > s1)
382 throw DataException("Error - lower index must less or equal upper index.");
383 return std::pair<int,int>(s0,s1);
384 }
385
386 DataArrayView::RegionLoopRangeType
387 getSliceRegionLoopRange(const DataArrayView::RegionType& region)
388 {
389 DataArrayView::RegionLoopRangeType region_loop_range(region.size());
390 for (int i=0;i<region.size();i++) {
391 if (region[i].first==region[i].second) {
392 region_loop_range[i].first=region[i].first;
393 region_loop_range[i].second=region[i].second+1;
394 } else {
395 region_loop_range[i].first=region[i].first;
396 region_loop_range[i].second=region[i].second;
397 }
398 }
399 return region_loop_range;
400 }
401
402 void
403 DataArrayView::copySlice(const DataArrayView& other,
404 const RegionLoopRangeType& region)
405 {
406 //
407 // version of copySlice that uses the default offsets
408 copySlice(m_offset,other,other.m_offset,region);
409 }
410
411 void
412 DataArrayView::copySlice(ValueType::size_type thisOffset,
413 const DataArrayView& other,
414 ValueType::size_type otherOffset,
415 const RegionLoopRangeType& region)
416 {
417
418 //
419 // Make sure views are not empty
420
421 EsysAssert(!isEmpty(),
422 "Error - this view is empty.");
423 EsysAssert(!other.isEmpty(),
424 "Error - other view is empty.");
425
426 //
427 // Check the view to be sliced from is compatible with the region to be sliced,
428 // and that the region to be sliced is compatible with this view:
429
430 EsysAssert(checkOffset(thisOffset),
431 "Error - offset incompatible with this view.");
432 EsysAssert(otherOffset+noValues()<=other.getData().size(),
433 "Error - offset incompatible with other view.");
434
435 EsysAssert(other.getRank()==region.size(),
436 "Error - slice not same rank as view to be sliced from.");
437
438 EsysAssert(noValues()==noValues(getResultSliceShape(region)),
439 "Error - slice shape not compatible shape for this view.");
440
441 //
442 // copy the values in the specified region of the other view into this view
443
444 // the following loops cannot be parallelised due to the numCopy counter
445 int numCopy=0;
446
447 switch (region.size()) {
448 case 0:
449 /* this case should never be encountered,
450 as python will never pass us an empty region.
451 here for completeness only, allows slicing of a scalar */
452 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
453 numCopy++;
454 break;
455 case 1:
456 for (int i=region[0].first;i<region[0].second;i++) {
457 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
458 numCopy++;
459 }
460 break;
461 case 2:
462 for (int j=region[1].first;j<region[1].second;j++) {
463 for (int i=region[0].first;i<region[0].second;i++) {
464 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
465 numCopy++;
466 }
467 }
468 break;
469 case 3:
470 for (int k=region[2].first;k<region[2].second;k++) {
471 for (int j=region[1].first;j<region[1].second;j++) {
472 for (int i=region[0].first;i<region[0].second;i++) {
473 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
474 numCopy++;
475 }
476 }
477 }
478 break;
479 case 4:
480 for (int l=region[3].first;l<region[3].second;l++) {
481 for (int k=region[2].first;k<region[2].second;k++) {
482 for (int j=region[1].first;j<region[1].second;j++) {
483 for (int i=region[0].first;i<region[0].second;i++) {
484 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];
485 numCopy++;
486 }
487 }
488 }
489 }
490 break;
491 default:
492 stringstream mess;
493 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
494 throw DataException(mess.str());
495 }
496 }
497
498 void
499 DataArrayView::copySliceFrom(const DataArrayView& other,
500 const RegionLoopRangeType& region_loop_range)
501 {
502 //
503 // version of copySliceFrom that uses the default offsets
504 copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
505 }
506
507 void
508 DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
509 const DataArrayView& other,
510 ValueType::size_type otherOffset,
511 const RegionLoopRangeType& region)
512 {
513
514 //
515 // Make sure views are not empty
516
517 EsysAssert(!isEmpty(),
518 "Error - this view is empty.");
519 EsysAssert(!other.isEmpty(),
520 "Error - other view is empty.");
521
522 //
523 // Check this view is compatible with the region to be sliced,
524 // and that the region to be sliced is compatible with the other view:
525
526 EsysAssert(other.checkOffset(otherOffset),
527 "Error - offset incompatible with other view.");
528 EsysAssert(thisOffset+other.noValues()<=getData().size(),
529 "Error - offset incompatible with this view.");
530
531 EsysAssert(getRank()==region.size(),
532 "Error - slice not same rank as this view.");
533
534 EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
535 "Error - slice shape not compatible shape for other view.");
536
537 //
538 // copy the values in the other view into the specified region of this view
539
540 // allow for case where other view is a scalar
541 if (other.getRank()==0) {
542
543 // the following loops cannot be parallelised due to the numCopy counter
544 int numCopy=0;
545
546 switch (region.size()) {
547 case 0:
548 /* this case should never be encountered,
549 as python will never pass us an empty region.
550 here for completeness only, allows slicing of a scalar */
551 (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
552 numCopy++;
553 break;
554 case 1:
555 for (int i=region[0].first;i<region[0].second;i++) {
556 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
557 numCopy++;
558 }
559 break;
560 case 2:
561 for (int j=region[1].first;j<region[1].second;j++) {
562 for (int i=region[0].first;i<region[0].second;i++) {
563 (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
564 numCopy++;
565 }
566 }
567 break;
568 case 3:
569 for (int k=region[2].first;k<region[2].second;k++) {
570 for (int j=region[1].first;j<region[1].second;j++) {
571 for (int i=region[0].first;i<region[0].second;i++) {
572 (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
573 numCopy++;
574 }
575 }
576 }
577 break;
578 case 4:
579 for (int l=region[3].first;l<region[3].second;l++) {
580 for (int k=region[2].first;k<region[2].second;k++) {
581 for (int j=region[1].first;j<region[1].second;j++) {
582 for (int i=region[0].first;i<region[0].second;i++) {
583 (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
584 numCopy++;
585 }
586 }
587 }
588 }
589 break;
590 default:
591 stringstream mess;
592 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
593 throw DataException(mess.str());
594 }
595
596 } else {
597
598 // the following loops cannot be parallelised due to the numCopy counter
599 int numCopy=0;
600
601 switch (region.size()) {
602 case 0:
603 /* this case should never be encountered,
604 as python will never pass us an empty region.
605 here for completeness only, allows slicing of a scalar */
606 (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
607 numCopy++;
608 break;
609 case 1:
610 for (int i=region[0].first;i<region[0].second;i++) {
611 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset+numCopy];
612 numCopy++;
613 }
614 break;
615 case 2:
616 for (int j=region[1].first;j<region[1].second;j++) {
617 for (int i=region[0].first;i<region[0].second;i++) {
618 (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset+numCopy];
619 numCopy++;
620 }
621 }
622 break;
623 case 3:
624 for (int k=region[2].first;k<region[2].second;k++) {
625 for (int j=region[1].first;j<region[1].second;j++) {
626 for (int i=region[0].first;i<region[0].second;i++) {
627 (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset+numCopy];
628 numCopy++;
629 }
630 }
631 }
632 break;
633 case 4:
634 for (int l=region[3].first;l<region[3].second;l++) {
635 for (int k=region[2].first;k<region[2].second;k++) {
636 for (int j=region[1].first;j<region[1].second;j++) {
637 for (int i=region[0].first;i<region[0].second;i++) {
638 (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset+numCopy];
639 numCopy++;
640 }
641 }
642 }
643 }
644 break;
645 default:
646 stringstream mess;
647 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
648 throw DataException(mess.str());
649 }
650
651 }
652
653 }
654
655 string
656 DataArrayView::toString(const string& suffix) const
657 {
658 EsysAssert(!isEmpty(),"Error - View is empty.");
659 stringstream temp;
660 string finalSuffix=suffix;
661 if (suffix.length() > 0) {
662 finalSuffix+=" ";
663 }
664 switch (getRank()) {
665 case 0:
666 temp << finalSuffix << (*this)();
667 break;
668 case 1:
669 for (int i=0;i<getShape()[0];i++) {
670 temp << "(" << i << ") " << finalSuffix << (*this)(i);
671 if (i!=(getShape()[0]-1)) {
672 temp << endl;
673 }
674 }
675 break;
676 case 2:
677 for (int i=0;i<getShape()[0];i++) {
678 for (int j=0;j<getShape()[1];j++) {
679 temp << "(" << i << "," << j << ") " << finalSuffix << (*this)(i,j);
680 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
681 temp << endl;
682 }
683 }
684 }
685 break;
686 case 3:
687 for (int i=0;i<getShape()[0];i++) {
688 for (int j=0;j<getShape()[1];j++) {
689 for (int k=0;k<getShape()[2];k++) {
690 temp << "(" << i << "," << j << "," << k << ") " << finalSuffix << (*this)(i,j,k);
691 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {
692 temp << endl;
693 }
694 }
695 }
696 }
697 break;
698 case 4:
699 for (int i=0;i<getShape()[0];i++) {
700 for (int j=0;j<getShape()[1];j++) {
701 for (int k=0;k<getShape()[2];k++) {
702 for (int l=0;l<getShape()[3];l++) {
703 temp << "(" << i << "," << j << "," << k << "," << l << ") " << finalSuffix << (*this)(i,j,k,l);
704 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {
705 temp << endl;
706 }
707 }
708 }
709 }
710 }
711 break;
712 default:
713 stringstream mess;
714 mess << "Error - (toString) Invalid rank: " << getRank();
715 throw DataException(mess.str());
716 }
717 return temp.str();
718 }
719
720 string
721 DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
722 {
723 stringstream temp;
724 temp << "(";
725 for (int i=0;i<shape.size();i++) {
726 temp << shape[i];
727 if (i < shape.size()-1) {
728 temp << ",";
729 }
730 }
731 temp << ")";
732 return temp.str();
733 }
734
735 void
736 DataArrayView::matMult(const DataArrayView& left,
737 const DataArrayView& right,
738 DataArrayView& result)
739 {
740
741 if (left.getRank()==0 || right.getRank()==0) {
742 stringstream temp;
743 temp << "Error - (matMult) Invalid for rank 0 objects.";
744 throw DataException(temp.str());
745 }
746
747 if (left.getShape()[left.getRank()-1] != right.getShape()[0]) {
748 stringstream temp;
749 temp << "Error - (matMult) Dimension: " << left.getRank()
750 << ", size: " << left.getShape()[left.getRank()-1]
751 << " of LHS and dimension: 1, size: " << right.getShape()[0]
752 << " of RHS don't match.";
753 throw DataException(temp.str());
754 }
755
756 int outputRank = left.getRank()+right.getRank()-2;
757
758 if (outputRank < 0) {
759 stringstream temp;
760 temp << "Error - (matMult) LHS and RHS cannot be multiplied "
761 << "as they have incompatable rank.";
762 throw DataException(temp.str());
763 }
764
765 if (outputRank != result.getRank()) {
766 stringstream temp;
767 temp << "Error - (matMult) Rank of result array is: "
768 << result.getRank()
769 << " it must be: " << outputRank;
770 throw DataException(temp.str());
771 }
772
773 for (int i=0; i<(left.getRank()-1); i++) {
774 if (left.getShape()[i] != result.getShape()[i]) {
775 stringstream temp;
776 temp << "Error - (matMult) Dimension: " << i
777 << " of LHS and result array don't match.";
778 throw DataException(temp.str());
779 }
780 }
781
782 for (int i=1; i<right.getRank(); i++) {
783 if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
784 stringstream temp;
785 temp << "Error - (matMult) Dimension: " << i
786 << ", size: " << right.getShape()[i]
787 << " of RHS and dimension: " << i+left.getRank()-1
788 << ", size: " << result.getShape()[i+left.getRank()-1]
789 << " of result array don't match.";
790 throw DataException(temp.str());
791 }
792 }
793
794 switch (left.getRank()) {
795
796 case 1:
797 switch (right.getRank()) {
798 case 1:
799 result()=0;
800 for (int i=0;i<left.getShape()[0];i++) {
801 result()+=left(i)*right(i);
802 }
803 break;
804 case 2:
805 for (int i=0;i<result.getShape()[0];i++) {
806 result(i)=0;
807 for (int j=0;j<right.getShape()[0];j++) {
808 result(i)+=left(j)*right(j,i);
809 }
810 }
811 break;
812 default:
813 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
814 throw DataException(temp.str());
815 break;
816 }
817 break;
818
819 case 2:
820 switch (right.getRank()) {
821 case 1:
822 result()=0;
823 for (int i=0;i<left.getShape()[0];i++) {
824 result(i)=0;
825 for (int j=0;j<left.getShape()[1];j++) {
826 result(i)+=left(i,j)*right(i);
827 }
828 }
829 break;
830 case 2:
831 for (int i=0;i<result.getShape()[0];i++) {
832 for (int j=0;j<result.getShape()[1];j++) {
833 result(i,j)=0;
834 for (int jR=0;jR<right.getShape()[0];jR++) {
835 result(i,j)+=left(i,jR)*right(jR,j);
836 }
837 }
838 }
839 break;
840 default:
841 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
842 throw DataException(temp.str());
843 break;
844 }
845 break;
846
847 default:
848 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << left.getRank();
849 throw DataException(temp.str());
850 break;
851 }
852
853 }
854
855 DataArrayView::ShapeType
856 DataArrayView::determineResultShape(const DataArrayView& left,
857 const DataArrayView& right)
858 {
859 ShapeType temp;
860 for (int i=0; i<(left.getRank()-1); i++) {
861 temp.push_back(left.getShape()[i]);
862 }
863 for (int i=1; i<right.getRank(); i++) {
864 temp.push_back(right.getShape()[i]);
865 }
866 return temp;
867 }
868
869 bool operator==(const DataArrayView& left, const DataArrayView& right)
870 {
871 //
872 // The views are considered equal if they have the same shape and each value
873 // of the view is the same.
874 bool result=(left.m_shape==right.m_shape);
875 if (result) {
876 for (int i=0;i<left.noValues();i++) {
877 result=result && (left.getData(i)==right.getData(i));
878 }
879 }
880 return result;
881 }
882
883 bool operator!=(const DataArrayView& left, const DataArrayView& right)
884 {
885 return (!operator==(left,right));
886 }
887
888 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26