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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1455 - (hide annotations)
Thu Feb 28 17:19:44 2008 UTC (11 years, 8 months ago) by phornby
File size: 29337 byte(s)
Merge of branches/windows_from_1431_trunk.

Revamp of the exception system.
Fix unused vars and signed/unsigned comparisons.
defined a macro THROW(ARG) in the system_dep.h's to
deal with the expectations of declarations on different architectures.

Details in the logs of branches/windows_from_1431_trunk.

pre-merge snapshot of the trunk in tags/trunk_at_1452


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26