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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 474 - (hide annotations)
Mon Jan 30 04:23:44 2006 UTC (13 years, 8 months ago) by jgs
File size: 26489 byte(s)
restructure escript source tree
move src/Data/* -> src
remove inc
modify #includes and cpppath settings accordingly

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26