/[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 1698 by jfenwick, Tue Aug 12 01:13:16 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;
257    } else {        } else {
258      throw DataException("Error - invalid offset specified.");           throw DataException("Error - invalid offset specified.");
259    }        }
260  }     }
261    
262  void     void
263  DataArrayView::incrOffset()     DataArrayView::incrOffset()
264  {     {
265    EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");        EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
266    if (checkOffset(m_offset+noValues())) {        if (checkOffset(m_offset+noValues())) {
267      m_offset=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");  
268        } else {        } else {
269           /* calculate values for slice range */           throw DataException("Error - Cannot increment offset.");
270           for (i=0;i<slice_rank;i++) {        }
271             out[i]=getSliceRange(key[i],getShape()[i]);     }
272    
273       bool
274       DataArrayView::checkOffset() const
275       {
276          return checkOffset(m_offset);
277       }
278    
279       bool
280       DataArrayView::checkOffset(ValueType::size_type offset) const
281       {
282          return (m_data->size() >= (offset+noValues()));
283       }
284    
285       DataTypes::ValueType&
286       DataArrayView::getData() const
287       {
288          EsysAssert(!isEmpty(),"Error - View is empty");
289          return *m_data;
290       }
291    
292       DataTypes::ValueType::reference
293       DataArrayView::getData(ValueType::size_type i) const
294       {
295          //
296          // don't do any index checking so allow one past the end of the
297          // vector providing the equivalent of end()
298          return (*m_data)[m_offset+i];
299       }
300    
301       DataTypes::ShapeType
302       DataArrayView::getResultSliceShape(const RegionType& region)
303       {
304          int dimSize;
305          ShapeType result;
306          RegionType::const_iterator i;
307          for (i=region.begin();i!=region.end();i++) {
308             dimSize=((i->second)-(i->first));
309             if (dimSize!=0) {
310                result.push_back(dimSize);
311           }           }
312        }        }
313    } else {        return result;
314        slice_rank=1;     }
315        if (slice_rank>this_rank) {  
316             throw DataException("Error - rank of slices does not match rank of slicee");     DataTypes::RegionType
317       DataArrayView::getSliceRegion(const boost::python::object& key) const
318       {
319          int slice_rank, i;
320          int this_rank=getRank();
321          RegionType out(this_rank);
322          /* allow for case where key is singular eg: [1], this implies we
323          want to generate a rank-1 dimension object, as opposed to eg: [1,2]
324          which implies we want to take a rank dimensional object with one
325          dimension of size 1 */
326          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 {        } else {
339             out[0]=getSliceRange(key,getShape()[0]);           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    for (i=slice_rank;i<this_rank;i++) {     }
351      out[i]=std::pair<int,int>(0,getShape()[i]);  
352    }     std::pair<int,int>
353    return out;     getSliceRange(const boost::python::object& key,
354  }                   const int shape)
355       {
356  std::pair<int,int>        /* default slice range is range of entire shape dimension */
357  getSliceRange(const boost::python::object& key,        int s0=0, s1=shape;;
358                const int shape)        extract<int> slice_int(key);
359  {        if (slice_int.check()) {
360    /* default slice range is range of entire shape dimension */           /* if the key is a single int set start=key and end=key */
361    int s0=0, s1=shape;;           /* in this case, we want to return a rank-1 dimension object from
362    extract<int> slice_int(key);           this object, taken from a single index value for one of this
363    if (slice_int.check()) {           object's dimensions */
364      /* if the key is a single int set start=key and end=key */           s0=slice_int();
365      /* 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;  
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             throw DataException("Error - slice index out of range.");
384          if (s0 == s1 && s1 >= shape)
385             throw DataException("Error - slice index out of range.");
386          if (s0 != s1 &&  s1>shape)
387             throw DataException("Error - slice index out of range.");
388          if (s0 > s1)
389             throw DataException("Error - lower index must less or equal upper index.");
390          return std::pair<int,int>(s0,s1);
391       }
392    
393       DataTypes::RegionLoopRangeType
394       getSliceRegionLoopRange(const DataTypes::RegionType& region)
395       {
396          DataTypes::RegionLoopRangeType region_loop_range(region.size());
397          unsigned int i;
398          for (i=0;i<region.size();i++) {
399             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                region_loop_range[i].first=region[i].first;
404                region_loop_range[i].second=region[i].second;
405             }
406        }        }
407      }        return region_loop_range;
408      return region_loop_range;     }
409  }  
410       void
411  void     DataArrayView::copySlice(const DataArrayView& other,
412  DataArrayView::copySlice(const DataArrayView& other,                              const RegionLoopRangeType& region)
413                           const RegionLoopRangeType& region)     {
414  {        //
415    //        // version of copySlice that uses the default offsets
416    // version of copySlice that uses the default offsets        copySlice(m_offset,other,other.m_offset,region);
417    copySlice(m_offset,other,other.m_offset,region);     }
418  }  
419       void
420  void     DataArrayView::copySlice(ValueType::size_type thisOffset,
421  DataArrayView::copySlice(ValueType::size_type thisOffset,                              const DataArrayView& other,
422                           const DataArrayView& other,                              ValueType::size_type otherOffset,
423                           ValueType::size_type otherOffset,                              const RegionLoopRangeType& region)
424                           const RegionLoopRangeType& region)     {
425  {  
426          //
427      //        // Make sure views are not empty
428      // Make sure views are not empty  
429          EsysAssert(!isEmpty(),
430      EsysAssert(!isEmpty(),                   "Error - this view is empty.");
431                 "Error - this view is empty.");        EsysAssert(!other.isEmpty(),
432      EsysAssert(!other.isEmpty(),                   "Error - other view is empty.");
433                 "Error - other view is empty.");  
434          //
435      //        // Check the view to be sliced from is compatible with the region to be sliced,
436      // 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:
437      // and that the region to be sliced is compatible with this view:  
438          EsysAssert(checkOffset(thisOffset),
439      EsysAssert(checkOffset(thisOffset),                   "Error - offset incompatible with this view.");
440                 "Error - offset incompatible with this view.");        EsysAssert(otherOffset+noValues()<=other.getData().size(),
441      EsysAssert(otherOffset+noValues()<=other.getData().size(),                   "Error - offset incompatible with other view.");
442                 "Error - offset incompatible with other view.");  
443          EsysAssert(other.getRank()==region.size(),
444      EsysAssert(other.getRank()==region.size(),                   "Error - slice not same rank as view to be sliced from.");
445                 "Error - slice not same rank as view to be sliced from.");  
446          EsysAssert(noValues()==noValues(getResultSliceShape(region)),
447      EsysAssert(noValues()==noValues(getResultSliceShape(region)),                   "Error - slice shape not compatible shape for this view.");
448                 "Error - slice shape not compatible shape for this view.");  
449          //
450      //        // copy the values in the specified region of the other view into this view
451      // copy the values in the specified region of the other view into this view  
452          // the following loops cannot be parallelised due to the numCopy counter
453      // the following loops cannot be parallelised due to the numCopy counter        int numCopy=0;
454      int numCopy=0;  
455          switch (region.size()) {
456      switch (region.size()) {        case 0:
457      case 0:           /* this case should never be encountered,
       /* 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.1698

  ViewVC Help
Powered by ViewVC 1.1.26