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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 113 - (hide annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/escript/src/Data/DataArrayView.cpp
File size: 26586 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26