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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1734 - (hide annotations)
Thu Aug 28 06:11:56 2008 UTC (11 years, 4 months ago) by jfenwick
File size: 30090 byte(s)
Added operations to Datatypes:
checkOffset
copySlice
copySliceFrom

Fixed some error reporting using EsysAssert.

Added two new c++ test suites:
DataTypesTest
DataMathsTest

Note that the test suite does not compile with dodebug=yes. There is an issue with linking one of the exception functions. I'm going to leave this 
until I have finished the rest of the work, perhaps Ken's scons changes will fix it.


1 jgs 82
2 ksteube 1312 /* $Id$ */
3 jgs 82
4 ksteube 1312 /*******************************************************
5 phornby 1455 *
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 ksteube 1312
16 jgs 478 #include "DataArrayView.h"
17     #include "DataException.h"
18    
19     #include <sstream>
20    
21     #include <boost/python/extract.hpp>
22    
23 jgs 82 using namespace std;
24     using namespace boost::python;
25 jfenwick 1697 using namespace escript::DataTypes;
26 jgs 82
27     namespace escript {
28    
29 phornby 1455 DataArrayView::DataArrayView():
30     m_shape(),
31     m_offset(0),
32     m_noValues(0)
33     {
34     m_data=0;
35     }
36 jgs 82
37 phornby 1455 DataArrayView::DataArrayView(ValueType& data,
38     const ShapeType& viewShape,
39     int offset):
40     m_data(&data),
41     m_shape(viewShape),
42     m_offset(offset),
43 jfenwick 1698 m_noValues(DataTypes::noValues(viewShape))
44 phornby 1455 {
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 jgs 82
64 phornby 1455 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 jgs 82
72 phornby 1628 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 phornby 1455 bool
83     DataArrayView::isEmpty() const
84     {
85 phornby 1628 return (m_data==NULL);
86 phornby 1455 }
87 jgs 82
88 phornby 1455 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 jgs 108
96 phornby 1455 EsysAssert((!isEmpty()&&checkShape(tempShape)),
97 jfenwick 1734 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape,m_shape));
98 jgs 82
99 phornby 1455 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 jgs 82 }
130 phornby 1455 }
131 jgs 102
132 phornby 1455 void
133     DataArrayView::copy(const DataArrayView& other)
134     {
135     copy(m_offset,other);
136     }
137 jgs 82
138 phornby 1455 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 jfenwick 1734 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
146 phornby 1455 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 jgs 82
153 phornby 1455 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 jfenwick 1734 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
162 phornby 1455 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 jgs 82
169 phornby 1455 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 jgs 82
183 phornby 1455 int
184     DataArrayView::getRank() const
185     {
186     return m_shape.size();
187     }
188 jgs 82
189 phornby 1455 const
190 jfenwick 1697 DataTypes::ShapeType&
191 phornby 1455 DataArrayView::getShape() const
192     {
193     return m_shape;
194     }
195 jgs 82
196 jfenwick 1698 // 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 jgs 102
222 phornby 1455 int
223     DataArrayView::noValues() const
224     {
225     return m_noValues;
226     }
227 jgs 82
228 phornby 1455 bool
229 jfenwick 1697 DataArrayView::checkShape(const DataTypes::ShapeType& other) const
230 phornby 1455 {
231     return (m_shape==other);
232     }
233 jgs 82
234 jfenwick 1724 // 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 jgs 108
245 jfenwick 1697 DataTypes::ValueType::size_type
246 phornby 1455 DataArrayView::getOffset() const
247     {
248     return m_offset;
249     }
250 jgs 108
251 phornby 1455 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 jgs 82
262 phornby 1455 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 jgs 108
273 phornby 1455 bool
274     DataArrayView::checkOffset() const
275     {
276     return checkOffset(m_offset);
277     }
278 jgs 108
279 phornby 1455 bool
280     DataArrayView::checkOffset(ValueType::size_type offset) const
281     {
282     return (m_data->size() >= (offset+noValues()));
283     }
284 jgs 108
285 jfenwick 1697 DataTypes::ValueType&
286 phornby 1455 DataArrayView::getData() const
287     {
288     EsysAssert(!isEmpty(),"Error - View is empty");
289     return *m_data;
290     }
291 jgs 108
292 jfenwick 1697 DataTypes::ValueType::reference
293 phornby 1455 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 jgs 108
301 jfenwick 1714 // 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 jgs 82
316 jfenwick 1714 // 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 jgs 102
352 jfenwick 1715 // 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 jgs 102
393 jfenwick 1715 // 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 jgs 108
410 phornby 1455 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 jgs 82
419 phornby 1455 void
420     DataArrayView::copySlice(ValueType::size_type thisOffset,
421     const DataArrayView& other,
422     ValueType::size_type otherOffset,
423     const RegionLoopRangeType& region)
424     {
425 jgs 102
426 phornby 1455 //
427     // Make sure views are not empty
428 jgs 102
429 phornby 1455 EsysAssert(!isEmpty(),
430     "Error - this view is empty.");
431     EsysAssert(!other.isEmpty(),
432     "Error - other view is empty.");
433 jgs 108
434 phornby 1455 //
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 jgs 102
438 phornby 1455 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 jgs 102
443 phornby 1455 EsysAssert(other.getRank()==region.size(),
444     "Error - slice not same rank as view to be sliced from.");
445 jgs 102
446 jfenwick 1734 EsysAssert(noValues()==DataTypes::noValues(getResultSliceShape(region)),
447 phornby 1455 "Error - slice shape not compatible shape for this view.");
448 jgs 102
449 phornby 1455 //
450     // copy the values in the specified region of the other view into this view
451 jgs 102
452 phornby 1455 // the following loops cannot be parallelised due to the numCopy counter
453     int numCopy=0;
454 jgs 108
455 phornby 1455 switch (region.size()) {
456     case 0:
457     /* this case should never be encountered,
458 jgs 108 as python will never pass us an empty region.
459     here for completeness only, allows slicing of a scalar */
460 phornby 1455 (*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 jgs 82 }
504 phornby 1455 }
505 jgs 102
506 phornby 1455 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 jgs 82
515 phornby 1455 void
516     DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
517     const DataArrayView& other,
518     ValueType::size_type otherOffset,
519     const RegionLoopRangeType& region)
520     {
521 jgs 102
522 phornby 1455 //
523     // Make sure views are not empty
524 jgs 102
525 phornby 1455 EsysAssert(!isEmpty(),
526     "Error - this view is empty.");
527     EsysAssert(!other.isEmpty(),
528     "Error - other view is empty.");
529 jgs 108
530 phornby 1455 //
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 jgs 102
534 phornby 1455 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 jgs 102
539 phornby 1455 EsysAssert(getRank()==region.size(),
540     "Error - slice not same rank as this view.");
541 jgs 102
542 jfenwick 1734 EsysAssert(other.getRank()==0 || other.noValues()==DataTypes::noValues(getResultSliceShape(region)),
543 phornby 1455 "Error - slice shape not compatible shape for other view.");
544 jgs 102
545 phornby 1455 //
546     // copy the values in the other view into the specified region of this view
547 jgs 102
548 phornby 1455 // allow for case where other view is a scalar
549     if (other.getRank()==0) {
550 jgs 108
551 phornby 1455 // the following loops cannot be parallelised due to the numCopy counter
552     int numCopy=0;
553 jgs 108
554 phornby 1455 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 jgs 108
604 phornby 1455 } else {
605 jgs 108
606 phornby 1455 // the following loops cannot be parallelised due to the numCopy counter
607     int numCopy=0;
608 jgs 108
609 phornby 1455 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 jgs 108
659 phornby 1455 }
660 jgs 108
661 phornby 1455 }
662 jgs 102
663 phornby 1455 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 jgs 82 }
672 phornby 1455 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 jgs 102 }
682 phornby 1455 }
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 jgs 102 }
725 phornby 1455 return temp.str();
726     }
727 jgs 82
728 jfenwick 1698 // 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 jgs 82
744 phornby 1455 void
745     DataArrayView::matMult(const DataArrayView& left,
746     const DataArrayView& right,
747     DataArrayView& result)
748     {
749 jgs 113
750 phornby 1455 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 jgs 82
756 phornby 1455 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 jgs 113
765 phornby 1455 int outputRank = left.getRank()+right.getRank()-2;
766 jgs 126
767 phornby 1455 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 jgs 82
774 phornby 1455 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 jgs 82
782 phornby 1455 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 jgs 82 }
790    
791 phornby 1455 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 jgs 82 }
802    
803 phornby 1455 switch (left.getRank()) {
804 jgs 113
805 jgs 82 case 1:
806 phornby 1455 switch (right.getRank()) {
807     case 1:
808 jgs 113 result()=0;
809     for (int i=0;i<left.getShape()[0];i++) {
810 phornby 1455 result()+=left(i)*right(i);
811 jgs 113 }
812     break;
813 phornby 1455 case 2:
814 jgs 113 for (int i=0;i<result.getShape()[0];i++) {
815 phornby 1455 result(i)=0;
816     for (int j=0;j<right.getShape()[0];j++) {
817     result(i)+=left(j)*right(j,i);
818     }
819 jgs 113 }
820     break;
821 phornby 1455 default:
822 jgs 113 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
823     throw DataException(temp.str());
824     break;
825 phornby 1455 }
826     break;
827 jgs 113
828 jgs 82 case 2:
829 phornby 1455 switch (right.getRank()) {
830     case 1:
831 jgs 113 result()=0;
832 phornby 1455 for (int i=0;i<left.getShape()[0];i++) {
833 jgs 113 result(i)=0;
834     for (int j=0;j<left.getShape()[1];j++) {
835 phornby 1455 result(i)+=left(i,j)*right(i);
836 jgs 113 }
837     }
838     break;
839 phornby 1455 case 2:
840 jgs 113 for (int i=0;i<result.getShape()[0];i++) {
841 phornby 1455 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 jgs 113 }
848     break;
849 phornby 1455 default:
850 jgs 113 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
851     throw DataException(temp.str());
852     break;
853 phornby 1455 }
854     break;
855 jgs 113
856 jgs 82 default:
857 phornby 1455 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << left.getRank();
858     throw DataException(temp.str());
859     break;
860     }
861 jgs 113
862 phornby 1455 }
863 jgs 82
864 jfenwick 1697 DataTypes::ShapeType
865 phornby 1455 DataArrayView::determineResultShape(const DataArrayView& left,
866     const DataArrayView& right)
867     {
868 phornby 1628 ShapeType result;
869 phornby 1455 for (int i=0; i<(left.getRank()-1); i++) {
870 phornby 1628 result.push_back(left.getShape()[i]);
871 phornby 1455 }
872     for (int i=1; i<right.getRank(); i++) {
873 phornby 1628 result.push_back(right.getShape()[i]);
874 phornby 1455 }
875 phornby 1628 return result;
876 phornby 1455 }
877 jgs 113
878 phornby 1455 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 jgs 82 }
889 phornby 1455 return result;
890     }
891 jgs 82
892 phornby 1455 bool operator!=(const DataArrayView& left, const DataArrayView& right)
893     {
894     return (!operator==(left,right));
895     }
896 jgs 82
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