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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 9 months ago) by phornby
Original Path: temp_trunk_copy/escript/src/DataArrayView.cpp
File size: 26114 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26