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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 478 - (hide annotations)
Tue Jan 31 02:21:49 2006 UTC (13 years, 8 months ago) by jgs
File size: 26645 byte(s)
rationalise #includes

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26