/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp
ViewVC logotype

Diff of /branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.1388  
changed lines
  Added in v.1696

  ViewVC Help
Powered by ViewVC 1.1.26