/[escript]/trunk-mpi-branch/escript/src/DataArrayView.cpp
ViewVC logotype

Contents of /trunk-mpi-branch/escript/src/DataArrayView.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1306 - (show annotations)
Tue Sep 18 05:51:09 2007 UTC (11 years, 11 months ago) by ksteube
File size: 26114 byte(s)
New Copyright in each .c .h .cpp and .py file

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26