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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 122 - (hide annotations)
Thu Jun 9 05:38:05 2005 UTC (14 years, 4 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/DataArrayView.cpp
File size: 26681 byte(s)
Merge of development branch back to main trunk on 2005-06-09

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26