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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 474 - (show annotations)
Mon Jan 30 04:23:44 2006 UTC (13 years, 8 months ago) by jgs
File size: 26489 byte(s)
restructure escript source tree
move src/Data/* -> src
remove inc
modify #includes and cpppath settings accordingly

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26