/[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 1734 by jfenwick, Thu Aug 28 06:11:56 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,m_shape));
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      }     }
 }  
   
 void  
 DataArrayView::copy(const DataArrayView& other)  
 {  
   copy(m_offset,other);  
 }  
   
 void  
 DataArrayView::copy(ValueType::size_type offset,  
                     const DataArrayView& other)  
 {  
     EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),  
            "Error - Couldn't copy due to insufficient storage.");  
     EsysAssert((checkShape(other.getShape())),  
            createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));  
     if (checkOffset(offset)) {  
       memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());  
     } else {  
       throw DataException("Error - invalid offset specified.");  
     }  
 }  
   
 void  
 DataArrayView::copy(ValueType::size_type thisOffset,  
                     const DataArrayView& other,  
                     ValueType::size_type otherOffset)  
 {  
     EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),  
            "Error - Couldn't copy due to insufficient storage.");  
     EsysAssert((checkShape(other.getShape())),  
            createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));  
     if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {  
       memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());  
     } else {  
       throw DataException("Error - invalid offset specified.");  
     }  
 }  
   
 void  
 DataArrayView::copy(ValueType::size_type offset,  
                     ValueType::value_type value)  
 {  
     EsysAssert((!isEmpty()&&checkOffset(offset)),  
            "Error - Couldn't copy due to insufficient storage.");  
     if (checkOffset(offset)) {  
       vector<double> temp(noValues(),value);  
       memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());  
     } else {  
       throw DataException("Error - invalid offset specified.");  
     }  
 }  
   
 int  
 DataArrayView::getRank() const  
 {  
   return m_shape.size();  
 }  
   
 const  
 DataArrayView::ShapeType&  
 DataArrayView::getShape() const  
 {  
   return m_shape;  
 }  
   
 int  
 DataArrayView::noValues(const ShapeType& shape)  
 {  
   ShapeType::const_iterator i;  
   //  
   // An empty shape vector means rank 0 which contains 1 value  
   int noValues=1;  
   for (i=shape.begin();i!=shape.end();i++) {  
     noValues*=(*i);  
   }  
   return noValues;  
 }  
   
 int  
 DataArrayView::noValues(const RegionLoopRangeType& region)  
 {  
   //  
   // An empty region vector means rank 0 which contains 1 value  
   int noValues=1;  
   for (int i=0;i<region.size();i++) {  
     noValues*=region[i].second-region[i].first;  
   }  
   return noValues;  
 }  
131    
132  int     void
133  DataArrayView::noValues() const     DataArrayView::copy(const DataArrayView& other)
134  {     {
135    return m_noValues;        copy(m_offset,other);
136  }     }
137    
138  bool     void
139  DataArrayView::checkShape(const DataArrayView::ShapeType& other) const     DataArrayView::copy(ValueType::size_type offset,
140  {                         const DataArrayView& other)
141    return (m_shape==other);     {
142  }        EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
143                     "Error - Couldn't copy due to insufficient storage.");
144  string        EsysAssert((checkShape(other.getShape())),
145  DataArrayView::createShapeErrorMessage(const string& messagePrefix,                   createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
146                                         const DataArrayView::ShapeType& other) const        if (checkOffset(offset)) {
147  {           memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
   stringstream temp;  
   temp << messagePrefix  
        << " This shape: " << shapeToString(m_shape)  
        << " Other shape: " << shapeToString(other);  
   return temp.str();  
 }  
   
 DataArrayView::ValueType::size_type  
 DataArrayView::getOffset() const  
 {  
   return m_offset;  
 }  
   
 void  
 DataArrayView::setOffset(ValueType::size_type offset)  
 {  
   EsysAssert((checkOffset(offset)), "Error - Invalid offset.");  
   if (checkOffset(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");  
148        } else {        } else {
149           /* 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]);  
          }  
150        }        }
151    } else {     }
152        slice_rank=1;  
153        if (slice_rank>this_rank) {     void
154             throw DataException("Error - rank of slices does not match rank of slicee");     DataArrayView::copy(ValueType::size_type thisOffset,
155                           const DataArrayView& other,
156                           ValueType::size_type otherOffset)
157       {
158          EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
159                     "Error - Couldn't copy due to insufficient storage.");
160          EsysAssert((checkShape(other.getShape())),
161                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
162          if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
163             memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
164        } else {        } else {
165             out[0]=getSliceRange(key,getShape()[0]);           throw DataException("Error - invalid offset specified.");
166        }        }
167    }     }
168    for (i=slice_rank;i<this_rank;i++) {  
169      out[i]=std::pair<int,int>(0,getShape()[i]);     void
170    }     DataArrayView::copy(ValueType::size_type offset,
171    return out;                         ValueType::value_type value)
172  }     {
173          EsysAssert((!isEmpty()&&checkOffset(offset)),
174  std::pair<int,int>                   "Error - Couldn't copy due to insufficient storage.");
175  getSliceRange(const boost::python::object& key,        if (checkOffset(offset)) {
176                const int shape)           vector<double> temp(noValues(),value);
177  {           memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
   /* default slice range is range of entire shape dimension */  
   int s0=0, s1=shape;;  
   extract<int> slice_int(key);  
   if (slice_int.check()) {  
     /* if the key is a single int set start=key and end=key */  
     /* in this case, we want to return a rank-1 dimension object from  
        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;  
178        } else {        } else {
179          region_loop_range[i].first=region[i].first;           throw DataException("Error - invalid offset specified.");
         region_loop_range[i].second=region[i].second;  
180        }        }
181      }     }
182      return region_loop_range;  
183  }     int
184       DataArrayView::getRank() const
185  void     {
186  DataArrayView::copySlice(const DataArrayView& other,        return m_shape.size();
187                           const RegionLoopRangeType& region)     }
188  {  
189    //     const
190    // version of copySlice that uses the default offsets     DataTypes::ShapeType&
191    copySlice(m_offset,other,other.m_offset,region);     DataArrayView::getShape() const
192  }     {
193          return m_shape;
194  void     }
195  DataArrayView::copySlice(ValueType::size_type thisOffset,  
196                           const DataArrayView& other,  //    int
197                           ValueType::size_type otherOffset,  //    DataArrayView::noValues(const ShapeType& shape)
198                           const RegionLoopRangeType& region)  //    {
199  {  //       ShapeType::const_iterator i;
200    //       //
201      //  //       // An empty shape vector means rank 0 which contains 1 value
202      // Make sure views are not empty  //       int noValues=1;
203    //       for (i=shape.begin();i!=shape.end();i++) {
204      EsysAssert(!isEmpty(),  //          noValues*=(*i);
205                 "Error - this view is empty.");  //       }
206      EsysAssert(!other.isEmpty(),  //       return noValues;
207                 "Error - other view is empty.");  //    }
208    //
209      //  //    int
210      // Check the view to be sliced from is compatible with the region to be sliced,  //    DataArrayView::noValues(const RegionLoopRangeType& region)
211      // and that the region to be sliced is compatible with this view:  //    {
212    //       //
213      EsysAssert(checkOffset(thisOffset),  //       // An empty region vector means rank 0 which contains 1 value
214                 "Error - offset incompatible with this view.");  //       int noValues=1;
215      EsysAssert(otherOffset+noValues()<=other.getData().size(),  //       unsigned int i;
216                 "Error - offset incompatible with other view.");  //       for (i=0;i<region.size();i++) {
217    //          noValues*=region[i].second-region[i].first;
218      EsysAssert(other.getRank()==region.size(),  //       }
219                 "Error - slice not same rank as view to be sliced from.");  //       return noValues;
220    //    }
     EsysAssert(noValues()==noValues(getResultSliceShape(region)),  
                "Error - slice shape not compatible shape for this view.");  
   
     //  
     // copy the values in the specified region of the other view into this view  
   
     // 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+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];  
       numCopy++;  
       break;  
     case 1:  
       for (int i=region[0].first;i<region[0].second;i++) {  
     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];  
         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+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];  
       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+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];  
         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+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];  
           numCopy++;  
         }  
       }  
     }  
       }  
       break;  
     default:  
       stringstream mess;  
       mess << "Error - (copySlice) Invalid slice region rank: " << region.size();  
       throw DataException(mess.str());  
     }  
 }  
221    
222  void     int
223  DataArrayView::copySliceFrom(const DataArrayView& other,     DataArrayView::noValues() const
224                               const RegionLoopRangeType& region_loop_range)     {
225  {        return m_noValues;
226      //     }
227      // version of copySliceFrom that uses the default offsets  
228      copySliceFrom(m_offset,other,other.m_offset,region_loop_range);     bool
229  }     DataArrayView::checkShape(const DataTypes::ShapeType& other) const
230       {
231  void        return (m_shape==other);
232  DataArrayView::copySliceFrom(ValueType::size_type thisOffset,     }
233                               const DataArrayView& other,  
234                               ValueType::size_type otherOffset,  //    string
235                               const RegionLoopRangeType& region)  //    DataArrayView::createShapeErrorMessage(const string& messagePrefix,
236  {  //                                           const DataTypes::ShapeType& other) const
237    //    {
238      //  //       stringstream temp;
239      // Make sure views are not empty  //       temp << messagePrefix
240    //            << " This shape: " << shapeToString(m_shape)
241      EsysAssert(!isEmpty(),  //            << " Other shape: " << shapeToString(other);
242                 "Error - this view is empty.");  //       return temp.str();
243      EsysAssert(!other.isEmpty(),  //    }
244                 "Error - other view is empty.");  
245       DataTypes::ValueType::size_type
246      //     DataArrayView::getOffset() const
247      // Check this view is compatible with the region to be sliced,     {
248      // and that the region to be sliced is compatible with the other view:        return m_offset;
249       }
250      EsysAssert(other.checkOffset(otherOffset),  
251                 "Error - offset incompatible with other view.");     void
252      EsysAssert(thisOffset+other.noValues()<=getData().size(),     DataArrayView::setOffset(ValueType::size_type offset)
253                 "Error - offset incompatible with this view.");     {
254          EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
255      EsysAssert(getRank()==region.size(),        if (checkOffset(offset)) {
256                 "Error - slice not same rank as this view.");           m_offset=offset;
257          } else {
258      EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),           throw DataException("Error - invalid offset specified.");
259                 "Error - slice shape not compatible shape for other view.");        }
260       }
     //  
     // copy the values in the other view into the specified region of this view  
   
     // allow for case where other view is a scalar  
     if (other.getRank()==0) {  
   
         // 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++;  
           break;  
         case 1:  
           for (int i=region[0].first;i<region[0].second;i++) {  
         (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];  
         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++;  
         }  
           }  
           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++;  
           }  
         }  
           }  
           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++;  
             }  
           }  
         }  
           }  
           break;  
         default:  
           stringstream mess;  
           mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();  
           throw DataException(mess.str());  
         }  
   
     } 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());  
         }  
261    
262      }     void
263       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 {
269             throw DataException("Error - Cannot increment offset.");
270          }
271       }
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    //       return result;
314    //    }
315    
316    //    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 {
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 {
367    //          /* if key is a pair extract begin and end values */
368    //          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    //    }
409    
410       void
411       DataArrayView::copySlice(const DataArrayView& other,
412                                const RegionLoopRangeType& region)
413       {
414          //
415          // version of copySlice that uses the default offsets
416          copySlice(m_offset,other,other.m_offset,region);
417       }
418    
419       void
420       DataArrayView::copySlice(ValueType::size_type thisOffset,
421                                const DataArrayView& other,
422                                ValueType::size_type otherOffset,
423                                const RegionLoopRangeType& region)
424       {
425    
426          //
427          // Make sure views are not empty
428    
429          EsysAssert(!isEmpty(),
430                     "Error - this view is empty.");
431          EsysAssert(!other.isEmpty(),
432                     "Error - other view is empty.");
433    
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()==DataTypes::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.
459             here for completeness only, allows slicing of a scalar */
460             (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
461             numCopy++;
462             break;
463          case 1:
464             for (int i=region[0].first;i<region[0].second;i++) {
465                (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
466                numCopy++;
467             }
468             break;
469          case 2:
470             for (int j=region[1].first;j<region[1].second;j++) {
471                for (int i=region[0].first;i<region[0].second;i++) {
472                   (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
473                   numCopy++;
474                }
475             }
476             break;
477          case 3:
478             for (int k=region[2].first;k<region[2].second;k++) {
479                for (int j=region[1].first;j<region[1].second;j++) {
480                   for (int i=region[0].first;i<region[0].second;i++) {
481                      (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
482                      numCopy++;
483                   }
484                }
485             }
486             break;
487          case 4:
488             for (int l=region[3].first;l<region[3].second;l++) {
489                for (int k=region[2].first;k<region[2].second;k++) {
490                   for (int j=region[1].first;j<region[1].second;j++) {
491                      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)];
493                         numCopy++;
494                      }
495                   }
496                }
497             }
498             break;
499          default:
500             stringstream mess;
501             mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
502             throw DataException(mess.str());
503          }
504       }
505    
506  string     void
507  DataArrayView::toString(const string& suffix) const     DataArrayView::copySliceFrom(const DataArrayView& other,
508  {                                  const RegionLoopRangeType& region_loop_range)
509      EsysAssert(!isEmpty(),"Error - View is empty.");     {
510      stringstream temp;        //
511      string finalSuffix=suffix;        // version of copySliceFrom that uses the default offsets
512      if (suffix.length() > 0) {        copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
513        finalSuffix+=" ";     }
514      }  
515      switch (getRank()) {     void
516      case 0:     DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
517        temp << finalSuffix << (*this)();                                  const DataArrayView& other,
518        break;                                  ValueType::size_type otherOffset,
519      case 1:                                  const RegionLoopRangeType& region)
520        for (int i=0;i<getShape()[0];i++) {     {
521      temp << finalSuffix << "(" << i << ") " << (*this)(i);  
522      if (i!=(getShape()[0]-1)) {        //
523        temp << endl;        // Make sure views are not empty
524      }  
525        }        EsysAssert(!isEmpty(),
526        break;                   "Error - this view is empty.");
527      case 2:        EsysAssert(!other.isEmpty(),
528        for (int i=0;i<getShape()[0];i++) {                   "Error - other view is empty.");
529      for (int j=0;j<getShape()[1];j++) {  
530        temp << finalSuffix << "(" << i << "," << j << ") " << (*this)(i,j);        //
531        if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {        // Check this view is compatible with the region to be sliced,
532          temp << endl;        // and that the region to be sliced is compatible with the other view:
533        }  
534      }        EsysAssert(other.checkOffset(otherOffset),
535        }                   "Error - offset incompatible with other view.");
536        break;        EsysAssert(thisOffset+other.noValues()<=getData().size(),
537      case 3:                   "Error - offset incompatible with this view.");
538        for (int i=0;i<getShape()[0];i++) {  
539      for (int j=0;j<getShape()[1];j++) {        EsysAssert(getRank()==region.size(),
540        for (int k=0;k<getShape()[2];k++) {                   "Error - slice not same rank as this view.");
541          temp << finalSuffix << "(" << i << "," << j << "," << k << ") " << (*this)(i,j,k);  
542          if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {        EsysAssert(other.getRank()==0 || other.noValues()==DataTypes::noValues(getResultSliceShape(region)),
543            temp << endl;                   "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
547        }  
548        break;        // allow for case where other view is a scalar
549      case 4:        if (other.getRank()==0) {
550        for (int i=0;i<getShape()[0];i++) {  
551      for (int j=0;j<getShape()[1];j++) {           // the following loops cannot be parallelised due to the numCopy counter
552        for (int k=0;k<getShape()[2];k++) {           int numCopy=0;
553          for (int l=0;l<getShape()[3];l++) {  
554            temp << finalSuffix << "(" << i << "," << j << "," << k << "," << l << ") " << (*this)(i,j,k,l);           switch (region.size()) {
555            if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {           case 0:
556              temp << endl;              /* this case should never be encountered,
557            }              as python will never pass us an empty region.
558              }              here for completeness only, allows slicing of a scalar */
559        }              (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
560      }              numCopy++;
561        }              break;
562        break;           case 1:
563      default:              for (int i=region[0].first;i<region[0].second;i++) {
564        stringstream mess;                 (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
565        mess << "Error - (toString) Invalid rank: " << getRank();                 numCopy++;
566        throw DataException(mess.str());              }
567      }              break;
568      return temp.str();           case 2:
569  }              for (int j=region[1].first;j<region[1].second;j++) {
570                   for (int i=region[0].first;i<region[0].second;i++) {
571  string                    (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
572  DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)                    numCopy++;
573  {                 }
574      stringstream temp;              }
575      temp << "(";              break;
576      for (int i=0;i<shape.size();i++) {           case 3:
577        temp << shape[i];              for (int k=region[2].first;k<region[2].second;k++) {
578        if (i < shape.size()-1) {                 for (int j=region[1].first;j<region[1].second;j++) {
579      temp << ",";                    for (int i=region[0].first;i<region[0].second;i++) {
580        }                       (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
581      }                       numCopy++;
582      temp << ")";                    }
583      return temp.str();                 }
584  }              }
585                break;
586  void           case 4:
587  DataArrayView::matMult(const DataArrayView& left,              for (int l=region[3].first;l<region[3].second;l++) {
588                         const DataArrayView& right,                 for (int k=region[2].first;k<region[2].second;k++) {
589                         DataArrayView& result)                    for (int j=region[1].first;j<region[1].second;j++) {
590  {                       for (int i=region[0].first;i<region[0].second;i++) {
591                            (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
592                            numCopy++;
593                         }
594                      }
595                   }
596                }
597                break;
598             default:
599                stringstream mess;
600                mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
601                throw DataException(mess.str());
602             }
603    
604      if (left.getRank()==0 || right.getRank()==0) {        } else {
       stringstream temp;  
       temp << "Error - (matMult) Invalid for rank 0 objects.";  
       throw DataException(temp.str());  
     }  
605    
606      if (left.getShape()[left.getRank()-1] != right.getShape()[0]) {           // the following loops cannot be parallelised due to the numCopy counter
607        stringstream temp;           int numCopy=0;
       temp << "Error - (matMult) Dimension: " << left.getRank()  
        << ", size: " << left.getShape()[left.getRank()-1]  
        << " of LHS and dimension: 1, size: " << right.getShape()[0]  
        << " of RHS don't match.";  
       throw DataException(temp.str());  
     }  
608    
609      int outputRank = left.getRank()+right.getRank()-2;           switch (region.size()) {
610             case 0:
611                /* this case should never be encountered,
612                as python will never pass us an empty region.
613                here for completeness only, allows slicing of a scalar */
614                (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
615                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      if (outputRank < 0) {        }
       stringstream temp;  
       temp << "Error - (matMult) LHS and RHS cannot be multiplied "  
        << "as they have incompatable rank.";  
       throw DataException(temp.str());  
     }  
660    
661      if (outputRank != result.getRank()) {     }
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) Rank of result array is: "        string finalSuffix=suffix;
669         << result.getRank()        if (suffix.length() > 0) {
670         << " it must be: " << outputRank;           finalSuffix+=" ";
671        throw DataException(temp.str());        }
672      }        switch (getRank()) {
673          case 0:
674      for (int i=0; i<(left.getRank()-1); i++) {           temp << finalSuffix << (*this)();
675        if (left.getShape()[i] != result.getShape()[i]) {           break;
676      stringstream temp;        case 1:
677      temp << "Error - (matMult) Dimension: " << i           for (int i=0;i<getShape()[0];i++) {
678           << " of LHS and result array don't match.";              temp << finalSuffix << "(" << i << ") " << (*this)(i);
679      throw DataException(temp.str());              if (i!=(getShape()[0]-1)) {
680        }                 temp << endl;
681      }              }
682             }
683      for (int i=1; i<right.getRank(); i++) {           break;
684        if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {        case 2:
685      stringstream temp;           for (int i=0;i<getShape()[0];i++) {
686      temp << "Error - (matMult) Dimension: " << i              for (int j=0;j<getShape()[1];j++) {
687           << ", size: " << right.getShape()[i]                 temp << finalSuffix << "(" << i << "," << j << ") " << (*this)(i,j);
688           << " of RHS and dimension: " << i+left.getRank()-1                 if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
689           << ", size: " << result.getShape()[i+left.getRank()-1]                    temp << endl;
690           << " of result array don't match.";                 }
691      throw DataException(temp.str());              }
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          for (int i=1; i<right.getRank(); i++) {
792             if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
793                stringstream temp;
794                temp << "Error - (matMult) Dimension: " << i
795                     << ", size: " << right.getShape()[i]
796                     << " of RHS and dimension: " << i+left.getRank()-1
797                     << ", size: " << result.getShape()[i+left.getRank()-1]
798                     << " of result array don't match.";
799                throw DataException(temp.str());
800             }
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.1734

  ViewVC Help
Powered by ViewVC 1.1.26