/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1724 - (show annotations)
Mon Aug 25 05:38:57 2008 UTC (10 years, 11 months ago) by jfenwick
File size: 30044 byte(s)
Branch commit

Moved createShapeErrorMessage() into DataTypes.h
Modified functions in DataAlgorithm.h to use non-DataArrayView accesses.

Added getVector() to each of DataTagged, DataConstant, DataExpanded - This method returns 
the underlying DataVector by reference/constant reference. Note that this method does not 
exist in DataAbstract so (at the momement) in order to pull the data from something you 
need to know what you are looking at. (Lower level access is still possible via double* 
though).

DataTagged now has a getOffsetForTag method and a getDefaultOffset method.

DataMaths.h - A new file containing the reductionOps etc from DataArrayView (but without 
requiring DAV).
This file requires significant commenting improvements.


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));
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()));
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()));
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()==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()==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