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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 757 - (hide annotations)
Mon Jun 26 13:12:56 2006 UTC (13 years, 3 months ago) by woo409
File size: 26248 byte(s)
+ Merge of intelc_win32 branch (revision 741:755) with trunk. Tested on iVEC altix (run_tests and py_tests all pass)

1 jgs 82 // $Id$
2    
3     /*
4 elspeth 615 ************************************************************
5     * Copyright 2006 by ACcESS MNRF *
6     * *
7     * http://www.access.edu.au *
8     * Primary Business: Queensland, Australia *
9     * Licensed under the Open Software License version 3.0 *
10     * http://www.opensource.org/licenses/osl-3.0.php *
11     * *
12     ************************************************************
13 jgs 82 */
14    
15 jgs 478 #include "DataArrayView.h"
16     #include "DataException.h"
17    
18     #include <sstream>
19    
20     #include <boost/python/extract.hpp>
21    
22 jgs 82 using namespace std;
23     using namespace boost::python;
24    
25     namespace escript {
26    
27     DataArrayView::DataArrayView():
28     m_shape(),
29 jgs 108 m_offset(0),
30     m_noValues(0)
31 jgs 82 {
32 jgs 122 m_data=0;
33 jgs 82 }
34    
35     DataArrayView::DataArrayView(ValueType& data,
36     const ShapeType& viewShape,
37     int offset):
38     m_data(&data),
39     m_shape(viewShape),
40     m_offset(offset),
41     m_noValues(noValues(viewShape))
42     {
43     //
44 jgs 108 // check the shape rank and size and throw an exception if a
45     // shape incompatible with the given data array is specified
46 jgs 82 if (viewShape.size() > m_maxRank) {
47     stringstream temp;
48 jgs 108 temp << "Error- Couldn't construct DataArrayView, invalid rank."
49     << "Input data rank: " << viewShape.size() << " "
50     << "Maximum rank allowed for DataArrayView: " << m_maxRank;
51 jgs 82 throw DataException(temp.str());
52     }
53 jgs 108 if (m_data->size() < (m_offset+noValues())) {
54 jgs 82 stringstream temp;
55 jgs 108 temp << "Error- Couldn't construct DataArrayView, insufficient data values. "
56     << "Shape requires: " << noValues()+m_offset << " values."
57     << "Data values available: " << m_data->size();
58 jgs 82 throw DataException(temp.str());
59     }
60     }
61    
62     DataArrayView::DataArrayView(const DataArrayView& other):
63     m_data(other.m_data),
64     m_offset(other.m_offset),
65     m_shape(other.m_shape),
66     m_noValues(other.m_noValues)
67     {
68     }
69    
70     bool
71     DataArrayView::isEmpty() const
72     {
73 jgs 122 return (m_data==0);
74 jgs 82 }
75    
76     void
77     DataArrayView::copy(const boost::python::numeric::array& value)
78     {
79 jgs 117 ShapeType tempShape;
80 jgs 108 for (int i=0; i<value.getrank(); i++) {
81 jgs 82 tempShape.push_back(extract<int>(value.getshape()[i]));
82     }
83 jgs 108
84 jgs 82 EsysAssert((!isEmpty()&&checkShape(tempShape)),
85     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
86    
87     if (value.getrank()==0) {
88     (*this)()=extract<double>(value[value.getshape()]);
89     } else if (value.getrank()==1) {
90 jgs 122 for (ValueType::size_type i=0;i<tempShape[0];i++) {
91 jgs 82 (*this)(i)=extract<double>(value[i]);
92     }
93     } else if (value.getrank()==2) {
94 jgs 122 for (ValueType::size_type i=0;i<tempShape[0];i++) {
95     for (ValueType::size_type j=0;j<tempShape[1];j++) {
96 jgs 82 (*this)(i,j)=extract<double>(value[i][j]);
97     }
98     }
99     } else if (value.getrank()==3) {
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     for (ValueType::size_type k=0;k<tempShape[2];k++) {
103 jgs 82 (*this)(i,j,k)=extract<double>(value[i][j][k]);
104     }
105     }
106     }
107     } else if (value.getrank()==4) {
108 jgs 122 for (ValueType::size_type i=0;i<tempShape[0];i++) {
109     for (ValueType::size_type j=0;j<tempShape[1];j++) {
110     for (ValueType::size_type k=0;k<tempShape[2];k++) {
111     for (ValueType::size_type l=0;l<tempShape[3];l++) {
112     (*this)(i,j,k,l)=extract<double>(value[i][j][k][l]);
113 jgs 82 }
114     }
115     }
116     }
117     }
118     }
119 jgs 102
120 jgs 82 void
121     DataArrayView::copy(const DataArrayView& other)
122     {
123 jgs 122 copy(m_offset,other);
124 jgs 82 }
125    
126     void
127     DataArrayView::copy(ValueType::size_type offset,
128     const DataArrayView& other)
129     {
130 jgs 108 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
131     "Error - Couldn't copy due to insufficient storage.");
132     EsysAssert((checkShape(other.getShape())),
133 jgs 82 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
134 jgs 108 if (checkOffset(offset)) {
135     memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
136     } else {
137     throw DataException("Error - invalid offset specified.");
138     }
139 jgs 82 }
140    
141     void
142     DataArrayView::copy(ValueType::size_type thisOffset,
143     const DataArrayView& other,
144     ValueType::size_type otherOffset)
145     {
146 jgs 108 EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
147     "Error - Couldn't copy due to insufficient storage.");
148     EsysAssert((checkShape(other.getShape())),
149 jgs 82 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
150 jgs 108 if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
151     memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
152     } else {
153     throw DataException("Error - invalid offset specified.");
154     }
155 jgs 82 }
156    
157     void
158     DataArrayView::copy(ValueType::size_type offset,
159     ValueType::value_type value)
160     {
161 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
162     "Error - Couldn't copy due to insufficient storage.");
163     if (checkOffset(offset)) {
164 jgs 122 vector<double> temp(noValues(),value);
165 jgs 108 memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
166     } else {
167     throw DataException("Error - invalid offset specified.");
168     }
169 jgs 82 }
170    
171     int
172     DataArrayView::getRank() const
173     {
174 jgs 122 return m_shape.size();
175 jgs 82 }
176    
177 jgs 117 const
178     DataArrayView::ShapeType&
179 jgs 82 DataArrayView::getShape() const
180     {
181 jgs 122 return m_shape;
182 jgs 82 }
183    
184     int
185     DataArrayView::noValues(const ShapeType& shape)
186     {
187 jgs 122 ShapeType::const_iterator i;
188     //
189     // An empty shape vector means rank 0 which contains 1 value
190     int noValues=1;
191     for (i=shape.begin();i!=shape.end();i++) {
192     noValues*=(*i);
193     }
194     return noValues;
195 jgs 82 }
196 jgs 108
197     int
198     DataArrayView::noValues(const RegionLoopRangeType& region)
199     {
200 jgs 122 //
201     // An empty region vector means rank 0 which contains 1 value
202     int noValues=1;
203     for (int i=0;i<region.size();i++) {
204     noValues*=region[i].second-region[i].first;
205     }
206     return noValues;
207 jgs 108 }
208 jgs 102
209 jgs 82 int
210     DataArrayView::noValues() const
211     {
212 jgs 122 return m_noValues;
213 jgs 82 }
214    
215     bool
216     DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
217     {
218 jgs 122 return (m_shape==other);
219 jgs 82 }
220    
221     string
222     DataArrayView::createShapeErrorMessage(const string& messagePrefix,
223     const DataArrayView::ShapeType& other) const
224     {
225 jgs 122 stringstream temp;
226     temp << messagePrefix
227     << " This shape: " << shapeToString(m_shape)
228     << " Other shape: " << shapeToString(other);
229     return temp.str();
230 jgs 82 }
231 jgs 108
232     DataArrayView::ValueType::size_type
233     DataArrayView::getOffset() const
234     {
235     return m_offset;
236     }
237    
238 jgs 82 void
239     DataArrayView::setOffset(ValueType::size_type offset)
240     {
241 jgs 122 EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
242     if (checkOffset(offset)) {
243     m_offset=offset;
244     } else {
245     throw DataException("Error - invalid offset specified.");
246     }
247 jgs 82 }
248    
249 jgs 108 void
250     DataArrayView::incrOffset()
251     {
252 jgs 122 EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
253     if (checkOffset(m_offset+noValues())) {
254     m_offset=m_offset+noValues();
255     } else {
256     throw DataException("Error - Cannot increment offset.");
257     }
258 jgs 108 }
259    
260     bool
261     DataArrayView::checkOffset() const
262     {
263 jgs 122 return checkOffset(m_offset);
264 jgs 108 }
265    
266     bool
267     DataArrayView::checkOffset(ValueType::size_type offset) const
268     {
269 jgs 122 return (m_data->size() >= (offset+noValues()));
270 jgs 108 }
271    
272     DataArrayView::ValueType&
273     DataArrayView::getData() const
274     {
275     EsysAssert(!isEmpty(),"Error - View is empty");
276     return *m_data;
277     }
278    
279     DataArrayView::ValueType::reference
280     DataArrayView::getData(ValueType::size_type i) const
281     {
282     //
283     // don't do any index checking so allow one past the end of the
284     // vector providing the equivalent of end()
285     return (*m_data)[m_offset+i];
286     }
287    
288 jgs 82 DataArrayView::ShapeType
289     DataArrayView::getResultSliceShape(const RegionType& region)
290     {
291 jgs 122 int dimSize;
292     RegionType::const_iterator i;
293     ShapeType result;
294     for (i=region.begin();i!=region.end();i++) {
295     dimSize=((i->second)-(i->first));
296     if (dimSize!=0) {
297     result.push_back(dimSize);
298 jgs 82 }
299 jgs 122 }
300     return result;
301 jgs 82 }
302    
303 jgs 102 DataArrayView::RegionType
304     DataArrayView::getSliceRegion(const boost::python::object& key) const
305     {
306 jgs 108 int slice_rank, i;
307 jgs 102 int this_rank=getRank();
308 jgs 117 RegionType out(this_rank);
309 jgs 102 /* allow for case where key is singular eg: [1], this implies we
310     want to generate a rank-1 dimension object, as opposed to eg: [1,2]
311     which implies we want to take a rank dimensional object with one
312     dimension of size 1 */
313     extract<tuple> key_tuple(key);
314     if (key_tuple.check()) {
315     slice_rank=extract<int> (key.attr("__len__")());
316 jgs 108 /* ensure slice is correctly dimensioned */
317     if (slice_rank>this_rank) {
318     throw DataException("Error - rank of slices does not match rank of slicee");
319     } else {
320     /* calculate values for slice range */
321     for (i=0;i<slice_rank;i++) {
322     out[i]=getSliceRange(key[i],getShape()[i]);
323     }
324     }
325 jgs 102 } else {
326     slice_rank=1;
327 jgs 108 if (slice_rank>this_rank) {
328     throw DataException("Error - rank of slices does not match rank of slicee");
329     } else {
330     out[0]=getSliceRange(key,getShape()[0]);
331     }
332 jgs 102 }
333     for (i=slice_rank;i<this_rank;i++) {
334     out[i]=std::pair<int,int>(0,getShape()[i]);
335     }
336     return out;
337     }
338    
339     std::pair<int,int>
340     getSliceRange(const boost::python::object& key,
341     const int shape)
342     {
343     /* default slice range is range of entire shape dimension */
344     int s0=0, s1=shape;;
345     extract<int> slice_int(key);
346     if (slice_int.check()) {
347     /* if the key is a single int set start=key and end=key */
348     /* in this case, we want to return a rank-1 dimension object from
349     this object, taken from a single index value for one of this
350     object's dimensions */
351     s0=slice_int();
352     s1=s0;
353     } else {
354     /* if key is a pair extract begin and end values */
355     extract<int> step(key.attr("step"));
356     if (step.check() && step()!=1) {
357     throw DataException("Error - Data does not support increments in slicing ");
358     } else {
359     extract<int> start(key.attr("start"));
360     if (start.check()) {
361     s0=start();
362     }
363     extract<int> stop(key.attr("stop"));
364     if (stop.check()) {
365     s1=stop();
366     }
367     }
368     }
369 jgs 108 if (s0 < 0)
370     throw DataException("Error - slice index out of range.");
371 woo409 757 if (s0 == s1 && s1 >= shape)
372 jgs 108 throw DataException("Error - slice index out of range.");
373 woo409 757 if (s0 != s1 && s1>shape)
374 jgs 108 throw DataException("Error - slice index out of range.");
375     if (s0 > s1)
376     throw DataException("Error - lower index must less or equal upper index.");
377 jgs 102 return std::pair<int,int>(s0,s1);
378     }
379    
380 jgs 108 DataArrayView::RegionLoopRangeType
381     getSliceRegionLoopRange(const DataArrayView::RegionType& region)
382     {
383     DataArrayView::RegionLoopRangeType region_loop_range(region.size());
384     for (int i=0;i<region.size();i++) {
385     if (region[i].first==region[i].second) {
386     region_loop_range[i].first=region[i].first;
387     region_loop_range[i].second=region[i].second+1;
388     } else {
389     region_loop_range[i].first=region[i].first;
390     region_loop_range[i].second=region[i].second;
391     }
392     }
393     return region_loop_range;
394     }
395    
396 jgs 82 void
397     DataArrayView::copySlice(const DataArrayView& other,
398 jgs 108 const RegionLoopRangeType& region)
399 jgs 82 {
400 jgs 122 //
401     // version of copySlice that uses the default offsets
402     copySlice(m_offset,other,other.m_offset,region);
403 jgs 82 }
404    
405     void
406     DataArrayView::copySlice(ValueType::size_type thisOffset,
407     const DataArrayView& other,
408     ValueType::size_type otherOffset,
409 jgs 108 const RegionLoopRangeType& region)
410 jgs 82 {
411 jgs 102
412 jgs 82 //
413 jgs 108 // Make sure views are not empty
414 jgs 102
415 jgs 108 EsysAssert(!isEmpty(),
416     "Error - this view is empty.");
417     EsysAssert(!other.isEmpty(),
418     "Error - other view is empty.");
419    
420 jgs 82 //
421 jgs 108 // Check the view to be sliced from is compatible with the region to be sliced,
422     // and that the region to be sliced is compatible with this view:
423 jgs 102
424 jgs 108 EsysAssert(checkOffset(thisOffset),
425     "Error - offset incompatible with this view.");
426     EsysAssert(otherOffset+noValues()<=other.getData().size(),
427     "Error - offset incompatible with other view.");
428 jgs 102
429 jgs 108 EsysAssert(other.getRank()==region.size(),
430     "Error - slice not same rank as view to be sliced from.");
431 jgs 102
432 jgs 108 EsysAssert(noValues()==noValues(getResultSliceShape(region)),
433     "Error - slice shape not compatible shape for this view.");
434 jgs 102
435     //
436 jgs 108 // copy the values in the specified region of the other view into this view
437 jgs 102
438 jgs 122 // the following loops cannot be parallelised due to the numCopy counter
439 jgs 82 int numCopy=0;
440 jgs 108
441     switch (region.size()) {
442 jgs 82 case 0:
443 jgs 108 /* this case should never be encountered,
444     as python will never pass us an empty region.
445     here for completeness only, allows slicing of a scalar */
446     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
447     numCopy++;
448 jgs 82 break;
449     case 1:
450 jgs 108 for (int i=region[0].first;i<region[0].second;i++) {
451     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
452     numCopy++;
453 jgs 82 }
454     break;
455     case 2:
456 jgs 108 for (int j=region[1].first;j<region[1].second;j++) {
457     for (int i=region[0].first;i<region[0].second;i++) {
458     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
459 jgs 102 numCopy++;
460 jgs 82 }
461     }
462     break;
463     case 3:
464 jgs 108 for (int k=region[2].first;k<region[2].second;k++) {
465     for (int j=region[1].first;j<region[1].second;j++) {
466     for (int i=region[0].first;i<region[0].second;i++) {
467     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
468 jgs 102 numCopy++;
469 jgs 82 }
470     }
471     }
472     break;
473     case 4:
474 jgs 108 for (int l=region[3].first;l<region[3].second;l++) {
475     for (int k=region[2].first;k<region[2].second;k++) {
476     for (int j=region[1].first;j<region[1].second;j++) {
477     for (int i=region[0].first;i<region[0].second;i++) {
478     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];
479 jgs 102 numCopy++;
480 jgs 82 }
481     }
482     }
483     }
484     break;
485     default:
486     stringstream mess;
487 jgs 108 mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
488 jgs 82 throw DataException(mess.str());
489     }
490     }
491 jgs 102
492 jgs 82 void
493     DataArrayView::copySliceFrom(const DataArrayView& other,
494 jgs 108 const RegionLoopRangeType& region_loop_range)
495 jgs 82 {
496     //
497 jgs 108 // version of copySliceFrom that uses the default offsets
498     copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
499 jgs 82 }
500    
501     void
502     DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
503     const DataArrayView& other,
504     ValueType::size_type otherOffset,
505 jgs 108 const RegionLoopRangeType& region)
506 jgs 82 {
507 jgs 102
508 jgs 82 //
509 jgs 108 // Make sure views are not empty
510 jgs 102
511 jgs 108 EsysAssert(!isEmpty(),
512     "Error - this view is empty.");
513     EsysAssert(!other.isEmpty(),
514     "Error - other view is empty.");
515    
516 jgs 82 //
517 jgs 108 // Check this view is compatible with the region to be sliced,
518     // and that the region to be sliced is compatible with the other view:
519 jgs 102
520 jgs 108 EsysAssert(other.checkOffset(otherOffset),
521     "Error - offset incompatible with other view.");
522     EsysAssert(thisOffset+other.noValues()<=getData().size(),
523     "Error - offset incompatible with this view.");
524 jgs 102
525 jgs 108 EsysAssert(getRank()==region.size(),
526     "Error - slice not same rank as this view.");
527 jgs 102
528 jgs 108 EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
529     "Error - slice shape not compatible shape for other view.");
530 jgs 102
531     //
532 jgs 108 // copy the values in the other view into the specified region of this view
533 jgs 102
534 jgs 108 // allow for case where other view is a scalar
535     if (other.getRank()==0) {
536    
537 jgs 122 // the following loops cannot be parallelised due to the numCopy counter
538 jgs 108 int numCopy=0;
539    
540     switch (region.size()) {
541     case 0:
542     /* this case should never be encountered,
543     as python will never pass us an empty region.
544     here for completeness only, allows slicing of a scalar */
545     (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
546 jgs 102 numCopy++;
547 jgs 108 break;
548     case 1:
549     for (int i=region[0].first;i<region[0].second;i++) {
550     (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
551 jgs 102 numCopy++;
552 jgs 108 }
553     break;
554     case 2:
555     for (int j=region[1].first;j<region[1].second;j++) {
556     for (int i=region[0].first;i<region[0].second;i++) {
557     (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
558 jgs 102 numCopy++;
559 jgs 82 }
560 jgs 108 }
561     break;
562     case 3:
563     for (int k=region[2].first;k<region[2].second;k++) {
564     for (int j=region[1].first;j<region[1].second;j++) {
565     for (int i=region[0].first;i<region[0].second;i++) {
566     (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
567     numCopy++;
568     }
569     }
570     }
571     break;
572     case 4:
573     for (int l=region[3].first;l<region[3].second;l++) {
574     for (int k=region[2].first;k<region[2].second;k++) {
575     for (int j=region[1].first;j<region[1].second;j++) {
576     for (int i=region[0].first;i<region[0].second;i++) {
577     (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
578     numCopy++;
579     }
580     }
581     }
582     }
583     break;
584     default:
585     stringstream mess;
586     mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
587     throw DataException(mess.str());
588     }
589    
590     } else {
591    
592 jgs 122 // the following loops cannot be parallelised due to the numCopy counter
593 jgs 108 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 jgs 126 if (left.getShape()[left.getRank()-1] != right.getShape()[0]) {
742 jgs 82 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 126 int outputRank = left.getRank()+right.getRank()-2;
751    
752 jgs 82 if (outputRank < 0) {
753     stringstream temp;
754 jgs 126 temp << "Error - (matMult) LHS and RHS cannot be multiplied "
755 jgs 82 << "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 jgs 126 for (int i=0; i<(left.getRank()-1); i++) {
768     if (left.getShape()[i] != result.getShape()[i]) {
769 jgs 82 stringstream temp;
770     temp << "Error - (matMult) Dimension: " << i
771 jgs 126 << " of LHS and result array don't match.";
772 jgs 82 throw DataException(temp.str());
773     }
774     }
775    
776 jgs 126 for (int i=1; i<right.getRank(); i++) {
777 jgs 82 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 jgs 126 << " of result array don't match.";
784 jgs 82 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 jgs 117 ShapeType temp;
854 jgs 113 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