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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (show annotations)
Wed Sep 17 01:45:46 2008 UTC (10 years, 7 months ago) by jfenwick
File size: 30090 byte(s)
Merged noarrayview branch onto trunk.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26