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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 126 - (show annotations)
Fri Jul 22 03:53:08 2005 UTC (13 years, 9 months ago) by jgs
File size: 26671 byte(s)
Merge of development branch back to main trunk on 2005-07-22

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26