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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1455 - (show annotations)
Thu Feb 28 17:19:44 2008 UTC (11 years, 5 months ago) by phornby
File size: 29337 byte(s)
Merge of branches/windows_from_1431_trunk.

Revamp of the exception system.
Fix unused vars and signed/unsigned comparisons.
defined a macro THROW(ARG) in the system_dep.h's to
deal with the expectations of declarations on different architectures.

Details in the logs of branches/windows_from_1431_trunk.

pre-merge snapshot of the trunk in tags/trunk_at_1452


1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #include "DataArrayView.h"
17 #include "DataException.h"
18
19 #include <sstream>
20
21 #include <boost/python/extract.hpp>
22
23 using namespace std;
24 using namespace boost::python;
25
26 namespace escript {
27
28 DataArrayView::DataArrayView():
29 m_shape(),
30 m_offset(0),
31 m_noValues(0)
32 {
33 m_data=0;
34 }
35
36 DataArrayView::DataArrayView(ValueType& data,
37 const ShapeType& viewShape,
38 int offset):
39 m_data(&data),
40 m_shape(viewShape),
41 m_offset(offset),
42 m_noValues(noValues(viewShape))
43 {
44 //
45 // check the shape rank and size and throw an exception if a
46 // shape incompatible with the given data array is specified
47 if (viewShape.size() > m_maxRank) {
48 stringstream temp;
49 temp << "Error- Couldn't construct DataArrayView, invalid rank."
50 << "Input data rank: " << viewShape.size() << " "
51 << "Maximum rank allowed for DataArrayView: " << m_maxRank;
52 throw DataException(temp.str());
53 }
54 if (m_data->size() < (m_offset+noValues())) {
55 stringstream temp;
56 temp << "Error- Couldn't construct DataArrayView, insufficient data values. "
57 << "Shape requires: " << noValues()+m_offset << " values."
58 << "Data values available: " << m_data->size();
59 throw DataException(temp.str());
60 }
61 }
62
63 DataArrayView::DataArrayView(const DataArrayView& other):
64 m_data(other.m_data),
65 m_offset(other.m_offset),
66 m_shape(other.m_shape),
67 m_noValues(other.m_noValues)
68 {
69 }
70
71 bool
72 DataArrayView::isEmpty() const
73 {
74 return (m_data==0);
75 }
76
77 void
78 DataArrayView::copy(const boost::python::numeric::array& value)
79 {
80 ShapeType tempShape;
81 for (int i=0; i<value.getrank(); i++) {
82 tempShape.push_back(extract<int>(value.getshape()[i]));
83 }
84
85 EsysAssert((!isEmpty()&&checkShape(tempShape)),
86 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
87
88 if (value.getrank()==0) {
89 (*this)()=extract<double>(value[value.getshape()]);
90 } else if (value.getrank()==1) {
91 for (ValueType::size_type i=0;i<tempShape[0];i++) {
92 (*this)(i)=extract<double>(value[i]);
93 }
94 } else if (value.getrank()==2) {
95 for (ValueType::size_type i=0;i<tempShape[0];i++) {
96 for (ValueType::size_type j=0;j<tempShape[1];j++) {
97 (*this)(i,j)=extract<double>(value[i][j]);
98 }
99 }
100 } else if (value.getrank()==3) {
101 for (ValueType::size_type i=0;i<tempShape[0];i++) {
102 for (ValueType::size_type j=0;j<tempShape[1];j++) {
103 for (ValueType::size_type k=0;k<tempShape[2];k++) {
104 (*this)(i,j,k)=extract<double>(value[i][j][k]);
105 }
106 }
107 }
108 } else if (value.getrank()==4) {
109 for (ValueType::size_type i=0;i<tempShape[0];i++) {
110 for (ValueType::size_type j=0;j<tempShape[1];j++) {
111 for (ValueType::size_type k=0;k<tempShape[2];k++) {
112 for (ValueType::size_type l=0;l<tempShape[3];l++) {
113 (*this)(i,j,k,l)=extract<double>(value[i][j][k][l]);
114 }
115 }
116 }
117 }
118 }
119 }
120
121 void
122 DataArrayView::copy(const DataArrayView& other)
123 {
124 copy(m_offset,other);
125 }
126
127 void
128 DataArrayView::copy(ValueType::size_type offset,
129 const DataArrayView& other)
130 {
131 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
132 "Error - Couldn't copy due to insufficient storage.");
133 EsysAssert((checkShape(other.getShape())),
134 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
135 if (checkOffset(offset)) {
136 memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
137 } else {
138 throw DataException("Error - invalid offset specified.");
139 }
140 }
141
142 void
143 DataArrayView::copy(ValueType::size_type thisOffset,
144 const DataArrayView& other,
145 ValueType::size_type otherOffset)
146 {
147 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
148 "Error - Couldn't copy due to insufficient storage.");
149 EsysAssert((checkShape(other.getShape())),
150 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
151 if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
152 memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
153 } else {
154 throw DataException("Error - invalid offset specified.");
155 }
156 }
157
158 void
159 DataArrayView::copy(ValueType::size_type offset,
160 ValueType::value_type value)
161 {
162 EsysAssert((!isEmpty()&&checkOffset(offset)),
163 "Error - Couldn't copy due to insufficient storage.");
164 if (checkOffset(offset)) {
165 vector<double> temp(noValues(),value);
166 memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
167 } else {
168 throw DataException("Error - invalid offset specified.");
169 }
170 }
171
172 int
173 DataArrayView::getRank() const
174 {
175 return m_shape.size();
176 }
177
178 const
179 DataArrayView::ShapeType&
180 DataArrayView::getShape() const
181 {
182 return m_shape;
183 }
184
185 int
186 DataArrayView::noValues(const ShapeType& shape)
187 {
188 ShapeType::const_iterator i;
189 //
190 // An empty shape vector means rank 0 which contains 1 value
191 int noValues=1;
192 for (i=shape.begin();i!=shape.end();i++) {
193 noValues*=(*i);
194 }
195 return noValues;
196 }
197
198 int
199 DataArrayView::noValues(const RegionLoopRangeType& region)
200 {
201 //
202 // An empty region vector means rank 0 which contains 1 value
203 int noValues=1;
204 unsigned int i;
205 for (i=0;i<region.size();i++) {
206 noValues*=region[i].second-region[i].first;
207 }
208 return noValues;
209 }
210
211 int
212 DataArrayView::noValues() const
213 {
214 return m_noValues;
215 }
216
217 bool
218 DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
219 {
220 return (m_shape==other);
221 }
222
223 string
224 DataArrayView::createShapeErrorMessage(const string& messagePrefix,
225 const DataArrayView::ShapeType& other) const
226 {
227 stringstream temp;
228 temp << messagePrefix
229 << " This shape: " << shapeToString(m_shape)
230 << " Other shape: " << shapeToString(other);
231 return temp.str();
232 }
233
234 DataArrayView::ValueType::size_type
235 DataArrayView::getOffset() const
236 {
237 return m_offset;
238 }
239
240 void
241 DataArrayView::setOffset(ValueType::size_type offset)
242 {
243 EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
244 if (checkOffset(offset)) {
245 m_offset=offset;
246 } else {
247 throw DataException("Error - invalid offset specified.");
248 }
249 }
250
251 void
252 DataArrayView::incrOffset()
253 {
254 EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
255 if (checkOffset(m_offset+noValues())) {
256 m_offset=m_offset+noValues();
257 } else {
258 throw DataException("Error - Cannot increment offset.");
259 }
260 }
261
262 bool
263 DataArrayView::checkOffset() const
264 {
265 return checkOffset(m_offset);
266 }
267
268 bool
269 DataArrayView::checkOffset(ValueType::size_type offset) const
270 {
271 return (m_data->size() >= (offset+noValues()));
272 }
273
274 DataArrayView::ValueType&
275 DataArrayView::getData() const
276 {
277 EsysAssert(!isEmpty(),"Error - View is empty");
278 return *m_data;
279 }
280
281 DataArrayView::ValueType::reference
282 DataArrayView::getData(ValueType::size_type i) const
283 {
284 //
285 // don't do any index checking so allow one past the end of the
286 // vector providing the equivalent of end()
287 return (*m_data)[m_offset+i];
288 }
289
290 DataArrayView::ShapeType
291 DataArrayView::getResultSliceShape(const RegionType& region)
292 {
293 int dimSize;
294 RegionType::const_iterator i;
295 ShapeType result;
296 for (i=region.begin();i!=region.end();i++) {
297 dimSize=((i->second)-(i->first));
298 if (dimSize!=0) {
299 result.push_back(dimSize);
300 }
301 }
302 return result;
303 }
304
305 DataArrayView::RegionType
306 DataArrayView::getSliceRegion(const boost::python::object& key) const
307 {
308 int slice_rank, i;
309 int this_rank=getRank();
310 RegionType out(this_rank);
311 /* allow for case where key is singular eg: [1], this implies we
312 want to generate a rank-1 dimension object, as opposed to eg: [1,2]
313 which implies we want to take a rank dimensional object with one
314 dimension of size 1 */
315 extract<tuple> key_tuple(key);
316 if (key_tuple.check()) {
317 slice_rank=extract<int> (key.attr("__len__")());
318 /* ensure slice is correctly dimensioned */
319 if (slice_rank>this_rank) {
320 throw DataException("Error - rank of slices does not match rank of slicee");
321 } else {
322 /* calculate values for slice range */
323 for (i=0;i<slice_rank;i++) {
324 out[i]=getSliceRange(key[i],getShape()[i]);
325 }
326 }
327 } else {
328 slice_rank=1;
329 if (slice_rank>this_rank) {
330 throw DataException("Error - rank of slices does not match rank of slicee");
331 } else {
332 out[0]=getSliceRange(key,getShape()[0]);
333 }
334 }
335 for (i=slice_rank;i<this_rank;i++) {
336 out[i]=std::pair<int,int>(0,getShape()[i]);
337 }
338 return out;
339 }
340
341 std::pair<int,int>
342 getSliceRange(const boost::python::object& key,
343 const int shape)
344 {
345 /* default slice range is range of entire shape dimension */
346 int s0=0, s1=shape;;
347 extract<int> slice_int(key);
348 if (slice_int.check()) {
349 /* if the key is a single int set start=key and end=key */
350 /* in this case, we want to return a rank-1 dimension object from
351 this object, taken from a single index value for one of this
352 object's dimensions */
353 s0=slice_int();
354 s1=s0;
355 } else {
356 /* if key is a pair extract begin and end values */
357 extract<int> step(key.attr("step"));
358 if (step.check() && step()!=1) {
359 throw DataException("Error - Data does not support increments in slicing ");
360 } else {
361 extract<int> start(key.attr("start"));
362 if (start.check()) {
363 s0=start();
364 }
365 extract<int> stop(key.attr("stop"));
366 if (stop.check()) {
367 s1=stop();
368 }
369 }
370 }
371 if (s0 < 0)
372 throw DataException("Error - slice index out of range.");
373 if (s0 == s1 && s1 >= shape)
374 throw DataException("Error - slice index out of range.");
375 if (s0 != s1 && s1>shape)
376 throw DataException("Error - slice index out of range.");
377 if (s0 > s1)
378 throw DataException("Error - lower index must less or equal upper index.");
379 return std::pair<int,int>(s0,s1);
380 }
381
382 DataArrayView::RegionLoopRangeType
383 getSliceRegionLoopRange(const DataArrayView::RegionType& region)
384 {
385 DataArrayView::RegionLoopRangeType region_loop_range(region.size());
386 unsigned int i;
387 for (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 // the following loops cannot be parallelised due to the numCopy counter
442 int numCopy=0;
443
444 switch (region.size()) {
445 case 0:
446 /* this case should never be encountered,
447 as python will never pass us an empty region.
448 here for completeness only, allows slicing of a scalar */
449 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
450 numCopy++;
451 break;
452 case 1:
453 for (int i=region[0].first;i<region[0].second;i++) {
454 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
455 numCopy++;
456 }
457 break;
458 case 2:
459 for (int j=region[1].first;j<region[1].second;j++) {
460 for (int i=region[0].first;i<region[0].second;i++) {
461 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
462 numCopy++;
463 }
464 }
465 break;
466 case 3:
467 for (int k=region[2].first;k<region[2].second;k++) {
468 for (int j=region[1].first;j<region[1].second;j++) {
469 for (int i=region[0].first;i<region[0].second;i++) {
470 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
471 numCopy++;
472 }
473 }
474 }
475 break;
476 case 4:
477 for (int l=region[3].first;l<region[3].second;l++) {
478 for (int k=region[2].first;k<region[2].second;k++) {
479 for (int j=region[1].first;j<region[1].second;j++) {
480 for (int i=region[0].first;i<region[0].second;i++) {
481 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];
482 numCopy++;
483 }
484 }
485 }
486 }
487 break;
488 default:
489 stringstream mess;
490 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
491 throw DataException(mess.str());
492 }
493 }
494
495 void
496 DataArrayView::copySliceFrom(const DataArrayView& other,
497 const RegionLoopRangeType& region_loop_range)
498 {
499 //
500 // version of copySliceFrom that uses the default offsets
501 copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
502 }
503
504 void
505 DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
506 const DataArrayView& other,
507 ValueType::size_type otherOffset,
508 const RegionLoopRangeType& region)
509 {
510
511 //
512 // Make sure views are not empty
513
514 EsysAssert(!isEmpty(),
515 "Error - this view is empty.");
516 EsysAssert(!other.isEmpty(),
517 "Error - other view is empty.");
518
519 //
520 // Check this view is compatible with the region to be sliced,
521 // and that the region to be sliced is compatible with the other view:
522
523 EsysAssert(other.checkOffset(otherOffset),
524 "Error - offset incompatible with other view.");
525 EsysAssert(thisOffset+other.noValues()<=getData().size(),
526 "Error - offset incompatible with this view.");
527
528 EsysAssert(getRank()==region.size(),
529 "Error - slice not same rank as this view.");
530
531 EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
532 "Error - slice shape not compatible shape for other view.");
533
534 //
535 // copy the values in the other view into the specified region of this view
536
537 // allow for case where other view is a scalar
538 if (other.getRank()==0) {
539
540 // the following loops cannot be parallelised due to the numCopy counter
541 int numCopy=0;
542
543 switch (region.size()) {
544 case 0:
545 /* this case should never be encountered,
546 as python will never pass us an empty region.
547 here for completeness only, allows slicing of a scalar */
548 (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
549 numCopy++;
550 break;
551 case 1:
552 for (int i=region[0].first;i<region[0].second;i++) {
553 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
554 numCopy++;
555 }
556 break;
557 case 2:
558 for (int j=region[1].first;j<region[1].second;j++) {
559 for (int i=region[0].first;i<region[0].second;i++) {
560 (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
561 numCopy++;
562 }
563 }
564 break;
565 case 3:
566 for (int k=region[2].first;k<region[2].second;k++) {
567 for (int j=region[1].first;j<region[1].second;j++) {
568 for (int i=region[0].first;i<region[0].second;i++) {
569 (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
570 numCopy++;
571 }
572 }
573 }
574 break;
575 case 4:
576 for (int l=region[3].first;l<region[3].second;l++) {
577 for (int k=region[2].first;k<region[2].second;k++) {
578 for (int j=region[1].first;j<region[1].second;j++) {
579 for (int i=region[0].first;i<region[0].second;i++) {
580 (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
581 numCopy++;
582 }
583 }
584 }
585 }
586 break;
587 default:
588 stringstream mess;
589 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
590 throw DataException(mess.str());
591 }
592
593 } else {
594
595 // the following loops cannot be parallelised due to the numCopy counter
596 int numCopy=0;
597
598 switch (region.size()) {
599 case 0:
600 /* this case should never be encountered,
601 as python will never pass us an empty region.
602 here for completeness only, allows slicing of a scalar */
603 (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
604 numCopy++;
605 break;
606 case 1:
607 for (int i=region[0].first;i<region[0].second;i++) {
608 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset+numCopy];
609 numCopy++;
610 }
611 break;
612 case 2:
613 for (int j=region[1].first;j<region[1].second;j++) {
614 for (int i=region[0].first;i<region[0].second;i++) {
615 (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset+numCopy];
616 numCopy++;
617 }
618 }
619 break;
620 case 3:
621 for (int k=region[2].first;k<region[2].second;k++) {
622 for (int j=region[1].first;j<region[1].second;j++) {
623 for (int i=region[0].first;i<region[0].second;i++) {
624 (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset+numCopy];
625 numCopy++;
626 }
627 }
628 }
629 break;
630 case 4:
631 for (int l=region[3].first;l<region[3].second;l++) {
632 for (int k=region[2].first;k<region[2].second;k++) {
633 for (int j=region[1].first;j<region[1].second;j++) {
634 for (int i=region[0].first;i<region[0].second;i++) {
635 (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset+numCopy];
636 numCopy++;
637 }
638 }
639 }
640 }
641 break;
642 default:
643 stringstream mess;
644 mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
645 throw DataException(mess.str());
646 }
647
648 }
649
650 }
651
652 string
653 DataArrayView::toString(const string& suffix) const
654 {
655 EsysAssert(!isEmpty(),"Error - View is empty.");
656 stringstream temp;
657 string finalSuffix=suffix;
658 if (suffix.length() > 0) {
659 finalSuffix+=" ";
660 }
661 switch (getRank()) {
662 case 0:
663 temp << finalSuffix << (*this)();
664 break;
665 case 1:
666 for (int i=0;i<getShape()[0];i++) {
667 temp << finalSuffix << "(" << i << ") " << (*this)(i);
668 if (i!=(getShape()[0]-1)) {
669 temp << endl;
670 }
671 }
672 break;
673 case 2:
674 for (int i=0;i<getShape()[0];i++) {
675 for (int j=0;j<getShape()[1];j++) {
676 temp << finalSuffix << "(" << i << "," << j << ") " << (*this)(i,j);
677 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
678 temp << endl;
679 }
680 }
681 }
682 break;
683 case 3:
684 for (int i=0;i<getShape()[0];i++) {
685 for (int j=0;j<getShape()[1];j++) {
686 for (int k=0;k<getShape()[2];k++) {
687 temp << finalSuffix << "(" << i << "," << j << "," << k << ") " << (*this)(i,j,k);
688 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {
689 temp << endl;
690 }
691 }
692 }
693 }
694 break;
695 case 4:
696 for (int i=0;i<getShape()[0];i++) {
697 for (int j=0;j<getShape()[1];j++) {
698 for (int k=0;k<getShape()[2];k++) {
699 for (int l=0;l<getShape()[3];l++) {
700 temp << finalSuffix << "(" << i << "," << j << "," << k << "," << l << ") " << (*this)(i,j,k,l);
701 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {
702 temp << endl;
703 }
704 }
705 }
706 }
707 }
708 break;
709 default:
710 stringstream mess;
711 mess << "Error - (toString) Invalid rank: " << getRank();
712 throw DataException(mess.str());
713 }
714 return temp.str();
715 }
716
717 string
718 DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
719 {
720 stringstream temp;
721 temp << "(";
722 unsigned int i;
723 for (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