escript  Revision_
DataVector.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16 
17 #ifndef __ESCRIPT_DATAVECTOR_H__
18 #define __ESCRIPT_DATAVECTOR_H__
19 
20 #include "system_dep.h"
21 #include "DataTypes.h"
22 #include "Assert.h"
23 #include "DataVectorAlt.h"
24 #include "DataVectorTaipan.h"
25 
26 // ensure that nobody else tries to instantiate the complex version
28 
29 
30 namespace escript {
31 
32 // Functions in DataTypes:: which manipulate DataVectors
33 namespace DataTypes
34 {
35 
36  // This is the main version we had
37  //typedef DataVectorTaipan DataVector;
40 
56  void
57  pointToStream(std::ostream& os, const RealVectorType::ElementType* data,const ShapeType& shape, int offset, bool needsep=true, const std::string& sep=",");
58 
74  void
75  pointToStream(std::ostream& os, const CplxVectorType::ElementType* data,const ShapeType& shape, int offset, bool needsep=true, const std::string& sep=",");
76 
85  std::string
86  pointToString(const RealVectorType& data,const ShapeType& shape, int offset, const std::string& prefix);
87 
88 
89  std::string
90  pointToString(const CplxVectorType& data,const ShapeType& shape, int offset, const std::string& prefix);
91 
101  void copyPoint(RealVectorType& dest, vec_size_type doffset, vec_size_type nvals, const RealVectorType& src, vec_size_type soffset);
102 
112  void copyPoint(CplxVectorType& dest, vec_size_type doffset, vec_size_type nvals, const CplxVectorType& src, vec_size_type soffset);
113 
119  void fillComplexFromReal(const RealVectorType& r, CplxVectorType& c);
120 
135  template <class VEC>
137  void
138  copySlice(VEC& left,
139  const ShapeType& leftShape,
140  typename VEC::size_type leftOffset,
141  const VEC& other,
142  const ShapeType& otherShape,
143  typename VEC::size_type otherOffset,
144  const RegionLoopRangeType& region)
145  {
146  //
147  // Make sure vectors are not empty
148 
149  ESYS_ASSERT(!left.size()==0, "left data is empty.");
150  ESYS_ASSERT(!other.size()==0, "other data is empty.");
151 
152  //
153  // Check the vector to be sliced from is compatible with the region to be sliced,
154  // and that the region to be sliced is compatible with this vector:
155  ESYS_ASSERT(checkOffset(leftOffset,left.size(),noValues(leftShape)),
156  "offset incompatible with this vector.");
157  ESYS_ASSERT(otherOffset+noValues(leftShape)<=other.size(),
158  "offset incompatible with other vector.");
159 
160  ESYS_ASSERT(getRank(otherShape)==region.size(),
161  "slice not same rank as vector to be sliced from.");
162 
163  ESYS_ASSERT(noValues(leftShape)==noValues(getResultSliceShape(region)),
164  "slice shape not compatible shape for this vector.");
165 
166  //
167  // copy the values in the specified region of the other vector into this vector
168 
169  // the following loops cannot be parallelised due to the numCopy counter
170  int numCopy=0;
171 
172  switch (region.size()) {
173  case 0:
174  /* this case should never be encountered,
175  as python will never pass us an empty region.
176  here for completeness only, allows slicing of a scalar */
177 // (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
178 
179  left[leftOffset+numCopy]=other[otherOffset];
180  numCopy++;
181  break;
182  case 1:
183  for (int i=region[0].first;i<region[0].second;i++) {
184  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
185  numCopy++;
186  }
187  break;
188  case 2:
189  for (int j=region[1].first;j<region[1].second;j++) {
190  for (int i=region[0].first;i<region[0].second;i++) {
191 /* (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
192  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
193  numCopy++;
194  }
195  }
196  break;
197  case 3:
198  for (int k=region[2].first;k<region[2].second;k++) {
199  for (int j=region[1].first;j<region[1].second;j++) {
200  for (int i=region[0].first;i<region[0].second;i++) {
201 // (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
202  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
203  numCopy++;
204  }
205  }
206  }
207  break;
208  case 4:
209  for (int l=region[3].first;l<region[3].second;l++) {
210  for (int k=region[2].first;k<region[2].second;k++) {
211  for (int j=region[1].first;j<region[1].second;j++) {
212  for (int i=region[0].first;i<region[0].second;i++) {
213 /* (*m_data)[leftOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
214  left[leftOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
215  numCopy++;
216  }
217  }
218  }
219  }
220  break;
221  default:
222  std::stringstream mess;
223  mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
224  throw DataException(mess.str());
225  }
226  }
227 
242  template<typename VEC>
244  void
245  copySliceFrom(VEC& left,
246  const ShapeType& leftShape,
247  typename VEC::size_type leftOffset,
248  const VEC& other,
249  const ShapeType& otherShape,
250  typename VEC::size_type otherOffset,
251  const RegionLoopRangeType& region)
252  {
253  //
254  // Make sure vectors are not empty
255 
256  ESYS_ASSERT(left.size()!=0, "this vector is empty.");
257  ESYS_ASSERT(other.size()!=0, "other vector is empty.");
258 
259  //
260  // Check this vector is compatible with the region to be sliced,
261  // and that the region to be sliced is compatible with the other vector:
262 
263  ESYS_ASSERT(checkOffset(otherOffset,other.size(),noValues(otherShape)),
264  "offset incompatible with other vector.");
265  ESYS_ASSERT(leftOffset+noValues(otherShape)<=left.size(),
266  "offset incompatible with this vector.");
267 
268  ESYS_ASSERT(getRank(leftShape)==region.size(),
269  "slice not same rank as this vector.");
270 
271  ESYS_ASSERT(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
272  "slice shape not compatible shape for other vector.");
273 
274  //
275  // copy the values in the other vector into the specified region of this vector
276 
277  // allow for case where other vector is a scalar
278  if (getRank(otherShape)==0) {
279 
280  // the following loops cannot be parallelised due to the numCopy counter
281  int numCopy=0;
282 
283  switch (region.size()) {
284  case 0:
285  /* this case should never be encountered,
286  as python will never pass us an empty region.
287  here for completeness only, allows slicing of a scalar */
288  //(*m_data)[leftOffset+relIndex()]=(*other.m_data)[otherOffset];
289  left[leftOffset]=other[otherOffset];
290  numCopy++;
291  break;
292  case 1:
293  for (int i=region[0].first;i<region[0].second;i++) {
294  left[leftOffset+getRelIndex(leftShape,i)]=other[otherOffset];
295  numCopy++;
296  }
297  break;
298  case 2:
299  for (int j=region[1].first;j<region[1].second;j++) {
300  for (int i=region[0].first;i<region[0].second;i++) {
301  left[leftOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
302  numCopy++;
303  }
304  }
305  break;
306  case 3:
307  for (int k=region[2].first;k<region[2].second;k++) {
308  for (int j=region[1].first;j<region[1].second;j++) {
309  for (int i=region[0].first;i<region[0].second;i++) {
310  left[leftOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset];
311  numCopy++;
312  }
313  }
314  }
315  break;
316  case 4:
317  for (int l=region[3].first;l<region[3].second;l++) {
318  for (int k=region[2].first;k<region[2].second;k++) {
319  for (int j=region[1].first;j<region[1].second;j++) {
320  for (int i=region[0].first;i<region[0].second;i++) {
321  left[leftOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
322  numCopy++;
323  }
324  }
325  }
326  }
327  break;
328  default:
329  std::stringstream mess;
330  mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
331  throw DataException(mess.str());
332  }
333 
334  } else {
335 
336  // the following loops cannot be parallelised due to the numCopy counter
337  int numCopy=0;
338 
339  switch (region.size()) {
340  case 0:
341  /* this case should never be encountered,
342  as python will never pass us an empty region.
343  here for completeness only, allows slicing of a scalar */
344  //(*m_data)[leftOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
345  left[leftOffset]=other[otherOffset+numCopy];
346  numCopy++;
347  break;
348  case 1:
349  for (int i=region[0].first;i<region[0].second;i++) {
350  left[leftOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
351  numCopy++;
352  }
353  break;
354  case 2:
355  for (int j=region[1].first;j<region[1].second;j++) {
356  for (int i=region[0].first;i<region[0].second;i++) {
357  left[leftOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
358  numCopy++;
359  }
360  }
361  break;
362  case 3:
363  for (int k=region[2].first;k<region[2].second;k++) {
364  for (int j=region[1].first;j<region[1].second;j++) {
365  for (int i=region[0].first;i<region[0].second;i++) {
366  left[leftOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset+numCopy];
367  numCopy++;
368  }
369  }
370  }
371  break;
372  case 4:
373  for (int l=region[3].first;l<region[3].second;l++) {
374  for (int k=region[2].first;k<region[2].second;k++) {
375  for (int j=region[1].first;j<region[1].second;j++) {
376  for (int i=region[0].first;i<region[0].second;i++) {
377  left[leftOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
378  numCopy++;
379  }
380  }
381  }
382  }
383  break;
384  default:
385  std::stringstream mess;
386  mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
387  throw DataException(mess.str());
388  }
389 
390  }
391 
392  }
393 }
394 
395 } // end of namespace
396 
397 #endif // __ESCRIPT_DATAVECTOR_H__
398 
bool checkOffset(vec_size_type offset, int size, int noval)
Definition: DataTypes.h:323
Definition: AbstractContinuousDomain.cpp:22
void copyPoint(RealVectorType &dest, vec_size_type doffset, vec_size_type nvals, const RealVectorType &src, vec_size_type soffset)
Copy a point from one vector to another. Note: This version does not check to see if shapes are the s...
std::string pointToString(const RealVectorType &data, const ShapeType &shape, int offset, const std::string &prefix)
Display a single value (with the specified shape) from the data.
vec_size_type getRelIndex(const DataTypes::ShapeType &shape, vec_size_type i)
Compute the offset (in 1D vector) of a given subscript with a shape.
Definition: DataTypes.h:233
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:42
std::vector< std::pair< int, int > > RegionLoopRangeType
Definition: DataTypes.h:44
DataTypes::ShapeType getResultSliceShape(const RegionType &region)
Determine the shape of the specified slice region.
Definition: DataTypes.cpp:172
void copySliceFrom(VEC &left, const ShapeType &leftShape, typename VEC::size_type leftOffset, const VEC &other, const ShapeType &otherShape, typename VEC::size_type otherOffset, const RegionLoopRangeType &region)
Copy data into a slice specified by the given region and offset in the left vector from the other vec...
Definition: DataVector.h:245
Definition: DataVectorAlt.h:36
void pointToStream(std::ostream &os, const RealVectorType::ElementType *data, const ShapeType &shape, int offset, bool needsep=true, const std::string &sep=",")
Display a single value (with the specified shape) from the data.
void copySlice(VEC &left, const ShapeType &leftShape, typename VEC::size_type leftOffset, const VEC &other, const ShapeType &otherShape, typename VEC::size_type otherOffset, const RegionLoopRangeType &region)
Copy a data slice specified by the given region and offset from the "other" vector into the "left" ve...
Definition: DataVector.h:138
T ElementType
Definition: DataVectorAlt.h:42
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:90
Definition: DataException.h:26
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:29
escript::DataTypes::DataVectorAlt< real_t > RealVectorType
Vector to store underlying data.
Definition: DataVector.h:38
escript::DataTypes::DataVectorAlt< cplx_t > CplxVectorType
Definition: DataVector.h:39
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:218
#define ESYS_ASSERT(a, b)
EsysAssert is a MACRO that will throw an exception if the boolean condition specified is false...
Definition: Assert.h:78
long vec_size_type
Definition: DataTypes.h:47
void fillComplexFromReal(const RealVectorType &r, CplxVectorType &c)
copy data from a real vector to a complex vector The complex vector will be resized as needed and any...