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

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

  ViewVC Help
Powered by ViewVC 1.1.26