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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 113 - (show annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 8 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/DataArrayView.cpp
File size: 26586 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26