MADNESS  0.10.1
tensoriter.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31 
32  $Id$
33 */
34 
35 
36 #ifndef MADNESS_TENSOR_TENSORITER_H__INCLUDED
37 #define MADNESS_TENSOR_TENSORITER_H__INCLUDED
38 
39 /// \file tensoriter.h
40 /// \brief Declares TensorIterator
41 
42 #include <madness/tensor/tensor.h>
44 
45 #include <iostream>
46 #include <algorithm>
47 #include <complex>
48 
49 #include <cmath>
50 
51 
52 namespace madness {
53 
54  template <class T> class Tensor;
55  template <class T> class SliceTensor;
56 
57  static const long default_jdim = 5551212; // Never a valid dimension num.
58 
59  /// Optimized iterator for tensors supporting unary, binary and ternary operations.
60  /// \ingroup tensor
61  template <class T, class Q = T, class R = T> class TensorIterator {
65  public:
66  T* _p0;
67  Q* _p1;
68  R* _p2;
69  long ndim;
70  long dimj;
71  long _s0;
72  long _s1;
73  long _s2;
79 
80  TensorIterator(const Tensor<T>* t0, const Tensor<Q>* t1=0, const Tensor<R>* t2=0,
81  long iterlevel=0,
82  bool optimize=true, bool fusedim=true,
83  long jdim=default_jdim);
84 
86 
87  inline bool operator == (const TensorIterator<T,Q,R>& a) const {
88  return _p0==a._p0;
89  }
90 
91  inline bool operator != (const TensorIterator<T,Q,R>& a) const {
92  return _p0!=a._p0;
93  }
94 
96  return this;
97  };
98 
99  inline T& operator*() const {
100  return *_p0;
101  };
102 
103  void reset();
104 
105  void reuse(const Tensor<T>* t0, const Tensor<Q>* t1=0, const Tensor<R>* t2=0);
106  };
107 
108  /// Constructor for general iterator to compose operations over up to three tensors
109 
110  /// \ingroup tensor
111  /// Macros have been defined in \c tensor_macros.h to take the pain out of
112  /// using iterators. The options have the following effects
113  /// - \c optimize reorders dimensions for optimal strides (only
114  /// applies if iterlevel=1). If jdim==default_jdim, all dimensions
115  /// are reordered for optimal stride. If jdim is not the default
116  /// value, then dimension jdim is excluded from the set of
117  /// dimensions being optimized.
118  ///
119  /// - \c fusedim concatenates contiguous dimensions into the inner
120  /// loop (only applies if iterlevel=1 and jdim=default_jdim).
121  ///
122  /// - \c iterlevel can have two values 0=elementwise, 1=vectorwise. Elementwise implies
123  /// that the iterator returns successive elements, and explicitly
124  /// in the expected order for the tensor (i.e., with the last index
125  /// varying fastest). Vectorwise implies that the user is responsible
126  /// for iterating over dimension jdim.
127  ///
128  /// - jdim --- if iterlevel=1, then jdim determines which dimension is
129  /// left for iteration with an explicit for loop. Negative values
130  /// implies jdim+=ndim following the convention of Slice (and
131  /// Python). If (optimize) the default is the fastest varying
132  /// dimension. If (!optimize) the default is the last dimension
133  /// (jdim=-1). If (fusedim) contiguous dimensions are fused into
134  /// this inner loop. Specifying a non-default value for jdim
135  /// disables fusedim and restricts optimization to reordering only
136  /// the exterior loops (so that the loop the user is iterating over
137  /// corresponds exactly to those in dimension jdim).
138  ///
139  /// \par During iteration:
140  ///
141  /// - \c ind[] will contain the current iteration indices .. BUT
142  /// if \c optimize=true , they will not necessarily be in the
143  /// order of those of the tensor. if fusedim is true, then
144  /// there may be fewer dimensions than the input tensor.
145  ///
146  /// - \c _p0, \c _p1, \c _p2 will point to the current elements of \c t0,t1,t2
147  /// noting that if \c iterlevel>0, the user must provide additional
148  /// \c for loops to iterate over the additional dimensions
149  ///
150  /// - \c stride0[], \c stride1[], \c stride2[] will contain the strides for each
151  /// (possibly reordered) dimension. If \c iterlevel=1, \c _s0,_s1,s2
152  /// contain the stride info for the dimension that the user is
153  /// responsible for iterating over.
154  ///
155  /// - \c dimj -> the size of the j'th dimension
156  template <class T, class Q, class R>
158  const Tensor<Q>* t1,
159  const Tensor<R>* t2,
160  long iterlevel,
161  bool optimize,
162  bool fusedim,
163  long jdim) {
164 
165 
166  if (!t0) {
167  // Used to indicate end of iteration.
168  _p0 = 0;
169  return;
170  }
171 
172  //std::printf("t0=%p t1=%p t2=%p optimize=%d fusedim=%d iterlevel=%ld jdim=%ld\n",
173  //t0,t1,t2,optimize,fusedim,iterlevel,jdim);
174 
175  TENSOR_ASSERT(iterlevel==0 || iterlevel==1,"invalid iteration level",iterlevel,t0);
176 
177  // First copy basic info over
178  ndim = t0->ndim();
179  _p0_save = _p0 = const_cast<T*>(t0->ptr()); // CONSTNESS VIOLATION
180  for (int i=0; i<ndim; ++i) {
181  dim[i] = t0->dim(i);
182  stride0[i] = t0->stride(i);
183  }
184  if (t1) {
185  TENSOR_ASSERT(t0->conforms(*t1),"first and second tensors do not conform",
186  0, t0);
187  _p1_save = _p1 = const_cast<Q*>(t1->ptr()); // CONSTNESS VIOLATION
188  for (int i=0; i<ndim; ++i) stride1[i] = t1->stride(i);
189  }
190  else {
191  _p1_save = _p1 = 0;
192  }
193  if (t2) {
194  TENSOR_ASSERT(t0->conforms(*t2),"first and third tensors do not conform",
195  0, t0);
196  _p2_save = _p2 = const_cast<R*>(t2->ptr()); // CONSTNESS VIOLATION
197  for (int i=0; i<ndim; ++i) stride2[i] = t2->stride(i);
198  }
199  else {
200  _p2_save = _p2 = 0;
201  }
202 
203  if (iterlevel == 0) {
204  // Iteration will include all dimensions
205  fusedim = false;
206  jdim = ndim;
207  dimj = 0;
208  _s0 = 0;
209  _s1 = 0;
210  _s2 = 0;
211  }
212  else if (iterlevel == 1) {
213  // Apply -ve indexing convention for dimensions
214  if (jdim < 0) jdim += ndim;
215 
216  // If permissible optimize the order of dimensions excluding
217  // any non-default value of jdim.
218  if (optimize) {
219  for (int i=0; i<ndim; ++i) {
220  if (i != jdim) {
221  for (int j=i; j<ndim; ++j) {
222  if (j != jdim) {
223  if (std::abs(stride0[i]) < std::abs(stride0[j])) {
224  std::swap(stride0[i],stride0[j]);
225  if (t1) std::swap(stride1[i],stride1[j]);
226  if (t2) std::swap(stride2[i],stride2[j]);
227  std::swap(dim[i],dim[j]);
228  }
229  }
230  }
231  }
232  //std::cout << "stride0[" << i << "]=" << stride0[i] << std::endl;
233  }
234  }
235 
236  // Iterations will exclude dimension jdim, default is last one
237  if (jdim == default_jdim) {
238  jdim = ndim-1;
239  }
240  else {
241  fusedim = false;
242  }
243  TENSOR_ASSERT(jdim>=0 && jdim < ndim, "invalid index for external iteration",
244  jdim, t0);
245  ndim--;
246 
247  // Stride and dimension info for the excluded dimension
248  _s0 = stride0[jdim];
249  if (t1) {
250  _s1 = stride1[jdim];
251  }
252  else {
253  _s1 = 0;
254  }
255  if (t2) {
256  _s2 = stride2[jdim];
257  }
258  else {
259  _s2 = 0;
260  }
261  dimj = dim[jdim];
262 
263  // Collapse stride and dimension info for remaining dimensions
264  for (int i=jdim+1; i<=ndim; ++i) { // GCC COMPILER WARNING ... but need <= ... why?
265  dim[i-1] = dim[i];
266  stride0[i-1] = stride0[i];
267  }
268  if (t1) {
269  for (int i=jdim+1; i<=ndim; ++i) stride1[i-1] = stride1[i]; // GCC COMPILER WARNING ... but need <= ... why?
270  }
271  if (t2) {
272  for (int i=jdim+1; i<=ndim; ++i) stride2[i-1] = stride2[i]; // GCC COMPILER WARNING ... but need <= ... why?
273  }
274 
275  if (fusedim) { // Only if jdim=default_jdim && iterlevel=1
276  if (t2) {
277  for (int i=ndim-1; i>=0; --i) {
278  if (_s0*dimj == stride0[i] &&
279  _s1*dimj == stride1[i] &&
280  _s2*dimj == stride2[i]) {
281  dimj *= dim[i];
282  ndim--;
283  }
284  else {
285  break;
286  }
287  }
288  }
289  else if (t1) {
290  for (int i=ndim-1; i>=0; --i) {
291  if (_s0*dimj == stride0[i] &&
292  _s1*dimj == stride1[i]) {
293  dimj *= dim[i];
294  ndim--;
295  }
296  else {
297  break;
298  }
299  }
300  }
301  else {
302  for (int i=ndim-1; i>=0; --i) {
303  if (_s0*dimj == stride0[i]) {
304  dimj *= dim[i];
305  ndim--;
306  }
307  else {
308  break;
309  }
310  }
311  }
312  }
313  }
314 
315  // Initialize indices for the counter Use TENSOR_MAXDIM so reference to
316  // optimized-away dimensions is vaguely meaningful.
317  for (int i=0; i<TENSOR_MAXDIM; ++i) ind[i] = 0;
318 
319  // std::printf("ndim=%ld dimj=%ld _s0=%ld _s1=%ld _s2=%ld _p0=%p _p1=%p _p2=%p\n",
320  // ndim,dimj,_s0,_s1,_s2,_p0,_p1,_p2);
321  // for (int i=0; i<ndim; ++i) {
322  // std::printf(" %d dim=%ld stride0=%ld stride1=%ld stride2=%ld\n",
323  // i,dim[i],stride0[i],stride1[i],stride2[i]);
324  // }
325  }
326 
327  template <class T, class Q, class R>
329  long d = ndim-1;
330  if (d<0 || _p0==0) {
331  _p0 = 0;
332  return *this;
333  }
334  while (ind[d] >= (dim[d] - 1)) {
335  _p0 -= ind[d] * stride0[d];
336  if (_p1) _p1 -= ind[d] * stride1[d];
337  if (_p2) _p2 -= ind[d] * stride2[d];
338  ind[d] = 0;
339  d--;
340  if (d < 0) {
341  _p0 = 0;
342  return *this;
343  }
344  }
345  _p0 += stride0[d];
346  if (_p1) _p1 += stride1[d];
347  if (_p2) _p2 += stride2[d];
348  ++(ind[d]);
349  return *this;
350  }
351 
352  /// Reset the iterator back to the start ...
353  template <class T, class Q, class R>
355  _p0 = _p0_save;
356  _p1 = _p1_save;
357  _p2 = _p2_save;
358  for (int i=0; i<TENSOR_MAXDIM; ++i) ind[i] = 0;
359  }
360 
361 
362  /// Reuse this iterator to iterate over other Tensors
363 
364  /// The point of this method is to optimize away the construction of a
365  /// TensorIterator when applying the same operation repeatedly to
366  /// multiple tensors with identical shapes & strides. We trust the
367  /// caller to ensure all shapes & strides match between the new
368  /// tensors and the ones used in the original constructor.
369  template <class T, class Q, class R>
371  const Tensor<Q> *t1,
372  const Tensor<R> *t2) {
373  _p0 = _p0_save = t0->ptr();
374  if (t1) _p1 = _p1_save = t1->ptr();
375  if (t2) _p2 = _p2_save = t2->ptr();
376  for (int i=0; i<TENSOR_MAXDIM; ++i) ind[i] = 0;
377  }
378 
379 
380 }
381 #endif // MADNESS_TENSOR_TENSORITER_H__INCLUDED
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
long stride(int i) const
Returns the stride associated with dimension i.
Definition: basetensor.h:150
long ndim() const
Returns the number of dimensions in the tensor.
Definition: basetensor.h:144
Definition: tensoriter.h:61
R * _p2_save
Definition: tensoriter.h:64
long stride0[TENSOR_MAXDIM]
Definition: tensoriter.h:76
bool operator==(const TensorIterator< T, Q, R > &a) const
Definition: tensoriter.h:87
T & operator*() const
Definition: tensoriter.h:99
long dim[TENSOR_MAXDIM]
Definition: tensoriter.h:74
long stride1[TENSOR_MAXDIM]
Definition: tensoriter.h:77
T * _p0_save
Definition: tensoriter.h:62
void reuse(const Tensor< T > *t0, const Tensor< Q > *t1=0, const Tensor< R > *t2=0)
Reuse this iterator to iterate over other Tensors.
Definition: tensoriter.h:370
long dimj
Definition: tensoriter.h:70
long _s1
Definition: tensoriter.h:72
bool operator!=(const TensorIterator< T, Q, R > &a) const
Definition: tensoriter.h:91
long _s0
Definition: tensoriter.h:71
Q * _p1
Definition: tensoriter.h:67
long ndim
Definition: tensoriter.h:69
TensorIterator< T, Q, R > & operator++()
Definition: tensoriter.h:328
void reset()
Reset the iterator back to the start ...
Definition: tensoriter.h:354
long ind[TENSOR_MAXDIM]
Definition: tensoriter.h:75
Q * _p1_save
Definition: tensoriter.h:63
TensorIterator< T, Q, R > * operator->()
Definition: tensoriter.h:95
R * _p2
Definition: tensoriter.h:68
T * _p0
Definition: tensoriter.h:66
long stride2[TENSOR_MAXDIM]
Definition: tensoriter.h:78
long _s2
Definition: tensoriter.h:73
A tensor is a multidimension array.
Definition: tensor.h:317
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
bool conforms(const Tensor< Q > &t) const
Test if *this and t conform.
Definition: tensor.h:1657
static const double R
Definition: csqrt.cc:46
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
TensorIterator(const Tensor< T > *t0, const Tensor< Q > *t1=0, const Tensor< R > *t2=0, long iterlevel=0, bool optimize=true, bool fusedim=true, long jdim=default_jdim)
Constructor for general iterator to compose operations over up to three tensors.
Definition: tensoriter.h:157
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
static const long default_jdim
Definition: tensoriter.h:57
void swap(Vector< T, N > &l, Vector< T, N > &r)
Swap the contents of two Vectors.
Definition: vector.h:497
static long abs(long a)
Definition: tensor.h:218
static const double a
Definition: nonlinschro.cc:118
double Q(double a)
Definition: relops.cc:20
Defines and implements most of Tensor.
#define TENSOR_MAXDIM
Definition: tensor_macros.h:194
Declares and implements TensorException.
#define TENSOR_ASSERT(condition, msg, value, t)
Definition: tensorexcept.h:130
void d()
Definition: test_sig.cc:79