MADNESS  0.10.1
tensor.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 #ifndef MADNESS_TENSOR_TENSOR_H__INCLUDED
33 #define MADNESS_TENSOR_TENSOR_H__INCLUDED
34 
35 #include <madness/madness_config.h>
36 #include <madness/misc/ran.h>
37 #include <madness/world/posixmem.h>
38 
39 #include <memory>
40 #include <complex>
41 #include <vector>
42 #include <array>
43 #include <cmath>
44 #include <cstdlib>
45 #include <cstddef>
46 
47 #include <madness/world/archive.h>
48 // #include <madness/world/print.h>
49 //
50 // typedef std::complex<float> float_complex;
51 // typedef std::complex<double> double_complex;
52 //
53 // // These probably have to be included in this order
54 // #include <madness/tensor/tensor_macros.h>
55 // #include <madness/tensor/type_data.h>
56 // #include <madness/tensor/slice.h>
57 // #include <madness/tensor/vector_factory.h>
59 #include <madness/tensor/aligned.h>
60 #include <madness/tensor/mxm.h>
63 
64 #ifdef ENABLE_GENTENSOR
65 #define HAVE_GENTENSOR 1
66 #else
67 #define HAVE_GENTENSOR 0
68 #endif
69 
70 
71 /*!
72  \file tensor.h
73  \brief Defines and implements most of Tensor
74  \ingroup tensor
75  \addtogroup tensor
76 
77  \par Introduction
78 
79  A tensor is a multi-dimensional array and does not incorporate any concepts
80  of covariance and contravariance.
81 
82  When a new tensor is created, the underlying data is also allocated.
83  E.g.,
84  \code
85  Tensor<double> a(3,4,5)
86  \endcode
87  creates a new 3-dimensional tensor and allocates a contiguous
88  block of 60 doubles which are initialized to zero. The dimensions
89  (numbered from the left starting at 0) are in C or row-major
90  order. Thus, for the tensor \c a , the stride between successive
91  elements of the right-most dimension is 1. For the middle
92  dimension it is 5. For the left-most dimension it is 20. Thus,
93  the loops
94  \code
95  for (i=0; i<3; ++i)
96  for (j=0; j<4; ++j)
97  for (k=0; k<5; ++k)
98  a(i,j,k) = ...
99  \endcode
100  will go sequentially (and thus efficiently) through memory.
101  If the dimensions have been reordered (e.g., with \c swapdim()
102  or \c map() ), or if the tensor is actually a slice of another
103  tensor, then the layout in memory may be more complex and
104  may not reflect a contiguous block of memory.
105 
106  Multiple tensors may be used to provide multiple identical or
107  distinct views of the same data. E.g., in the following
108  \code
109  Tensor<double> a(2,3); // A new tensor initialized to zero
110  Tensor<double> b = a;
111  \endcode
112  \c a and \c b provide identical views of the same data, thus
113  \code
114  b(1,2) = 99;
115  cout << a(1,2) << endl; // Outputs 99
116  cout << b(1,2) << endl; // Outputs 99
117  \endcode
118 
119  \par Shallow copy and assignment
120 
121  It is important to appreciate that the views and the data are
122  quite independent. In particular, the default copy constructor
123  and assignment operations only copy the tensor (the view) and not
124  the data --- <em> i.e., the copy constructor and assigment operations
125  only take shallow copies</em>. This is for both consistency and
126  efficiency. Thus, assigning one tensor to another generates another
127  view of the same data, replacing any previous view and not moving
128  or copying any of the data.
129  E.g.,
130  \code
131  Tensor<double> a(2,3); // A new tensor initialized to zero
132  Tensor<double> c(3,3,3); // Another new tensor
133  Tensor<double> b = a; // b is a view of the same data as a
134  a = c; // a is now a view of c's data
135  b = c // b is now also a view of c's data and the
136  // data allocated originally for a is freed
137  \endcode
138  The above example also illustrates how reference counting is used
139  to keep track of the underlying data. Once there are no views
140  of the data, it is automatically freed.
141 
142  There are only two ways to actually copy the underlying data. A
143  new, complete, and contiguous copy of a tensor and its data may be
144  generated with the \c copy() function. Or, to copy data from one tensor
145  into the data viewed by another tensor, you must use a Slice.
146 
147  \par Indexing
148 
149  One dimensional tensors (i.e., vectors) may be indexed using
150  either square brackets (e.g., \c v[i] ) or round brackets (e.g.,
151  \c v(i) ). All higher-dimensional tensors must use only round
152  brackets (e.g., \c t(i,j,k) ). This is due to C++'s restriction
153  that the indexing operator (\c [] ) can only have one argument.
154  The indexing operation should generate efficient code.
155 
156  For the sake of efficiency, no bounds checking is performed by
157  default by most single element indexing operations. Checking can
158  be enabled at compile time by defining \c -DTENSOR_BOUNDS_CHECKING for
159  application files including \c tensor.h. The MADNESS configure script
160  has the option \c --enable-tensor-bound-checking to define the macro
161  in \c madness_config.h . The general indexing
162  operation that takes a \c std::vector<long> index and all slicing
163  operations always perform bounds checking. To make indexing with
164  checking a bit easier, a factory function has been provided for
165  vectors ... but note you need to explicitly use longs as the
166  index.
167  \code
168  Tensor<long> a(7,7,7);
169  a(3,4,5) += 1; // OK ... adds 1 to element (3,4,5)
170  a(3,4,9) += 1; // BAD ... undetected out-of-bounds access
171  a(vector_factory(3L,4L,9L)) += 1; // OK ... out-bounds access will
172  // be detected at runtime.
173  \endcode
174 
175  \par Slicing
176 
177  Slices generate sub-tensors --- i.e., views of patches of the
178  data. E.g., to refer to all but the first and last elements in
179  each dimension of a matrix use
180  \code
181  a(Slice(1,-2),Slice(1,-2))
182  \endcode
183  Or to view odd elements in each dimension
184  \code
185  a(Slice(0,-1,2),Slice(0,-1,2))
186  \endcode
187  A slice or patch of a
188  tensor behaves exactly like a tensor \em except for assignment.
189  When a slice is assigned to, the data is copied with the
190  requirement that the source and destinations agree in size and
191  shape (i.e., they conform). Thus, to copy the all of the data
192  from a to b,
193  \code
194  Tensor<double> a(3,4), b(3,4);
195  a = 1; // Set all elements of a to 1
196  b = 2; // Set all elements of b to 2
197  a(Slice(0,-1,1),Slice(0,-1,1)) = b; // Copy all data from b to a
198  a(_,_) = b(_,_); // Copy all data from b to a
199  a(___) = b(___); // Copy all data from b to a
200  a(Slice(1,2),Slice(1,2)) = b; // Error, do not conform
201  \endcode
202  Special slice values \c _ ,\c _reverse, and \c ___ have
203  been defined to refer to all elements in a dimension, all
204  elements in a dimension but reversed, and all elements in all
205  dimensions, respectively.
206 
207  \par Iteration and algorithms
208 
209  See tensor_macros.h for documentation on the easiest mechanisms for iteration over
210  elements of tensors and tips for optimization. See \c TensorIterator for
211  the most general form of iteration.
212 */
213 
214 
215 #ifndef HAVE_STD_ABS_LONG
216 #ifndef HAVE_STD_LABS
217 namespace std {
218  static long abs(long a) {
219  return a>=0 ? a : -a;
220  }
221 }
222 #else
223 namespace std {
224  static long abs(long a) {
225  return std::labs(a);
226  }
227 }
228 #endif
229 #endif
230 
231 
232 namespace madness {
233 #define IS_ODD(n) ((n)&0x1)
234 #define IS_UNALIGNED(p) (((unsigned long)(p))&0x7)
235 
236 
237  /// For real types return value, for complex return conjugate
238  template <typename Q, bool iscomplex>
240  static Q op(const Q& coeff) {
241  return coeff;
242  }
243  };
244 
245  /// For real types return value, for complex return conjugate
246  template <typename Q>
247  struct conditional_conj_struct<Q,true> {
248  static Q op(const Q& coeff) {
249  return conj(coeff);
250  }
251  };
252 
253  /// For real types return value, for complex return conjugate
254  template <typename Q>
255  Q conditional_conj(const Q& coeff) {
257  }
258 
259  namespace detail {
260  template <typename T> T mynorm(T t) {
261  return t*t;
262  }
263  template <typename T=int> double mynorm(int t) {
264  return double(t)*double(t);
265  }
266 
267  template <typename T> T mynorm(std::complex<T> t) {
268  return std::norm(t);
269  }
270  }
271 
272  template <class T> class SliceTensor;
273 
274 
275 
276  //#define TENSOR_USE_SHARED_ALIGNED_ARRAY
277 #ifdef TENSOR_USE_SHARED_ALIGNED_ARRAY
278 #define TENSOR_SHARED_PTR detail::SharedAlignedArray
279  // this code has been tested and seems to work correctly on all
280  // tests and moldft ... however initial testing indicated a hard
281  // to measure improvement in thread scaling (from 20 to 60 threads
282  // on cn-mem) and no change in the modest thread count (20)
283  // execution time with tbballoc. hence we are not presently using
284  // it but will test more in the future with different allocators
285  namespace detail {
286  // Minimal ref-counted array with data+counter in one alloc
287  template <typename T> class SharedAlignedArray {
288  T* volatile p;
289  AtomicInt* cnt;
290  void dec() {if (p && ((*cnt)-- == 1)) {free(p); p = 0;}}
291  void inc() {if (p) (*cnt)++;}
292  public:
293  SharedAlignedArray() : p(0), cnt(0) {}
294  T* allocate(std::size_t size, unsigned int alignment) {
295  std::size_t offset = (size*sizeof(T)-1)/sizeof(AtomicInt) + 1; // Where the counter will be
296  std::size_t nbyte = (offset+1)*sizeof(AtomicInt);
297  if (posix_memalign((void **) &p, alignment, nbyte)) throw 1;
298  cnt = (AtomicInt*)(p) + offset;
299  *cnt = 1;
300  return p;
301  }
302  SharedAlignedArray<T>& operator=(const SharedAlignedArray<T>& other) {
303  if (this != &other) {dec(); p = other.p; cnt = other.cnt; inc();}
304  return *this;
305  }
306  void reset() {dec(); p = 0;}
307  ~SharedAlignedArray() {dec();}
308  };
309  }
310 #else
311 #define TENSOR_SHARED_PTR std::shared_ptr
312 #endif
313 
314  /// A tensor is a multidimension array
315 
316  /// \ingroup tensor
317  template <class T> class Tensor : public BaseTensor {
318  template <class U> friend class SliceTensor;
319 
320  protected:
321  T* MADNESS_RESTRICT _p;
322  TENSOR_SHARED_PTR <T> _shptr;
323 
324  void allocate(long nd, const long d[], bool dozero) {
326  if (nd < 0) {
327  _p = 0;
328  _shptr.reset();
329  _size = 0;
330  _ndim = -1;
331  return;
332  }
333 
334  TENSOR_ASSERT(nd>0 && nd <= TENSOR_MAXDIM,"invalid ndim in new tensor", nd, 0);
335  // sanity check ... 2GB in doubles
336  for (int i=0; i<nd; ++i) {
337  TENSOR_ASSERT(d[i]>=0 && d[i]<268435456, "invalid dimension size in new tensor",d[i],0);
338  }
339  set_dims_and_size(nd, d);
340  if (_size) {
341  TENSOR_ASSERT(_size>=0 && _size<268435456, "invalid size in new tensor",_size,0);
342  try {
343 #if HAVE_IBMBGP
344 #define TENSOR_ALIGNMENT 16
345 #elif HAVE_IBMBGQ
346 #define TENSOR_ALIGNMENT 32
347 #elif MADNESS_HAVE_AVX2
348 /* 32B alignment is best for performance according to
349  * http://www.nas.nasa.gov/hecc/support/kb/haswell-processors_492.html */
350 #define TENSOR_ALIGNMENT 32
351 #elif MADNESS_HAVE_AVX512
352 /* One can infer from the AVX2 case that 64B alignment helps with 512b SIMD. */
353 #define TENSOR_ALIGNMENT 64
354 #else
355 // Typical cache line size
356 #define TENSOR_ALIGNMENT 64
357 #endif
358 
359 #ifdef TENSOR_USE_SHARED_ALIGNED_ARRAY
360  _p = _shptr.allocate(_size, TENSOR_ALIGNMENT);
361 #elif defined WORLD_GATHER_MEM_STATS
362  _p = new T[_size];
363  _shptr = std::shared_ptr<T>(_p);
364 #else
365  if (posix_memalign((void **) &_p, TENSOR_ALIGNMENT, sizeof(T)*_size)) throw 1;
366  _shptr.reset(_p, &free);
367 #endif
368  }
369  catch (...) {
370  // Ideally use if constexpr here but want headers C++14 for cuda compatibility
371 MADNESS_PRAGMA_GCC(diagnostic push)
372 MADNESS_PRAGMA_GCC(diagnostic ignored "-Warray-bounds")
373  std::printf("new failed nd=%ld type=%ld size=%ld\n", nd, id(), _size);
374  std::printf(" %ld %ld %ld %ld %ld %ld\n",
375  d[0], d[1], d[2], d[3], d[4], d[5]);
376 MADNESS_PRAGMA_GCC(diagnostic pop)
377  TENSOR_EXCEPTION("new failed",_size,this);
378  }
379  //std::printf("allocated %p [%ld] %ld\n", _p, size, p.use_count());
380  if (dozero) {
381  //T zero = 0; for (long i=0; i<_size; ++i) _p[i] = zero;
382  // or
383 #ifdef HAVE_MEMSET
384  memset((void *) _p, 0, _size*sizeof(T));
385 #else
387 #endif
388  }
389  }
390  else {
391  _p = 0;
392  _shptr.reset();
393  }
394  }
395 
396  // Free memory and restore default constructor state
397  void deallocate() {
398  _p = 0;
399  _shptr.reset();
400  _size = 0;
401  _ndim = -1;
402  }
403 
404  public:
405  /// C++ typename of this tensor.
406  typedef T type;
407 
408  /// C++ typename of the real type associated with a complex type.
410 
411  /// C++ typename of the floating point type associated with scalar real type
413 
414  /// Default constructor does not allocate any data and sets ndim=-1, size=0, _p=0, and id.
415  Tensor() : _p(0) {
417  }
418 
419  /// Copy constructor is shallow (same as assignment)
420 
421  /// \em Caveat \em emptor: The shallow copy constructor has many virtues but
422  /// enables you to violate constness with simple code such as
423  /// \code
424  /// const Tensor<double> a(5);
425  /// Tensor<double> b(a);
426  /// b[1] = 3; // a[1] is now also 3
427  /// \endcode
428  Tensor(const Tensor<T>& t) {
430  *this = t;
431  }
432 
433  /// Assignment is shallow (same as copy constructor)
434 
435  /// \em Caveat \em emptor: The shallow assignment has many virtues but
436  /// enables you to violate constness with simple code such as
437  /// \code
438  /// const Tensor<double> a(5);
439  /// Tensor<double> b;
440  /// b = a;
441  /// b[1] = 3; // a[1] is now also 3
442  /// \endcode
444  if (this != &t) {
445  _p = t._p;
446  _shptr = t._shptr;
447  _size = t._size;
448  _ndim = t._ndim;
449 MADNESS_PRAGMA_GCC(diagnostic push)
450 MADNESS_PRAGMA_GCC(diagnostic ignored "-Wmaybe-uninitialized")
451  for (int i=0; i<TENSOR_MAXDIM; ++i) {
452  _dim[i] = t._dim[i];
453  _stride[i] = t._stride[i];
454  }
455 MADNESS_PRAGMA_GCC(diagnostic pop)
456  }
457  return *this;
458  }
459 
460 
461  /// Type conversion makes a deep copy
462  template <class Q> operator Tensor<Q>() const { // type conv => deep copy
463  Tensor<Q> result = Tensor<Q>(this->_ndim,this->_dim,false);
464  BINARY_OPTIMIZED_ITERATOR(Q, result, const T, (*this), *_p0 = (Q)(*_p1));
465  return result;
466  }
467 
468 
469  /// Create and zero new 1-d tensor
470 
471  /// @param[in] d0 Size of dimension 0
472  explicit Tensor(long d0) : _p(0) {
473  _dim[0] = d0;
474  allocate(1, _dim, true);
475  }
476 
477  /// Create and zero new 2-d tensor
478 
479  /// @param[in] d0 Size of dimension 0
480  /// @param[in] d1 Size of dimension 1
481  explicit Tensor(long d0, long d1) : _p(0) {
482  _dim[0] = d0; _dim[1] = d1;
483  allocate(2, _dim, true);
484  }
485 
486  /// Create and zero new 3-d tensor
487 
488  /// @param[in] d0 Size of dimension 0
489  /// @param[in] d1 Size of dimension 1
490  /// @param[in] d2 Size of dimension 2
491  explicit Tensor(long d0, long d1, long d2) : _p(0) {
492  _dim[0] = d0; _dim[1] = d1; _dim[2] = d2;
493  allocate(3, _dim, true);
494  }
495 
496  /// Create and zero new 4-d tensor
497 
498  /// @param[in] d0 Size of dimension 0
499  /// @param[in] d1 Size of dimension 1
500  /// @param[in] d2 Size of dimension 2
501  /// @param[in] d3 Size of dimension 3
502  explicit Tensor(long d0, long d1, long d2, long d3) : _p(0) {
503  _dim[0] = d0; _dim[1] = d1; _dim[2] = d2; _dim[3] = d3;
504  allocate(4, _dim, true);
505  }
506 
507  /// Create and zero new 5-d tensor
508 
509  /// @param[in] d0 Size of dimension 0
510  /// @param[in] d1 Size of dimension 1
511  /// @param[in] d2 Size of dimension 2
512  /// @param[in] d3 Size of dimension 3
513  /// @param[in] d4 Size of dimension 4
514  explicit Tensor(long d0, long d1, long d2, long d3, long d4) : _p(0) {
515  _dim[0] = d0; _dim[1] = d1; _dim[2] = d2; _dim[3] = d3; _dim[4] = d4;
516  allocate(5, _dim, true);
517  }
518 
519  /// Create and zero new 6-d tensor
520 
521  /// @param[in] d0 Size of dimension 0
522  /// @param[in] d1 Size of dimension 1
523  /// @param[in] d2 Size of dimension 2
524  /// @param[in] d3 Size of dimension 3
525  /// @param[in] d4 Size of dimension 4
526  /// @param[in] d5 Size of dimension 5
527  explicit Tensor(long d0, long d1, long d2, long d3, long d4, long d5) {
528  _dim[0] = d0; _dim[1] = d1; _dim[2] = d2; _dim[3] = d3; _dim[4] = d4; _dim[5] = d5;
529  allocate(6, _dim, true);
530  }
531 
532  /// Create and optionally zero new n-d tensor. This is the most general constructor.
533 
534  /// @param[in] d Vector containing size of each dimension, number of dimensions inferred from vector size.
535  /// @param[in] dozero If true (default) the tensor is initialized to zero
536  explicit Tensor(const std::vector<long>& d, bool dozero=true) : _p(0) {
537  allocate(d.size(), d.size() ? &(d[0]) : 0, dozero);
538  }
539 
540  /// Politically incorrect general constructor.
541 
542  /// @param[in] nd Number of dimensions
543  /// @param[in] d Size of each dimension
544  /// @param[in] dozero If true (default) the tensor is initialized to zero
545  explicit Tensor(long nd, const long d[], bool dozero=true) : _p(0) {
546  allocate(nd,d,dozero);
547  }
548 
549  /// Inplace fill tensor with scalar
550 
551  /// @param[in] x Value used to fill tensor via assigment
552  /// @return %Reference to this tensor
554  UNARY_OPTIMIZED_ITERATOR(T,(*this),*_p0 = x);
555  return *this;
556  }
557 
558  /// Inplace fill with a scalar (legacy name)
559 
560  /// @param[in] x Value used to fill tensor via assigment
561  /// @return %Reference to this tensor
563  *this = x;
564  return *this;
565  }
566 
567  /// Inplace addition of two tensors
568 
569  /// @param[in] t Conforming tensor to be added in-place to this tensor
570  /// @return %Reference to this tensor
571  template <typename Q>
573  BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 += *_p1);
574  return *this;
575  }
576 
577  /// Inplace subtraction of two tensors
578 
579  /// @param[in] t Conforming tensor to be subtracted in-place from this tensor
580  /// @return %Reference to this tensor
581  template <typename Q>
583  BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 -= *_p1);
584  return *this;
585  }
586 
587  /// Addition of two tensors to produce a new tensor
588 
589  /// @param[in] t Conforming tensor to be added out-of-place to this tensor
590  /// @return New tensor
591  template <typename Q>
593  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
594  Tensor<resultT> result(_ndim,_dim,false);
595  TERNARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), const Q, t, *_p0 = *_p1 + *_p2);
596  return result;
597  }
598 
599  /// Subtraction of two tensors to produce a new tensor
600 
601  /// @param[in] t Conforming tensor to be subtracted out-of-place from this tensor
602  /// @return New tensor
603  template <typename Q>
605  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
606  Tensor<resultT> result(_ndim,_dim,false);
607  TERNARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), const Q, t, *_p0 = *_p1 - *_p2);
608  return result;
609  }
610 
611  /// Multiplication of tensor by a scalar of a supported type to produce a new tensor
612 
613  /// @param[in] x Scalar value
614  /// @return New tensor
615  template <typename Q>
617  operator*(const Q& x) const {
618  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
619  Tensor<resultT> result(_ndim,_dim,false);
620  BINARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), *_p0 = *_p1 * x);
621  return result;
622  }
623 
624  /// Divide tensor by a scalar of a supported type to produce a new tensor
625 
626  /// @param[in] x Scalar value
627  /// @return New tensor
628  template <typename Q>
630  operator/(const Q& x) const {
631  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
632  Tensor<resultT> result(_ndim,_dim);
633  BINARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), *_p0 = *_p1 / x);
634  return result;
635  }
636 
637  /// Add a scalar of the same type to all elements of a tensor producing a new tensor
638 
639  /// @param[in] x Scalar value
640  /// @return New tensor
641  template <typename Q>
643  operator+(const Q& x) const {
644  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
645  Tensor<resultT> result(_ndim,_dim);
646  BINARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), *_p0 = *_p1 + x);
647  return result;
648  }
649 
650  /// Subtract a scalar of the same type from all elements producing a new tensor
651 
652  /// @param[in] x Scalar value
653  /// @return New tensor
654  template <typename Q>
656  operator-(const Q& x) const {
657  return (*this) + (-x);
658  }
659 
660  /// Unary negation producing a new tensor
661 
662  /// @return New tensor
664  Tensor<T> result = Tensor<T>(_ndim,_dim,false);
665  BINARY_OPTIMIZED_ITERATOR(T, result, const T, (*this), *(_p0) = - (*_p1));
666  return result;
667  }
668 
669  /// Inplace multiplication by scalar of supported type
670 
671  /// @param[in] x Scalar value
672  /// @return %Reference to this tensor
673  template <typename Q>
675  operator*=(const Q& x) {
676  UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 *= x);
677  return *this;
678  }
679 
680  /// Inplace multiplication by scalar of supported type (legacy name)
681 
682  /// @param[in] x Scalar value
683  /// @return %Reference to this tensor
684  template <typename Q>
686  scale(Q x) {
687  return (*this)*=x;
688  }
689 
690  /// Inplace increment by scalar of supported type
691 
692  /// @param[in] x Scalar value
693  /// @return %Reference to this tensor
694  template <typename Q>
696  operator+=(const Q& x) {
697  UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 += x);
698  return *this;
699  }
700 
701  /// Inplace decrement by scalar of supported type
702 
703  /// @param[in] x Scalar value
704  /// @return %Reference to this tensor
705  template <typename Q>
707  operator-=(const Q& x) {
708  UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 -= x);
709  return *this;
710  }
711 
712 
713  /// Inplace complex conjugate
714 
715  /// @return %Reference to this tensor
717  UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 = conditional_conj(*_p0));
718  return *this;
719  }
720 
721  /// Inplace fill with random values ( \c [0,1] for floats, \c [0,MAXSIZE] for integers)
722 
723  /// @return %Reference to this tensor
725  if (iscontiguous()) {
726  madness::RandomVector<T>(size(), ptr());
727  }
728  else {
729  UNARY_OPTIMIZED_ITERATOR(T,(*this), *_p0 = madness::RandomValue<T>());
730  }
731  return *this;
732  }
733 
734  /// Inplace fill with the index of each element
735 
736  /// Each element is assigned it's logical index according to this loop structure
737  /// \code
738  /// Tensor<float> t(5,6,7,...)
739  /// long index=0;
740  /// for (long i=0; i<_dim[0]; ++i)
741  /// for (long j=0; j<_dim[1]; ++j)
742  /// for (long k=0; k<_dim[2]; ++k)
743  /// ...
744  /// tensor(i,j,k,...) = index++
745  /// \endcode
746  ///
747  /// @return %Reference to this tensor
749  long count = 0;
750  UNARY_UNOPTIMIZED_ITERATOR(T,(*this), *_p0 = count++); // Fusedim would be OK
751  return *this;
752  }
753 
754  /// Inplace set elements of \c *this less than \c x in absolute magnitude to zero.
755 
756  /// @param[in] x Scalar value
757  /// @return %Reference to this tensor
758  Tensor<T>& screen(double x) {
759  T zero = 0;
760  UNARY_OPTIMIZED_ITERATOR(T,(*this), if (std::abs(*_p0)<x) *_p0=zero);
761  return *this;
762  }
763 
764 
765  /// Return true if bounds checking was enabled at compile time
766 
767  /// @return True if bounds checking was enabled at compile time
768  static bool bounds_checking() {
769 #ifdef TENSOR_BOUNDS_CHECKING
770  return true;
771 #else
772  return false;
773 #endif
774  }
775 
776  /// 1-d indexing operation using \c [] \em without bounds checking.
777 
778  /// @param[in] i index for dimension 0
779  /// @return %Reference to element
780  T& operator[](long i) {
781 #ifdef TENSOR_BOUNDS_CHECKING
782  TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
783 #endif
784  return _p[i*_stride[0]];
785  }
786 
787  /// 1-d indexing operation using \c [] \em without bounds checking.
788 
789  /// @param[in] i index for dimension 0
790  /// @return %Reference to element
791  const T& operator[](long i) const {
792 #ifdef TENSOR_BOUNDS_CHECKING
793  TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
794 #endif
795  return _p[i*_stride[0]];
796  }
797 
798  /// 1-d indexing operation \em without bounds checking.
799 
800  /// @param[in] i index for dimension 0
801  /// @return %Reference to element
802  T& operator()(long i) {
803 #ifdef TENSOR_BOUNDS_CHECKING
804  TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
805 #endif
806  return _p[i*_stride[0]];
807  }
808 
809  /// 1-d indexing operation \em without bounds checking.
810 
811  /// @param[in] i index for dimension 0
812  /// @return %Reference to element
813  const T& operator()(long i) const {
814 #ifdef TENSOR_BOUNDS_CHECKING
815  TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
816 #endif
817  return _p[i*_stride[0]];
818  }
819 
820  /// 2-d indexing operation \em without bounds checking.
821 
822  /// @param[in] i index for dimension 0
823  /// @param[in] j index for dimension 1
824  /// @return %Reference to element
825  T& operator()(long i, long j) {
826 #ifdef TENSOR_BOUNDS_CHECKING
827  TENSOR_ASSERT(i>=0 && i<_dim[0],"2d bounds check failed dim=0",i,this);
828  TENSOR_ASSERT(j>=0 && j<_dim[1],"2d bounds check failed dim=1",j,this);
829 #endif
830  return _p[i*_stride[0]+j*_stride[1]];
831  }
832 
833  /// 2-d indexing operation \em without bounds checking.
834 
835  /// @param[in] i index for dimension 0
836  /// @param[in] j index for dimension 1
837  /// @return %Reference to element
838  const T& operator()(long i, long j) const {
839 #ifdef TENSOR_BOUNDS_CHECKING
840  TENSOR_ASSERT(i>=0 && i<_dim[0],"2d bounds check failed dim=0",i,this);
841  TENSOR_ASSERT(j>=0 && j<_dim[1],"2d bounds check failed dim=1",j,this);
842 #endif
843  return _p[i*_stride[0]+j*_stride[1]];
844  }
845 
846  /// 3-d indexing operation \em without bounds checking.
847 
848  /// @param[in] i index for dimension 0
849  /// @param[in] j index for dimension 1
850  /// @param[in] k index for dimension 2
851  /// @return %Reference to element
852  T& operator()(long i, long j, long k) {
853 #ifdef TENSOR_BOUNDS_CHECKING
854  TENSOR_ASSERT(i>=0 && i<_dim[0],"3d bounds check failed dim=0",i,this);
855  TENSOR_ASSERT(j>=0 && j<_dim[1],"3d bounds check failed dim=1",j,this);
856  TENSOR_ASSERT(k>=0 && k<_dim[2],"3d bounds check failed dim=2",k,this);
857 #endif
858  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]];
859  }
860 
861  /// 3-d indexing operation \em without bounds checking.
862 
863  /// @param[in] i index for dimension 0
864  /// @param[in] j index for dimension 1
865  /// @param[in] k index for dimension 2
866  /// @return %Reference to element
867  const T& operator()(long i, long j, long k) const {
868 #ifdef TENSOR_BOUNDS_CHECKING
869  TENSOR_ASSERT(i>=0 && i<_dim[0],"3d bounds check failed dim=0",i,this);
870  TENSOR_ASSERT(j>=0 && j<_dim[1],"3d bounds check failed dim=1",j,this);
871  TENSOR_ASSERT(k>=0 && k<_dim[2],"3d bounds check failed dim=2",k,this);
872 #endif
873  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]];
874  }
875 
876  /// 4-d indexing operation \em without bounds checking.
877 
878  /// @param[in] i index for dimension 0
879  /// @param[in] j index for dimension 1
880  /// @param[in] k index for dimension 2
881  /// @param[in] l index for dimension 3
882  /// @return %Reference to element
883  T& operator()(long i, long j, long k, long l) {
884 #ifdef TENSOR_BOUNDS_CHECKING
885  TENSOR_ASSERT(i>=0 && i<_dim[0],"4d bounds check failed dim=0",i,this);
886  TENSOR_ASSERT(j>=0 && j<_dim[1],"4d bounds check failed dim=1",j,this);
887  TENSOR_ASSERT(k>=0 && k<_dim[2],"4d bounds check failed dim=2",k,this);
888  TENSOR_ASSERT(l>=0 && l<_dim[3],"4d bounds check failed dim=3",l,this);
889 #endif
890  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
891  l*_stride[3]];
892  }
893 
894  /// 4-d indexing operation \em without bounds checking.
895 
896  /// @param[in] i index for dimension 0
897  /// @param[in] j index for dimension 1
898  /// @param[in] k index for dimension 2
899  /// @param[in] l index for dimension 3
900  /// @return %Reference to element
901  const T& operator()(long i, long j, long k, long l) const {
902 #ifdef TENSOR_BOUNDS_CHECKING
903  TENSOR_ASSERT(i>=0 && i<_dim[0],"4d bounds check failed dim=0",i,this);
904  TENSOR_ASSERT(j>=0 && j<_dim[1],"4d bounds check failed dim=1",j,this);
905  TENSOR_ASSERT(k>=0 && k<_dim[2],"4d bounds check failed dim=2",k,this);
906  TENSOR_ASSERT(l>=0 && l<_dim[3],"4d bounds check failed dim=3",l,this);
907 #endif
908  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
909  l*_stride[3]];
910  }
911 
912  /// 5-d indexing operation \em without bounds checking.
913 
914  /// @param[in] i index for dimension 0
915  /// @param[in] j index for dimension 1
916  /// @param[in] k index for dimension 2
917  /// @param[in] l index for dimension 3
918  /// @param[in] m index for dimension 4
919  /// @return %Reference to element
920  T& operator()(long i, long j, long k, long l, long m) {
921 #ifdef TENSOR_BOUNDS_CHECKING
922  TENSOR_ASSERT(i>=0 && i<_dim[0],"5d bounds check failed dim=0",i,this);
923  TENSOR_ASSERT(j>=0 && j<_dim[1],"5d bounds check failed dim=1",j,this);
924  TENSOR_ASSERT(k>=0 && k<_dim[2],"5d bounds check failed dim=2",k,this);
925  TENSOR_ASSERT(l>=0 && l<_dim[3],"5d bounds check failed dim=3",l,this);
926  TENSOR_ASSERT(m>=0 && m<_dim[4],"5d bounds check failed dim=4",m,this);
927 #endif
928  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
929  l*_stride[3]+m*_stride[4]];
930  }
931 
932  /// 5-d indexing operation \em without bounds checking.
933 
934  /// @param[in] i index for dimension 0
935  /// @param[in] j index for dimension 1
936  /// @param[in] k index for dimension 2
937  /// @param[in] l index for dimension 3
938  /// @param[in] m index for dimension 4
939  /// @return %Reference to element
940  const T& operator()(long i, long j, long k, long l, long m) const {
941 #ifdef TENSOR_BOUNDS_CHECKING
942  TENSOR_ASSERT(i>=0 && i<_dim[0],"5d bounds check failed dim=0",i,this);
943  TENSOR_ASSERT(j>=0 && j<_dim[1],"5d bounds check failed dim=1",j,this);
944  TENSOR_ASSERT(k>=0 && k<_dim[2],"5d bounds check failed dim=2",k,this);
945  TENSOR_ASSERT(l>=0 && l<_dim[3],"5d bounds check failed dim=3",l,this);
946  TENSOR_ASSERT(m>=0 && m<_dim[4],"5d bounds check failed dim=4",m,this);
947 #endif
948  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
949  l*_stride[3]+m*_stride[4]];
950  }
951 
952  /// 6-d indexing operation \em without bounds checking.
953 
954  /// @param[in] i index for dimension 0
955  /// @param[in] j index for dimension 1
956  /// @param[in] k index for dimension 2
957  /// @param[in] l index for dimension 3
958  /// @param[in] m index for dimension 4
959  /// @param[in] n index for dimension 5
960  /// @return %Reference to element
961  T& operator()(long i, long j, long k, long l, long m, long n) {
962 #ifdef TENSOR_BOUNDS_CHECKING
963  TENSOR_ASSERT(i>=0 && i<_dim[0],"6d bounds check failed dim=0",i,this);
964  TENSOR_ASSERT(j>=0 && j<_dim[1],"6d bounds check failed dim=1",j,this);
965  TENSOR_ASSERT(k>=0 && k<_dim[2],"6d bounds check failed dim=2",k,this);
966  TENSOR_ASSERT(l>=0 && l<_dim[3],"6d bounds check failed dim=3",l,this);
967  TENSOR_ASSERT(m>=0 && m<_dim[4],"6d bounds check failed dim=4",m,this);
968  TENSOR_ASSERT(n>=0 && n<_dim[5],"6d bounds check failed dim=5",n,this);
969 #endif
970  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
971  l*_stride[3]+m*_stride[4]+n*_stride[5]];
972  }
973 
974  /// 6-d indexing operation \em without bounds checking.
975 
976  /// @param[in] i index for dimension 0
977  /// @param[in] j index for dimension 1
978  /// @param[in] k index for dimension 2
979  /// @param[in] l index for dimension 3
980  /// @param[in] m index for dimension 4
981  /// @param[in] n index for dimension 5
982  /// @return %Reference to element
983  const T& operator()(long i, long j, long k, long l, long m, long n) const {
984 #ifdef TENSOR_BOUNDS_CHECKING
985  TENSOR_ASSERT(i>=0 && i<_dim[0],"6d bounds check failed dim=0",i,this);
986  TENSOR_ASSERT(j>=0 && j<_dim[1],"6d bounds check failed dim=1",j,this);
987  TENSOR_ASSERT(k>=0 && k<_dim[2],"6d bounds check failed dim=2",k,this);
988  TENSOR_ASSERT(l>=0 && l<_dim[3],"6d bounds check failed dim=3",l,this);
989  TENSOR_ASSERT(m>=0 && m<_dim[4],"6d bounds check failed dim=4",m,this);
990  TENSOR_ASSERT(n>=0 && n<_dim[5],"6d bounds check failed dim=5",n,this);
991 #endif
992  return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
993  l*_stride[3]+m*_stride[4]+n*_stride[5]];
994  }
995 
996  /// Politically incorrect general indexing operation \em without bounds checking.
997 
998  /// @param[in] ind Array containing index for each dimension
999  /// @return %Reference to element
1000  T& operator()(const long ind[]) {
1001  long offset = 0;
1002  for (int d=0; d<_ndim; ++d) {
1003  long i = ind[d];
1004 #ifdef TENSOR_BOUNDS_CHECKING
1005  TENSOR_ASSERT(i>=0 && i<_dim[0],"non-PC general indexing bounds check failed dim=",d,this);
1006 #endif
1007  offset += i*_stride[d];
1008  }
1009  return _p[offset];
1010  }
1011 
1012  /// Politically incorrect general indexing operation \em without bounds checking.
1013 
1014  /// @param[in] ind Array containing index for each dimension
1015  /// @return %Reference to element
1016  const T& operator()(const long ind[]) const {
1017  long offset = 0;
1018  for (int d=0; d<_ndim; ++d) {
1019  long i = ind[d];
1020 #ifdef TENSOR_BOUNDS_CHECKING
1021  TENSOR_ASSERT(i>=0 && i<_dim[0],"non-PC general indexing bounds check failed dim=",d,this);
1022 #endif
1023  offset += i*_stride[d];
1024  }
1025  return _p[offset];
1026  }
1027 
1028  /// General indexing operation \em with bounds checking.
1029 
1030  /// @param[in] ind Vector containing index for each dimension
1031  /// @return %Reference to element
1032  T& operator()(const std::vector<long> ind) {
1033  TENSOR_ASSERT(ind.size()>=(unsigned int) _ndim,"invalid number of dimensions",ind.size(),this);
1034  long index=0;
1035  for (long d=0; d<_ndim; ++d) {
1036  TENSOR_ASSERT(ind[d]>=0 && ind[d]<_dim[d],"out-of-bounds access",ind[d],this);
1037  index += ind[d]*_stride[d];
1038  }
1039  return _p[index];
1040  }
1041 
1042  /// General indexing operation \em with bounds checking.
1043 
1044  /// @param[in] ind Vector containing index for each dimension
1045  /// @return %Reference to element
1046  const T& operator()(const std::vector<long> ind) const {
1047  TENSOR_ASSERT(ind.size()>=(unsigned int) _ndim,"invalid number of dimensions",ind.size(),this);
1048  long index=0;
1049  for (long d=0; d<_ndim; ++d) {
1050  TENSOR_ASSERT(ind[d]>=0 && ind[d]<_dim[d],"out-of-bounds access",ind[d],this);
1051  index += ind[d]*_stride[d];
1052  }
1053  return _p[index];
1054  }
1055 
1056  /// General slicing operation
1057 
1058  /// @param[in] s Vector containing slice for each dimension
1059  /// @return SliceTensor viewing patch of original tensor
1060  SliceTensor<T> operator()(const std::vector<Slice>& s) {
1061  TENSOR_ASSERT(s.size()>=(unsigned)(this->ndim()), "invalid number of dimensions",
1062  this->ndim(),this);
1063  return SliceTensor<T>(*this,&(s[0]));
1064  }
1065 
1066  /// General slicing operation (const)
1067 
1068  /// @param[in] s Vector containing slice for each dimension
1069  /// @return Constant Tensor viewing patch of original tensor
1070  const Tensor<T> operator()(const std::vector<Slice>& s) const {
1071  TENSOR_ASSERT(s.size()>=(unsigned)(this->ndim()), "invalid number of dimensions",
1072  this->ndim(),this);
1073  return SliceTensor<T>(*this,&(s[0]));
1074  }
1075 
1076  /// General slicing operation
1077 
1078  /// @param[in] s array containing slice for each dimension
1079  /// @return SliceTensor viewing patch of original tensor
1080  SliceTensor<T> operator()(const std::array<Slice,TENSOR_MAXDIM>& s) {
1081  return SliceTensor<T>(*this,s);
1082  }
1083 
1084  /// General slicing operation (const)
1085 
1086  /// @param[in] s array containing slice for each dimension
1087  /// @return Constant Tensor viewing patch of original tensor
1088  const Tensor<T> operator()(const std::array<Slice,TENSOR_MAXDIM>& s) const {
1089  return SliceTensor<T>(*this,s);
1090  }
1091 
1092  /// Return a 1d SliceTensor that views the specified range of the 1d Tensor
1093 
1094  /// @return SliceTensor viewing patch of original tensor
1096  TENSOR_ASSERT(this->ndim()==1,"invalid number of dimensions",
1097  this->ndim(),this);
1098  Slice s[1] = {s0};
1099  return SliceTensor<T>(*this,s);
1100  }
1101 
1102  /// Return a 1d SliceTensor that views the specified range of the 1d Tensor
1103 
1104  /// @return Constant Tensor viewing patch of original tensor: \f$ R(*,*,\ldots) \rightarrow I(*,*,\ldots) \f$
1105  const Tensor<T> operator()(const Slice& s0) const {
1106  TENSOR_ASSERT(this->ndim()==1,"invalid number of dimensions",
1107  this->ndim(),this);
1108  Slice s[1] = {s0};
1109  return SliceTensor<T>(*this,s);
1110  }
1111 
1112  /// Return a 1d SliceTensor that views the specified range of the 2d Tensor
1113 
1114  /// @return SliceTensor viewing patch of original tensor: \f$ R(*) \rightarrow I(i,*) \f$
1115  SliceTensor<T> operator()(long i, const Slice& s1) {
1116  TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
1117  this->ndim(),this);
1118  Slice s[2] = {Slice(i,i,0),s1};
1119  return SliceTensor<T>(*this,s);
1120  }
1121 
1122  /// Return a 1d SliceTensor that views the specified range of the 2d Tensor
1123 
1124  /// @return Constant Tensor viewing patch of original tensor
1125  const Tensor<T> operator()(long i, const Slice& s1) const {
1126  TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
1127  this->ndim(),this);
1128  Slice s[2] = {Slice(i,i,0),s1};
1129  return SliceTensor<T>(*this,s);
1130  }
1131 
1132  /// Return a 1d SliceTensor that views the specified range of the 2d Tensor
1133 
1134  /// @return SliceTensor viewing patch of original tensor
1135  SliceTensor<T> operator()(const Slice& s0, long j) {
1136  TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
1137  this->ndim(),this);
1138  Slice s[2] = {s0,Slice(j,j,0)};
1139  return SliceTensor<T>(*this,s);
1140  }
1141 
1142  /// Return a 1d constant Tensor that views the specified range of the 2d Tensor
1143 
1144  /// @return Constant Tensor viewing patch of original tensor
1145  const Tensor<T> operator()(const Slice& s0, long j) const {
1146  TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
1147  this->ndim(),this);
1148  Slice s[2] = {s0,Slice(j,j,0)};
1149  return SliceTensor<T>(*this,s);
1150  }
1151 
1152  /// Return a 2d SliceTensor that views the specified range of the 2d Tensor
1153 
1154  /// @return SliceTensor viewing patch of original tensor
1155  SliceTensor<T> operator()(const Slice& s0, const Slice& s1) {
1156  TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
1157  this->ndim(),this);
1158  Slice s[2] = {s0,s1};
1159  return SliceTensor<T>(*this,s);
1160  }
1161 
1162  /// Return a 2d constant Tensor that views the specified range of the 2d Tensor
1163 
1164  /// @return Constant Tensor viewing patch of original tensor
1165  const Tensor<T> operator()(const Slice& s0, const Slice& s1) const {
1166  TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
1167  this->ndim(),this);
1168  Slice s[2] = {s0,s1};
1169  return SliceTensor<T>(*this,s);
1170  }
1171 
1172  /// Return a 3d SliceTensor that views the specified range of the 3d Tensor
1173 
1174  /// @return SliceTensor viewing patch of original tensor
1175  SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2) {
1176  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1177  this->ndim(),this);
1178  Slice s[3] = {s0,s1,s2};
1179  return SliceTensor<T>(*this,s);
1180  }
1181 
1182  /// Return a 3d constant Tensor that views the specified range of the 3d Tensor
1183 
1184  /// @return Constant Tensor viewing patch of original tensor
1185  const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2) const {
1186  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1187  this->ndim(),this);
1188  Slice s[3] = {s0,s1,s2};
1189  return SliceTensor<T>(*this,s);
1190  }
1191 
1192  /// Return a 2d SliceTensor that views the specified range of the 3d Tensor
1193 
1194  /// @return SliceTensor viewing patch of original tensor
1195  SliceTensor<T> operator()(long i, const Slice& s1, const Slice& s2) {
1196  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1197  this->ndim(),this);
1198  Slice s[3] = {Slice(i,i,0),s1,s2};
1199  return SliceTensor<T>(*this,s);
1200  }
1201 
1202  /// Return a 2d constant Tensor that views the specified range of the 3d Tensor
1203 
1204  /// @return Constant Tensor viewing patch of original tensor
1205  const Tensor<T> operator()(long i, const Slice& s1, const Slice& s2) const {
1206  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1207  this->ndim(),this);
1208  Slice s[3] = {Slice(i,i,0),s1,s2};
1209  return SliceTensor<T>(*this,s);
1210  }
1211 
1212  /// Return a 2d SliceTensor that views the specified range of the 3d Tensor
1213 
1214  /// @return SliceTensor viewing patch of original tensor
1215  SliceTensor<T> operator()(const Slice& s0, long j, const Slice& s2) {
1216  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1217  this->ndim(),this);
1218  Slice s[3] = {s0,Slice(j,j,0),s2};
1219  return SliceTensor<T>(*this,s);
1220  }
1221 
1222  /// Return a 2d constant Tensor that views the specified range of the 3d Tensor
1223 
1224  /// @return Constant Tensor viewing patch of original tensor
1225  const Tensor<T> operator()(const Slice& s0, long j, const Slice& s2) const {
1226  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1227  this->ndim(),this);
1228  Slice s[3] = {s0,Slice(j,j,0),s2};
1229  return SliceTensor<T>(*this,s);
1230  }
1231 
1232  /// Return a 2d SliceTensor that views the specified range of the 3d Tensor
1233 
1234  /// @return SliceTensor viewing patch of original tensor
1235  SliceTensor<T> operator()(const Slice& s0, const Slice& s1, long k) {
1236  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1237  this->ndim(),this);
1238  Slice s[3] = {s0,s1,Slice(k,k,0)};
1239  return SliceTensor<T>(*this,s);
1240  }
1241 
1242  /// Return a 2d constant Tensor that views the specified range of the 3d Tensor
1243 
1244  /// @return Constant Tensor viewing patch of original tensor
1245  const Tensor<T> operator()(const Slice& s0, const Slice& s1, long k) const {
1246  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1247  this->ndim(),this);
1248  Slice s[3] = {s0,s1,Slice(k,k,0)};
1249  return SliceTensor<T>(*this,s);
1250  }
1251 
1252  /// Return a 1d SliceTensor that views the specified range of the 3d Tensor
1253 
1254  /// @return SliceTensor viewing patch of original tensor
1255  SliceTensor<T> operator()(long i, long j, const Slice& s2) {
1256  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1257  this->ndim(),this);
1258  Slice s[3] = {Slice(i,i,0),Slice(j,j,0),s2};
1259  return SliceTensor<T>(*this,s);
1260  }
1261 
1262  /// Return a 1d constant Tensor that views the specified range of the 3d Tensor
1263 
1264  /// @return Constant Tensor viewing patch of original tensor
1265  const Tensor<T> operator()(long i, long j, const Slice& s2) const {
1266  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1267  this->ndim(),this);
1268  Slice s[3] = {Slice(i,i,0),Slice(j,j,0),s2};
1269  return SliceTensor<T>(*this,s);
1270  }
1271 
1272  /// Return a 1d SliceTensor that views the specified range of the 3d Tensor
1273 
1274  /// @return SliceTensor viewing patch of original tensor
1275  SliceTensor<T> operator()(long i, const Slice& s1, long k) {
1276  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1277  this->ndim(),this);
1278  Slice s[3] = {Slice(i,i,0),s1,Slice(k,k,0)};
1279  return SliceTensor<T>(*this,s);
1280  }
1281 
1282  /// Return a 1d constant Tensor that views the specified range of the 3d Tensor
1283 
1284  /// @return Constant Tensor viewing patch of original tensor
1285  const Tensor<T> operator()(long i, const Slice& s1, long k) const {
1286  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1287  this->ndim(),this);
1288  Slice s[3] = {Slice(i,i,0),s1,Slice(k,k,0)};
1289  return SliceTensor<T>(*this,s);
1290  }
1291 
1292  /// Return a 1d SliceTensor that views the specified range of the 3d Tensor
1293 
1294  /// @return SliceTensor viewing patch of original tensor
1295  SliceTensor<T> operator()(const Slice& s0, long j, long k) {
1296  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1297  this->ndim(),this);
1298  Slice s[3] = {s0,Slice(j,j,0),Slice(k,k,0)};
1299  return SliceTensor<T>(*this,s);
1300  }
1301 
1302  /// Return a 1d constant Tensor that views the specified range of the 3d Tensor
1303 
1304  /// @return Constant Tensor viewing patch of original tensor
1305  const Tensor<T> operator()(const Slice& s0, long j, long k) const {
1306  TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
1307  this->ndim(),this);
1308  Slice s[3] = {s0,Slice(j,j,0),Slice(k,k,0)};
1309  return SliceTensor<T>(*this,s);
1310  }
1311 
1312  /// Return a 1-4d SliceTensor that views the specified range of the 4d Tensor
1313 
1314  /// @return SliceTensor viewing patch of original tensor
1315  SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
1316  const Slice& s3) {
1317  TENSOR_ASSERT(this->ndim()==4,"invalid number of dimensions",
1318  this->ndim(),this);
1319  Slice s[4] = {s0,s1,s2,s3};
1320  return SliceTensor<T>(*this,s);
1321  }
1322 
1323  /// Return a 1-4d constant Tensor that views the specified range of the 4d Tensor
1324 
1325  /// @return Constant Tensor viewing patch of original tensor
1326  const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
1327  const Slice& s3) const {
1328  TENSOR_ASSERT(this->ndim()==4,"invalid number of dimensions",
1329  this->ndim(),this);
1330  Slice s[4] = {s0,s1,s2,s3};
1331  return SliceTensor<T>(*this,s);
1332  }
1333 
1334  /// Return a 1-5d SliceTensor that views the specified range of the 5d Tensor
1335 
1336  /// @return SliceTensor viewing patch of original tensor
1337  SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
1338  const Slice& s3, const Slice& s4) {
1339  TENSOR_ASSERT(this->ndim()==5,"invalid number of dimensions",
1340  this->ndim(),this);
1341  Slice s[5] = {s0,s1,s2,s3,s4};
1342  return SliceTensor<T>(*this,s);
1343  }
1344 
1345  /// Return a 1-5d constant Tensor that views the specified range of the 5d Tensor
1346 
1347  /// @return Constant Tensor viewing patch of original tensor
1348  const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
1349  const Slice& s3, const Slice& s4) const {
1350  TENSOR_ASSERT(this->ndim()==5,"invalid number of dimensions",
1351  this->ndim(),this);
1352  Slice s[5] = {s0,s1,s2,s3,s4};
1353  return SliceTensor<T>(*this,s);
1354  }
1355 
1356  /// Return a 1-6d SliceTensor that views the specified range of the 6d Tensor
1357 
1358  /// @return SliceTensor viewing patch of original tensor
1359  SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
1360  const Slice& s3, const Slice& s4, const Slice& s5) {
1361  TENSOR_ASSERT(this->ndim()==6,"invalid number of dimensions",
1362  this->ndim(),this);
1363  Slice s[6] = {s0,s1,s2,s3,s4,s5};
1364  return SliceTensor<T>(*this,s);
1365  }
1366 
1367 
1368  /// Return a 1-6d constant Tensor that views the specified range of the 6d Tensor
1369 
1370  /// @return Constant Tensor viewing patch of original tensor
1371  const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
1372  const Slice& s3, const Slice& s4, const Slice& s5) const {
1373  TENSOR_ASSERT(this->ndim()==6,"invalid number of dimensions",
1374  this->ndim(),this);
1375  Slice s[6] = {s0,s1,s2,s3,s4,s5};
1376  return SliceTensor<T>(*this,s);
1377  }
1378 
1379  /// Returns new view/tensor reshaping size/number of dimensions to conforming tensor
1380 
1381  /// @param[in] ndimnew Number of dimensions in the result
1382  /// @param[in] d Array containing size of each new dimension
1383  /// @return New tensor (viewing same underlying data as the original but with different shape)
1384  Tensor<T> reshape(int ndimnew, const long* d) {
1385  Tensor<T> result(*this);
1386  result.reshape_inplace(ndimnew,d);
1387  return result;
1388  }
1389 
1390  /// Returns new view/tensor reshaping size/number of dimensions to conforming tensor
1391 
1392  /// @param[in] ndimnew Number of dimensions in the result
1393  /// @param[in] d Array containing size of each new dimension
1394  /// @return New tensor (viewing same underlying data as the original but with different shape)
1395  const Tensor<T> reshape(int ndimnew, const long* d) const {
1396  Tensor<T> result(*const_cast<Tensor<T>*>(this));
1397  result.reshape_inplace(ndimnew,d);
1398  return result;
1399  }
1400 
1401  /// Returns new view/tensor reshaping size/number of dimensions to conforming tensor
1402 
1403  /// @param[in] d Array containing size of each new dimension
1404  /// @return New tensor (viewing same underlying data as the original but with different shape)
1405  Tensor<T> reshape(const std::vector<long>& d) {
1406  return reshape(d.size(), d.size() ? &d[0] : 0);
1407  }
1408 
1409  /// Returns new view/tensor reshaping size/number of dimensions to conforming tensor
1410 
1411  /// @param[in] d Array containing size of each new dimension
1412  /// @return New tensor (viewing same underlying data as the original but with different shape)
1413  const Tensor<T> reshape(const std::vector<long>& d) const {
1414  return reshape(d.size(), d.size() ? &d[0] : 0);
1415  }
1416 
1417  /// Returns new view/tensor rehapings to conforming 1-d tensor with given dimension
1418 
1419  /// @param[in] dim0 Size of new dimension 0
1420  /// @return New tensor (viewing same underlying data as the original but with different shape)
1421  Tensor<T> reshape(long dim0) {
1422  long d[1] = {dim0};
1423  return reshape(1,d);
1424  }
1425  /// Returns new view/tensor rehapings to conforming 1-d tensor with given dimension
1426 
1427  /// @param[in] dim0 Size of new dimension 0
1428  /// @return New tensor (viewing same underlying data as the original but with different shape)
1429  const Tensor<T> reshape(long dim0) const {
1430  long d[1] = {dim0};
1431  return reshape(1,d);
1432  }
1433 
1434  /// Returns new view/tensor rehaping to conforming 2-d tensor with given dimensions
1435 
1436  /// @param[in] dim0 Size of new dimension 0
1437  /// @param[in] dim1 Size of new dimension 1
1438  /// @return New tensor (viewing same underlying data as the original but with different shape)
1439  Tensor<T> reshape(long dim0, long dim1) {
1440  long d[2] = {dim0,dim1};
1441  return reshape(2,d);
1442  }
1443 
1444  /// Returns new view/tensor rehaping to conforming 2-d tensor with given dimensions
1445 
1446  /// @param[in] dim0 Size of new dimension 0
1447  /// @param[in] dim1 Size of new dimension 1
1448  /// @return New tensor (viewing same underlying data as the original but with different shape)
1449  const Tensor<T> reshape(long dim0, long dim1) const {
1450  long d[2] = {dim0,dim1};
1451  return reshape(2,d);
1452  }
1453 
1454  /// Returns new view/tensor rehaping to conforming 3-d tensor with given dimensions
1455 
1456  /// @param[in] dim0 Size of new dimension 0
1457  /// @param[in] dim1 Size of new dimension 1
1458  /// @param[in] dim2 Size of new dimension 2
1459  /// @return New tensor (viewing same underlying data as the original but with different shape)
1460  Tensor<T> reshape(long dim0, long dim1, long dim2) {
1461  long d[3] = {dim0,dim1,dim2};
1462  return reshape(3,d);
1463  }
1464 
1465  /// Returns new view/tensor rehaping to conforming 3-d tensor with given dimensions
1466 
1467  /// @param[in] dim0 Size of new dimension 0
1468  /// @param[in] dim1 Size of new dimension 1
1469  /// @param[in] dim2 Size of new dimension 2
1470  /// @return New tensor (viewing same underlying data as the original but with different shape)
1471  const Tensor<T> reshape(long dim0, long dim1, long dim2) const {
1472  long d[3] = {dim0,dim1,dim2};
1473  return reshape(3,d);
1474  }
1475 
1476  /// Returns new view/tensor rehaping to conforming 4-d tensor with given dimensions
1477 
1478  /// @param[in] dim0 Size of new dimension 0
1479  /// @param[in] dim1 Size of new dimension 1
1480  /// @param[in] dim2 Size of new dimension 2
1481  /// @param[in] dim3 Size of new dimension 3
1482  /// @return New tensor (viewing same underlying data as the original but with different shape)
1483  Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3) {
1484  long d[4] = {dim0,dim1,dim2,dim3};
1485  return reshape(4,d);
1486  }
1487 
1488  /// Returns new view/tensor rehaping to conforming 4-d tensor with given dimensions
1489 
1490  /// @param[in] dim0 Size of new dimension 0
1491  /// @param[in] dim1 Size of new dimension 1
1492  /// @param[in] dim2 Size of new dimension 2
1493  /// @param[in] dim3 Size of new dimension 3
1494  /// @return New tensor (viewing same underlying data as the original but with different shape)
1495  const Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3) const {
1496  long d[4] = {dim0,dim1,dim2,dim3};
1497  return reshape(4,d);
1498  }
1499 
1500  /// Returns new view/tensor rehaping to conforming 5-d tensor with given dimensions
1501 
1502  /// @param[in] dim0 Size of new dimension 0
1503  /// @param[in] dim1 Size of new dimension 1
1504  /// @param[in] dim2 Size of new dimension 2
1505  /// @param[in] dim3 Size of new dimension 3
1506  /// @param[in] dim4 Size of new dimension 4
1507  /// @return New tensor (viewing same underlying data as the original but with different shape)
1508  Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4) {
1509  long d[5] = {dim0,dim1,dim2,dim3,dim4};
1510  return reshape(5,d);
1511  }
1512 
1513  /// Returns new view/tensor rehaping to conforming 5-d tensor with given dimensions
1514 
1515  /// @param[in] dim0 Size of new dimension 0
1516  /// @param[in] dim1 Size of new dimension 1
1517  /// @param[in] dim2 Size of new dimension 2
1518  /// @param[in] dim3 Size of new dimension 3
1519  /// @param[in] dim4 Size of new dimension 4
1520  /// @return New tensor (viewing same underlying data as the original but with different shape)
1521  const Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4) const {
1522  long d[5] = {dim0,dim1,dim2,dim3,dim4};
1523  return reshape(5,d);
1524  }
1525 
1526  /// Returns new view/tensor rehaping to conforming 6-d tensor with given dimensions
1527 
1528  /// @param[in] dim0 Size of new dimension 0
1529  /// @param[in] dim1 Size of new dimension 1
1530  /// @param[in] dim2 Size of new dimension 2
1531  /// @param[in] dim3 Size of new dimension 3
1532  /// @param[in] dim4 Size of new dimension 4
1533  /// @param[in] dim5 Size of new dimension 5
1534  /// @return New tensor (viewing same underlying data as the original but with different shape)
1535  Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4, long dim5) {
1536  long d[6] = {dim0,dim1,dim2,dim3,dim4,dim5};
1537  return reshape(6,d);
1538  }
1539 
1540  /// Returns new view/tensor rehaping to conforming 6-d tensor with given dimensions
1541 
1542  /// @param[in] dim0 Size of new dimension 0
1543  /// @param[in] dim1 Size of new dimension 1
1544  /// @param[in] dim2 Size of new dimension 2
1545  /// @param[in] dim3 Size of new dimension 3
1546  /// @param[in] dim4 Size of new dimension 4
1547  /// @param[in] dim5 Size of new dimension 5
1548  /// @return New tensor (viewing same underlying data as the original but with different shape)
1549  const Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4, long dim5) const {
1550  long d[6] = {dim0,dim1,dim2,dim3,dim4,dim5};
1551  return reshape(6,d);
1552  }
1553 
1554  /// Returns new view/tensor rehshaping to flat (1-d) tensor
1556  long d[1] = {_size};
1557  return reshape(1,d);
1558  }
1559 
1560  /// Returns new view/tensor rehshaping to flat (1-d) tensor
1561  const Tensor<T> flat() const {
1562  long d[1] = {_size};
1563  return reshape(1,d);
1564  }
1565 
1566  /// Returns new view/tensor splitting dimension \c i as \c dimi0*dimi1 to produce conforming d+1 dimension tensor
1567 
1568  /// @return New tensor (viewing same underlying data as the original but with additional dimensions)
1569  Tensor<T> splitdim(long i, long dimi0, long dimi1) {
1570  Tensor<T> result(*this);
1571  result.splitdim_inplace(i, dimi0, dimi1);
1572  return result;
1573  }
1574 
1575  /// Returns new view/tensor splitting dimension \c i as \c dimi0*dimi1 to produce conforming d+1 dimension tensor
1576 
1577  /// @return New tensor (viewing same underlying data as the original but with additional dimensions)
1578  const Tensor<T> splitdim(long i, long dimi0, long dimi1) const {
1579  Tensor<T> result(*const_cast<Tensor<T>*>(this));
1580  result.splitdim_inplace(i, dimi0, dimi1);
1581  return result;
1582  }
1583 
1584  /// Returns new view/tensor fusing contiguous dimensions \c i and \c i+1
1585 
1586  /// @return New tensor (viewing same underlying data as the original but with fewer dimensions)
1587  Tensor<T> fusedim(long i) {
1588  Tensor<T> result(*this);
1589  result.fusedim_inplace(i);
1590  return result;
1591  }
1592 
1593  /// Returns new view/tensor fusing contiguous dimensions \c i and \c i+1
1594 
1595  /// @return New tensor (viewing same underlying data as the original but with fewer dimensions)
1596  const Tensor<T> fusedim(long i) const {
1597  Tensor<T> result(*const_cast<Tensor<T>*>(this));
1598  result.fusedim_inplace(i);
1599  return result;
1600  }
1601 
1602  /// Returns new view/tensor swaping dimensions \c i and \c j
1603 
1604  /// @return New tensor (viewing same underlying data as the original but with reordered dimensions)
1605  Tensor<T> swapdim(long idim, long jdim) {
1606  Tensor<T> result(*this);
1607  result.swapdim_inplace(idim, jdim);
1608  return result;
1609  }
1610 
1611  /// Returns new view/tensor swaping dimensions \c i and \c j
1612 
1613  /// @return New tensor (viewing same underlying data as the original but with reordered dimensions)
1614  const Tensor<T> swapdim(long idim, long jdim) const {
1615  Tensor<T> result(*const_cast<Tensor<T>*>(this));
1616  result.swapdim_inplace(idim, jdim);
1617  return result;
1618  }
1619 
1620  /// Returns new view/tensor permuting the dimensions
1621 
1622  /// @param[in] map Old dimension i becomes new dimension \c map[i]
1623  /// @return New tensor (viewing same underlying data as the original but with reordered dimensions)
1624  Tensor<T> mapdim(const std::vector<long>& map) {
1625  Tensor<T> result(*this);
1626  result.mapdim_inplace(map);
1627  return result;
1628  }
1629 
1630  /// Returns new view/tensor permuting the dimensions
1631 
1632  /// @return New tensor (viewing same underlying data as the original but with reordered dimensions)
1633  const Tensor<T> mapdim(const std::vector<long>& map) const {
1634  Tensor<T> result(*const_cast<Tensor<T>*>(this));
1635  result.mapdim_inplace(map);
1636  return result;
1637  }
1638 
1639 
1640  /// Returns new view/tensor cycling the sub-dimensions `(start,...,end)` with `shift` steps
1641  Tensor<T> cycledim(long nshift, long start, long end) {
1642  Tensor<T> result(*this);
1643  result.cycledim_inplace(nshift, start, end);
1644  return result;
1645  }
1646 
1647 
1648  /// Returns new view/tensor cycling the sub-dimensions `(start,...,end)` with `shift` steps
1649  const Tensor<T> cycledim(long nshift, long start, long end) const {
1650  Tensor<T> result(*const_cast<Tensor<T>*>(this));
1651  result.cycledim_inplace(nshift, start, end);
1652  return result;
1653  }
1654 
1655 
1656  /// Test if \c *this and \c t conform.
1657  template <class Q> bool conforms(const Tensor<Q>& t) const {
1658  return BaseTensor::conforms(&t);
1659  }
1660 
1661  /// Returns the sum of all elements of the tensor
1662  T sum() const {
1663  T result = 0;
1664  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result += *_p0);
1665  return result;
1666  }
1667 
1668  /// Returns the sum of the squares of the elements
1669  T sumsq() const {
1670  T result = 0;
1671  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result += (*_p0) * (*_p0));
1672  return result;
1673  }
1674 
1675  /// Return the product of all elements of the tensor
1676  T product() const {
1677  T result = 1;
1678  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result *= *_p0);
1679  return result;
1680  }
1681 
1682  /// Return the minimum value (and if ind is non-null, its index) in the Tensor
1683  T min(long* ind=0) const {
1684  T result = *(this->_p);
1685  if (ind) {
1686  for (long i=0; i<_ndim; ++i) ind[i]=0;
1687  long nd = _ndim-1;
1688  UNARY_UNOPTIMIZED_ITERATOR(const T,(*this),
1689  if (result > *_p0) {
1690  result = *_p0;
1691  for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
1692  ind[nd] = _j;
1693  }
1694  );
1695  }
1696  else {
1697  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::min<T>(result,*_p0));
1698  }
1699  return result;
1700  }
1701 
1702  /// Return the maximum value (and if ind is non-null, its index) in the Tensor
1703  T max(long* ind=0) const {
1704  T result = *(this->_p);
1705  if (ind) {
1706  for (long i=0; i<_ndim; ++i) ind[i]=0;
1707  long nd = _ndim-1;
1708  UNARY_UNOPTIMIZED_ITERATOR(const T,(*this),
1709  if (result < *_p0) {
1710  result = *_p0;
1711  for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
1712  ind[nd] = _j;
1713  }
1714  );
1715  }
1716  else {
1717  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::max<T>(result,*_p0));
1718  }
1719  return result;
1720  }
1721 
1722  // For complex types, this next group returns the appropriate real type
1723  // For real types, the same type as T is returned (type_data.h)
1724 
1725  /// Returns the Frobenius norm of the tensor
1727  float_scalar_type result = 0;
1728  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result += ::madness::detail::mynorm(*_p0));
1729  return (float_scalar_type) std::sqrt(result);
1730  }
1731 
1732  /// Return the absolute minimum value (and if ind is non-null, its index) in the Tensor
1733  scalar_type absmin(long *ind = 0) const {
1734  scalar_type result = std::abs(*(this->_p));
1735  if (ind) {
1736  for (long i=0; i<_ndim; ++i) ind[i]=0;
1737  long nd = _ndim-1;
1738  UNARY_UNOPTIMIZED_ITERATOR(const T,(*this),
1739  scalar_type absval = std::abs(*_p0);
1740  if (result > absval) {
1741  result = absval;
1742  for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
1743  ind[nd] = _j;
1744  }
1745  );
1746  }
1747  else {
1748  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::min<scalar_type>(result,std::abs(*_p0)));
1749  }
1750  return result;
1751  }
1752 
1753  /// Return the absolute maximum value (and if ind is non-null, its index) in the Tensor
1754  scalar_type absmax(long *ind = 0) const {
1755  scalar_type result = std::abs(*(this->_p));
1756  if (ind) {
1757  for (long i=0; i<_ndim; ++i) ind[i]=0;
1758  long nd = _ndim-1;
1759  UNARY_UNOPTIMIZED_ITERATOR(T,(*this),
1760  scalar_type absval = std::abs(*_p0);
1761  if (result < absval) {
1762  result = absval;
1763  for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
1764  ind[nd] = _j;
1765  }
1766  );
1767  }
1768  else {
1769  UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::max<scalar_type>(result,std::abs(*_p0)));
1770  }
1771  return result;
1772  }
1773 
1774 
1775  /// Return the trace of two tensors (no complex conjugate invoked)
1776  T trace(const Tensor<T>& t) const {
1777  T result = 0;
1778  BINARY_OPTIMIZED_ITERATOR(const T,(*this),const T,t,result += (*_p0)*(*_p1));
1779  return result;
1780  }
1781 
1782  /// Return the trace of two tensors with complex conjugate of the leftmost (i.e., this)
1783  template <class Q>
1784  TENSOR_RESULT_TYPE(T,Q) trace_conj(const Tensor<Q>& t) const {
1785  TENSOR_RESULT_TYPE(T,Q) result = 0;
1786  BINARY_OPTIMIZED_ITERATOR(const T,(*this),const Q,t,result += conditional_conj(*_p0)*(*_p1));
1787  return result;
1788  }
1789 
1790  /// Inplace apply a unary function to each element of the tensor
1791  template <typename opT>
1793  UNARY_OPTIMIZED_ITERATOR(T,(*this),*_p0=op(*_p0));
1794  return *this;
1795  }
1796 
1797  /// Inplace multiply by corresponding elements of argument Tensor
1798  Tensor<T>& emul(const Tensor<T>& t) {
1799  BINARY_OPTIMIZED_ITERATOR(T,(*this),const T,t,*_p0 *= *_p1);
1800  return *this;
1801  }
1802 
1803  /// Inplace generalized saxpy ... this = this*alpha + other*beta
1805  if (iscontiguous() && t.iscontiguous()) {
1806  T* MADNESS_RESTRICT a = ptr();
1807  const T* MADNESS_RESTRICT b = t.ptr();
1808  if (alpha == T(1.0)) {
1809  for (long i=0; i<_size; ++i) a[i] += b[i]*beta;
1810  }
1811  else {
1812  for (long i=0; i<_size; ++i) a[i] = a[i]*alpha + b[i]*beta;
1813  }
1814  }
1815  else {
1816  //BINARYITERATOR(T,(*this),T,t, (*_p0) = alpha*(*_p0) + beta*(*_p1));
1817  BINARY_OPTIMIZED_ITERATOR(T,(*this),const T,t, (*_p0) = alpha*(*_p0) + beta*(*_p1));
1818  //ITERATOR((*this),(*this)(IND) = alpha*(*this)(IND) + beta*t(IND));
1819  }
1820  return *this;
1821  }
1822 
1823  /// Returns a pointer to the internal data
1824  T* ptr() {
1825  return _p;
1826  }
1827 
1828  /// Returns a pointer to the internal data
1829  const T* ptr() const {
1830  return _p;
1831  }
1832 
1833  /// Returns a pointer to the base class
1835  return static_cast<BaseTensor*>(this);
1836  }
1837 
1838  /// Returns a pointer to the base class
1839  const BaseTensor* base() const {
1840  return static_cast<BaseTensor*>(this);
1841  }
1842 
1843  /// Return iterator over single tensor
1845  bool optimize=true,
1846  bool fusedim=true,
1847  long jdim=default_jdim) const {
1848  return TensorIterator<T>(this,(const Tensor<T>*) 0, (const Tensor<T>*) 0,
1849  iterlevel, optimize, fusedim, jdim);
1850  }
1851 
1852  /// Return iterator over two tensors
1853  template <class Q>
1855  long iterlevel=0,
1856  bool optimize=true,
1857  bool fusedim=true,
1858  long jdim=default_jdim) const {
1859  return TensorIterator<T,Q>(this,&q,(const Tensor<T>*) 0,
1860  iterlevel, optimize, fusedim, jdim);
1861  }
1862 
1863  /// Return iterator over three tensors
1864  template <class Q, class R>
1866  const Tensor<R>& r,
1867  long iterlevel=0,
1868  bool optimize=true,
1869  bool fusedim=true,
1870  long jdim=default_jdim) const {
1871  return TensorIterator<T,Q,R>(this,&q,&r,
1872  iterlevel, optimize, fusedim, jdim);
1873  }
1874 
1875  /// End point for forward iteration
1876  const TensorIterator<T>& end() const {
1877  static TensorIterator<T> theend(0,0,0,0,0,0);
1878  return theend;
1879  }
1880 
1881  virtual ~Tensor() {}
1882 
1883  /// Frees all memory and resests to state of default constructor
1884  void clear() {deallocate();}
1885 
1886  bool has_data() const {return size()!=0;};
1887 
1888  };
1889 
1890  template <class T>
1891  std::ostream& operator << (std::ostream& out, const Tensor<T>& t);
1892 
1893 
1894  namespace archive {
1895  /// Serialize a tensor
1896  template <class Archive, typename T>
1897  struct ArchiveStoreImpl< Archive, Tensor<T> > {
1898  static void store(const Archive& s, const Tensor<T>& t) {
1899  if (t.iscontiguous()) {
1900  s & t.size() & t.id();
1901  if (t.size()) s & t.ndim() & wrap(t.dims(),TENSOR_MAXDIM) & wrap(t.ptr(),t.size());
1902  }
1903  else {
1904  s & copy(t);
1905  }
1906  };
1907  };
1908 
1909 
1910  /// Deserialize a tensor ... existing tensor is replaced
1911  template <class Archive, typename T>
1912  struct ArchiveLoadImpl< Archive, Tensor<T> > {
1913  static void load(const Archive& s, Tensor<T>& t) {
1914  long sz = 0l, id = 0l;
1915  s & sz & id;
1916  if (id != t.id()) throw "type mismatch deserializing a tensor";
1917  if (sz) {
1918  long _ndim = 0l, _dim[TENSOR_MAXDIM];
1919  s & _ndim & wrap(_dim,TENSOR_MAXDIM);
1920  t = Tensor<T>(_ndim, _dim, false);
1921  if (sz != t.size()) throw "size mismatch deserializing a tensor";
1922  s & wrap(t.ptr(), t.size());
1923  }
1924  else {
1925  t = Tensor<T>();
1926  }
1927  };
1928  };
1929 
1930  }
1931 
1932  /// The class defines tensor op scalar ... here define scalar op tensor.
1933 
1934  /// \ingroup tensor
1935  template <typename T, typename Q>
1936  typename IsSupported < TensorTypeData<Q>, Tensor<T> >::type
1937  operator+(Q x, const Tensor<T>& t) {
1938  return t+x;
1939  }
1940 
1941  /// The class defines tensor op scalar ... here define scalar op tensor.
1942 
1943  /// \ingroup tensor
1944  template <typename T, typename Q>
1945  typename IsSupported < TensorTypeData<Q>, Tensor<T> >::type
1946  operator*(const Q& x, const Tensor<T>& t) {
1947  return t*x;
1948  }
1949 
1950  /// The class defines tensor op scalar ... here define scalar op tensor.
1951 
1952  /// \ingroup tensor
1953  template <typename T, typename Q>
1954  typename IsSupported < TensorTypeData<Q>, Tensor<T> >::type
1955  operator-(Q x, const Tensor<T>& t) {
1956  return (-t)+=x;
1957  }
1958 
1959  /// Returns a new contiguous tensor that is a deep copy of the input
1960 
1961  /// \ingroup tensor
1962  /// @result Returns a new contiguous tensor that is a deep copy of the input
1963  template <class T> Tensor<T> copy(const Tensor<T>& t) {
1964  if (t.size()) {
1965  Tensor<T> result = Tensor<T>(t.ndim(),t.dims(),false);
1966  BINARY_OPTIMIZED_ITERATOR(T, result, const T, t, *_p0 = *_p1);
1967  return result;
1968  }
1969  else {
1970  return Tensor<T>();
1971  }
1972  }
1973 
1974  /// Returns a new contiguous tensor of type Q that is a deep copy of the input
1975 
1976  /// \ingroup tensor
1977  /// @result Returns a new contiguous tensor that is a deep copy of the input
1978  template <class Q, class T>
1980  if (t.size()) {
1981  Tensor<Q> result = Tensor<Q>(t.ndim(),t.dims(),false);
1982  BINARY_OPTIMIZED_ITERATOR(Q, result, const T, t, *_p0 = *_p1);
1983  return result;
1984  }
1985  else {
1986  return Tensor<Q>();
1987  }
1988  }
1989 
1990 
1991  /// Transforms one dimension of the tensor t by the matrix c, returns new contiguous tensor
1992 
1993  /// \ingroup tensor
1994  /// \code
1995  /// transform_dir(t,c,1) = r(i,j,k,...) = sum(j') t(i,j',k,...) * c(j',j)
1996  /// \endcode
1997  /// @param[in] t Tensor to transform (size of dimension to be transformed must match size of first dimension of \c c )
1998  /// @param[in] c Matrix used for the transformation
1999  /// @param[in] axis Dimension (or axis) to be transformed
2000  /// @result Returns a new, contiguous tensor
2001  template <class T, class Q>
2003  if (axis == 0) {
2004  return inner(c,t,0,axis);
2005  }
2006  else if (axis == t.ndim()-1) {
2007  return inner(t,c,axis,0);
2008  }
2009  else {
2010  return copy(inner(t,c,axis,0).cycledim(1,axis, -1)); // Copy to make contiguous
2011  }
2012  }
2013 
2014  /// Returns a new deep copy of the transpose of the input tensor
2015 
2016  /// \ingroup tensor
2017  template <class T>
2019  TENSOR_ASSERT(t.ndim() == 2, "transpose requires a matrix", t.ndim(), &t);
2020  return copy(t.swapdim(0,1));
2021  }
2022 
2023  /// Returns a new deep copy of the complex conjugate transpose of the input tensor
2024 
2025  /// \ingroup tensor
2026  template <class T>
2028  TENSOR_ASSERT(t.ndim() == 2, "conj_transpose requires a matrix", t.ndim(), &t);
2029  return conj(t.swapdim(0,1));
2030  }
2031 
2032  /// Indexing a non-constant tensor with slices returns a SliceTensor
2033 
2034  /// \ingroup tensor
2035  /// A slice tensor differs from a tensor only in that assignment
2036  /// causes the data to be copied rather than entire new copy
2037  /// generated. You will usually not instantiate one except as a
2038  /// temporary produced by indexing a tensor with slice and then
2039  /// assigning it back to a tensor, or performing some other
2040  /// operation and discarding.
2041  template <class T> class SliceTensor : public Tensor<T> {
2042  private:
2044 
2045  public:
2046 
2047  // delegating constructor
2048  SliceTensor(const Tensor<T>& t, const std::array<Slice,TENSOR_MAXDIM> s)
2049  : SliceTensor(t,s.data()) {}
2050 
2051 
2052  SliceTensor(const Tensor<T>& t, const Slice s[])
2053  : Tensor<T>(const_cast<Tensor<T>&>(t)) //!!!!!!!!!!!
2054  {
2055  // C++ standard says class derived from parameterized base class cannot
2056  // directly access the base class elements ... must explicitly reference.
2057 
2058  long nd = 0, size=1;
2059  for (long i=0; i<t._ndim; ++i) {
2060  long start=s[i].start, end=s[i].end, step=s[i].step;
2061  //std::printf("%ld input start=%ld end=%ld step=%ld\n",
2062  //i, start, end, step);
2063  if (start < 0) start += this->_dim[i];
2064  if (end < 0) end += this->_dim[i];
2065  long len = end-start+1;
2066  if (step) len /= step; // Rounds len towards zero
2067 
2068  // if input length is not exact multiple of step, round end towards start
2069  // for the same behaviour of for (i=start; i<=end; i+=step);
2070  end = start + (len-1)*step;
2071 
2072  //std::printf("%ld munged start=%ld end=%ld step=%ld len=%ld _dim=%ld\n",
2073  // i, start, end, step, len, this->_dim[i]);
2074 
2075  TENSOR_ASSERT(start>=0 && start<this->_dim[i],"slice start invalid",start,this);
2076  TENSOR_ASSERT(end>=0 && end<this->_dim[i],"slice end invalid",end,this);
2077  TENSOR_ASSERT(len>0,"slice length must be non-zero",len,this);
2078 
2079  this->_p += start * t._stride[i];
2080 
2081  if (step) {
2082  size *= len;
2083  this->_dim[nd] = len;
2084  this->_stride[nd] = step * t._stride[i];
2085  ++nd;
2086  }
2087  }
2088  //For Python interface need to be able to return a scalar inside a tensor with nd=0
2089  //TENSOR_ASSERT(nd>0,"slicing produced a scalar, but cannot return one",nd,this);
2090  for (long i=nd; i<TENSOR_MAXDIM; ++i) { // So can iterate over missing dimensions
2091  this->_dim[i] = 1;
2092  this->_stride[i] = 0;
2093  }
2094 
2095  this->_ndim = nd;
2096  this->_size = size;
2097  }
2098 
2100  BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 = (T)(*_p1));
2101  return *this;
2102  }
2103 
2104  template <class Q>
2106  BINARY_OPTIMIZED_ITERATOR(T, (*this), const Q, t, *_p0 = (T)(*_p1));
2107  return *this;
2108  }
2109 
2111  BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 = (T)(*_p1));
2112  return *this;
2113  }
2114 
2115  template <class Q>
2117  BINARY_OPTIMIZED_ITERATOR(T, (*this), const Q, t, *_p0 = (T)(*_p1));
2118  return *this;
2119  }
2120 
2122  UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 = t);
2123  return *this;
2124  }
2125 
2126  virtual ~SliceTensor() {}; // Tensor<T> destructor does enough
2127  };
2128 
2129 
2130  // Specializations for complex types
2131  template<> float_complex Tensor<float_complex>::min(long* ind) const ;
2132  template<> double_complex Tensor<double_complex>::min(long* ind) const ;
2133  template<> float_complex Tensor<float_complex>::max(long* ind) const ;
2134  template<> double_complex Tensor<double_complex>::max(long* ind) const ;
2135 
2136  // Stream stuff
2137 
2138  /// Print (for human consumption) a tensor to the stream
2139 
2140  /// \ingroup tensor
2141  template <class T>
2142  std::ostream& operator << (std::ostream& s, const Tensor<T>& t) {
2143  if (t.size() == 0) {
2144  s << "[empty tensor]\n";
2145  return s;
2146  }
2147 
2148  long maxdim = 0;
2149  long index_width = 0;
2150  for (int i = 0; i<(t.ndim()-1); ++i) {
2151  if (maxdim < t.dim(i)) maxdim = t.dim(i);
2152  }
2153  if (maxdim < 10)
2154  index_width = 1;
2155  else if (maxdim < 100)
2156  index_width = 2;
2157  else if (maxdim < 1000)
2158  index_width = 3;
2159  else if (maxdim < 10000)
2160  index_width = 4;
2161  else
2162  index_width = 6;
2163 
2164  std::ios::fmtflags oldflags = s.setf(std::ios::scientific);
2165  long oldprec = s.precision();
2166  long oldwidth = s.width();
2167 
2168  // C++ formatted IO is worse than Fortran !!
2169  for (TensorIterator<T> iter=t.unary_iterator(1,false,false); iter!=t.end(); ++iter) {
2170  const T* p = iter._p0;
2171  long inc = iter._s0;
2172  long dimj = iter.dimj;
2173  s.unsetf(std::ios::scientific);
2174  s << '[';
2175  for (long i=0; i<iter.ndim; ++i) {
2176  s.width(index_width);
2177  s << iter.ind[i];
2178  //if (i != iter.ndim)
2179  s << ",";
2180  }
2181  s << "*]";
2182 //flo s.setf(std::ios::scientific);
2183  s.setf(std::ios::fixed);
2184  for (long j=0; j<dimj; ++j, p+=inc) {
2185 //flo s.precision(4);
2186  s << " ";
2187  s.precision(8);
2188  s.width(12);
2189  s << *p;
2190  }
2191  s.unsetf(std::ios::scientific);
2192  s << std::endl;
2193  }
2194  s.setf(oldflags,std::ios::floatfield);
2195  s.precision(oldprec);
2196  s.width(oldwidth);
2197 
2198  return s;
2199  }
2200 
2201 
2202  /// Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...)
2203 
2204  /// \ingroup tensor
2205  template <class T>
2206  Tensor<T> outer(const Tensor<T>& left, const Tensor<T>& right) {
2207  long nd = left.ndim() + right.ndim();
2208  TENSOR_ASSERT(nd <= TENSOR_MAXDIM, "too many dimensions in result",
2209  nd, 0);
2210  long d[TENSOR_MAXDIM];
2211  for (long i = 0; i < left.ndim(); ++i) d[i] = left.dim(i);
2212  for (long i = 0; i < right.ndim(); ++i) d[i + left.ndim()] = right.dim(i);
2213  Tensor<T> result(nd, d, false);
2214  outer_result(left,right,result);
2215  return result;
2216  }
2217 
2218  /// Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...)
2219 
2220  /// accumulate into result, no allocation is performed
2221  template<class T>
2222  void outer_result(const Tensor<T>& left, const Tensor<T>& right, Tensor<T>& result) {
2223  TENSOR_ASSERT(left.ndim() + right.ndim() == result.ndim(),"inconsisten dimension in outer_resultn",
2224  result.ndim(),0);
2225  T *ptr = result.ptr();
2226  TensorIterator<T> iter=right.unary_iterator(1,false,true);
2227  for (TensorIterator<T> p=left.unary_iterator(); p!=left.end(); ++p) {
2228  T val1 = *p;
2229  // Cannot reorder dimensions, but can fuse contiguous dimensions
2230  for (iter.reset(); iter._p0; ++iter) {
2231  long dimj = iter.dimj;
2232  T* _p0 = iter._p0;
2233  long Tstride = iter._s0;
2234  for (long _j=0; _j<dimj; ++_j, _p0+=Tstride) {
2235  *ptr++ = val1 * (*_p0);
2236  }
2237  }
2238  }
2239  }
2240 
2241 
2242  /// Inner product ... result(i,j,...,p,q,...) = sum(z) left(i,j,...,z)*right(z,p,q,...)
2243 
2244  /// \ingroup tensor
2245  /// By default it contracts the last dimension of the left tensor and
2246  /// the first dimension of the right tensor. These defaults can be
2247  /// changed by specifying \c k0 and \c k1 , the index to contract in
2248  /// the left and right side tensors, respectively. The defaults
2249  /// correspond to (\c k0=-1 and \c k1=0 ).
2250  template <class T, class Q>
2251  Tensor<TENSOR_RESULT_TYPE(T,Q)> inner(const Tensor<T>& left, const Tensor<Q>& right,
2252  long k0=-1, long k1=0) {
2253  if (k0 < 0) k0 += left.ndim();
2254  if (k1 < 0) k1 += right.ndim();
2255  long nd = left.ndim() + right.ndim() - 2;
2256  TENSOR_ASSERT(nd!=0, "result is a scalar but cannot return one ... use dot",
2257  nd, &left);
2258 
2259 
2260  TENSOR_ASSERT(left.dim(k0) == right.dim(k1),"common index must be same length",
2261  right.dim(k1), &left);
2262 
2263  TENSOR_ASSERT(nd > 0 && nd <= TENSOR_MAXDIM,
2264  "invalid number of dimensions in the result", nd,0);
2265 
2266  long d[TENSOR_MAXDIM];
2267 
2268  long base=0;
2269  for (long i=0; i<k0; ++i) d[i] = left.dim(i);
2270  for (long i=k0+1; i<left.ndim(); ++i) d[i-1] = left.dim(i);
2271  base = left.ndim()-1;
2272  for (long i=0; i<k1; ++i) d[i+base] = right.dim(i);
2273  base--;
2274  for (long i=k1+1; i<right.ndim(); ++i) d[i+base] = right.dim(i);
2275 
2276  Tensor<TENSOR_RESULT_TYPE(T,Q)> result(nd,d);
2277 
2278  inner_result(left,right,k0,k1,result);
2279 
2280  return result;
2281  }
2282 
2283  /// Accumulate inner product into user provided, contiguous, correctly sized result tensor
2284 
2285  /// \ingroup tensor
2286  /// This routine may be used to optimize away the tensor constructor
2287  /// of the result tensor in inner loops when the result tensor may be
2288  /// reused or accumulated into. If the user calls this routine
2289  /// directly very little checking is done since it is intended as an
2290  /// optimization for small tensors. As far as the result goes, the
2291  /// caller is completely responsible for providing a contiguous tensor
2292  /// that has the correct dimensions and is appropriately initialized.
2293  /// The inner product is accumulated into result.
2294  template <class T, class Q>
2295  void inner_result(const Tensor<T>& left, const Tensor<Q>& right,
2296  long k0, long k1, Tensor< TENSOR_RESULT_TYPE(T,Q) >& result) {
2297 
2298  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
2299  // Need to include explicit optimizations for common special cases
2300  // E.g., contiguous, matrix-matrix, and 3d-tensor*matrix
2301 
2302  resultT* ptr = result.ptr();
2303 
2304  if (k0 < 0) k0 += left.ndim();
2305  if (k1 < 0) k1 += right.ndim();
2306 
2307  if (left.iscontiguous() && right.iscontiguous()) {
2308  if (k0==0 && k1==0) {
2309  // c[i,j] = a[k,i]*b[k,j] ... collapsing extra indices to i & j
2310  long dimk = left.dim(k0);
2311  long dimj = right.stride(0);
2312  long dimi = left.stride(0);
2313  mTxm(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
2314  return;
2315  }
2316  else if (k0==(left.ndim()-1) && k1==(right.ndim()-1)) {
2317  // c[i,j] = a[i,k]*b[j,k] ... collapsing extra indices to i & j
2318  long dimk = left.dim(k0);
2319  long dimi = left.size()/dimk;
2320  long dimj = right.size()/dimk;
2321  mxmT(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
2322  return;
2323  }
2324  else if (k0==0 && k1==(right.ndim()-1)) {
2325  // c[i,j] = a[k,i]*b[j,k] ... collapsing extra indices to i & j
2326  long dimk = left.dim(k0);
2327  long dimi = left.stride(0);
2328  long dimj = right.size()/dimk;
2329  mTxmT(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
2330  return;
2331  }
2332  else if (k0==(left.ndim()-1) && k1==0) {
2333  // c[i,j] = a[i,k]*b[k,j] ... collapsing extra indices to i & j
2334  long dimk = left.dim(k0);
2335  long dimi = left.size()/dimk;
2336  long dimj = right.stride(0);
2337  mxm(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
2338  return;
2339  }
2340  }
2341 
2342  long dimj = left.dim(k0);
2343  TensorIterator<Q> iter1=right.unary_iterator(1,false,false,k1);
2344 
2345  for (TensorIterator<T> iter0=left.unary_iterator(1,false,false,k0);
2346  iter0._p0; ++iter0) {
2347  T* MADNESS_RESTRICT xp0 = iter0._p0;
2348  long s0 = iter0._s0;
2349  for (iter1.reset(); iter1._p0; ++iter1) {
2350  T* MADNESS_RESTRICT p0 = xp0;
2351  Q* MADNESS_RESTRICT p1 = iter1._p0;
2352  long s1 = iter1._s0;
2353  resultT sum = 0;
2354  for (long j=0; j<dimj; ++j,p0+=s0,p1+=s1) {
2355  sum += (*p0) * (*p1);
2356  }
2357  *ptr++ += sum;
2358  }
2359  }
2360  }
2361 
2362  /// Transform all dimensions of the tensor t by the matrix c
2363 
2364  /// \ingroup tensor
2365  /// Often used to transform all dimensions from one basis to another
2366  /// \code
2367  /// result(i,j,k...) <-- sum(i',j', k',...) t(i',j',k',...) c(i',i) c(j',j) c(k',k) ...
2368  /// \endcode
2369  /// The input dimensions of \c t must all be the same and agree with
2370  /// the first dimension of \c c . The dimensions of \c c may differ in
2371  /// size. If the dimensions of \c c are the same, and the operation
2372  /// is being performed repeatedly, then you might consider calling \c
2373  /// fast_transform instead which enables additional optimizations and
2374  /// can eliminate all constructor overhead and improve cache locality.
2375  ///
2376  template <class T, class Q>
2378  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
2379  TENSOR_ASSERT(c.ndim() == 2,"second argument must be a matrix",c.ndim(),&c);
2380  if (c.dim(0)==c.dim(1) && t.iscontiguous() && c.iscontiguous()) {
2381  Tensor<resultT> result(t.ndim(),t.dims(),false);
2382  Tensor<resultT> work(t.ndim(),t.dims(),false);
2383  return fast_transform(t, c, result, work);
2384  }
2385  else {
2386  Tensor<resultT> result = t;
2387  for (long i=0; i<t.ndim(); ++i) {
2388  result = inner(result,c,0,0);
2389  }
2390  return result;
2391  }
2392  }
2393 
2394  /// Transform all dimensions of the tensor t by distinct matrices c
2395 
2396  /// \ingroup tensor
2397  /// Similar to transform but each dimension is transformed with a
2398  /// distinct matrix.
2399  /// \code
2400  /// result(i,j,k...) <-- sum(i',j', k',...) t(i',j',k',...) c[0](i',i) c[1](j',j) c[2](k',k) ...
2401  /// \endcode
2402  /// The first dimension of the matrices c must match the corresponding
2403  /// dimension of t.
2404  template <class T, class Q>
2406  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
2407  Tensor<resultT> result = t;
2408  for (long i=0; i<t.ndim(); ++i) {
2409  result = inner(result,c[i],0,0);
2410  }
2411  return result;
2412  }
2413 
2414  /// Restricted but heavily optimized form of transform()
2415 
2416  /// \ingroup tensor
2417  /// Both dimensions of \c c must be the same and match all dimensions
2418  /// of the input tensor \c t. All tensors must be contiguous.
2419  ///
2420  /// Performs the same operation as \c transform but it requires
2421  /// that the caller pass in workspace and a preallocated result,
2422  /// hoping that that both can be reused. If the result and
2423  /// workspace are reused between calls, then no tensor
2424  /// constructors need be called and cache locality should be
2425  /// improved. By passing in the workspace, this routine is kept
2426  /// thread safe.
2427  ///
2428  /// The input, result and workspace tensors must be distinct.
2429  ///
2430  /// All input tensors must be contiguous and fastest execution
2431  /// will result if all dimensions are approriately aligned and
2432  /// multiples of the underlying vector length. The workspace and
2433  /// the result must be of the same size as the input \c t . The
2434  /// result tensor need not be initialized before calling
2435  /// fast_transform.
2436  ///
2437  /// \code
2438  /// result(i,j,k,...) <-- sum(i',j', k',...) t(i',j',k',...) c(i',i) c(j',j) c(k',k) ...
2439  /// \endcode
2440  ///
2441  /// The input dimensions of \c t must all be the same .
2442  template <class T, class Q>
2444  Tensor< TENSOR_RESULT_TYPE(T,Q) >& workspace) {
2445  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
2446  const Q *pc=c.ptr();
2447  resultT *t0=workspace.ptr(), *t1=result.ptr();
2448  if (t.ndim()&1) {
2449  t0 = result.ptr();
2450  t1 = workspace.ptr();
2451  }
2452 
2453  long dimj = c.dim(1);
2454  long dimi = 1;
2455  for (int n=1; n<t.ndim(); ++n) dimi *= dimj;
2456 
2457 #if HAVE_IBMBGQ
2458  long nij = dimi*dimj;
2459  if (IS_UNALIGNED(pc) || IS_UNALIGNED(t0) || IS_UNALIGNED(t1)) {
2460  for (long i=0; i<nij; ++i) t0[i] = 0.0;
2461  mTxm(dimi, dimj, dimj, t0, t.ptr(), pc);
2462  for (int n=1; n<t.ndim(); ++n) {
2463  for (long i=0; i<nij; ++i) t1[i] = 0.0;
2464  mTxm(dimi, dimj, dimj, t1, t0, pc);
2465  std::swap(t0,t1);
2466  }
2467  }
2468  else {
2469  mTxmq_padding(dimi, dimj, dimj, dimj, t0, t.ptr(), pc);
2470  for (int n=1; n<t.ndim(); ++n) {
2471  mTxmq_padding(dimi, dimj, dimj, dimj, t1, t0, pc);
2472  std::swap(t0,t1);
2473  }
2474  }
2475 #else
2476  // Now assume no restriction on the use of mtxmq
2477  mTxmq(dimi, dimj, dimj, t0, t.ptr(), pc);
2478  for (int n=1; n<t.ndim(); ++n) {
2479  mTxmq(dimi, dimj, dimj, t1, t0, pc);
2480  std::swap(t0,t1);
2481  }
2482 #endif
2483 
2484  return result;
2485  }
2486 
2487  /// Return a new tensor holding the absolute value of each element of t
2488 
2489  /// \ingroup tensor
2490  template <class T>
2492  typedef typename Tensor<T>::scalar_type scalar_type;
2493  Tensor<scalar_type> result(t.ndim(),t.dims(),false);
2494  BINARY_OPTIMIZED_ITERATOR(scalar_type,result,const T,t,*_p0 = std::abs(*_p1));
2495  return result;
2496  }
2497 
2498  /// Return a new tensor holding the argument of each element of t (complex types only)
2499 
2500  /// \ingroup tensor
2501  template <class T>
2503  typedef typename Tensor<T>::scalar_type scalar_type;
2504  Tensor<scalar_type> result(t.ndim(),t.dims(),false);
2505  BINARY_OPTIMIZED_ITERATOR(scalar_type,result,T,t,*_p0 = std::arg(*_p1));
2506  return result;
2507  }
2508 
2509  /// Return a new tensor holding the real part of each element of t (complex types only)
2510 
2511  /// \ingroup tensor
2512  template <class T>
2514  typedef typename Tensor<T>::scalar_type scalar_type;
2515  Tensor<scalar_type> result(t.ndim(),t.dims(),false);
2516  BINARY_OPTIMIZED_ITERATOR(scalar_type,result,const T,t,*_p0 = std::real(*_p1));
2517  return result;
2518  }
2519 
2520  /// Return a new tensor holding the imaginary part of each element of t (complex types only)
2521 
2522  /// \ingroup tensor
2523  template <class T>
2525  typedef typename Tensor<T>::scalar_type scalar_type;
2526  Tensor<scalar_type> result(t.ndim(),t.dims(),false);
2527  BINARY_OPTIMIZED_ITERATOR(scalar_type,result,const T,t,*_p0 = std::imag(*_p1));
2528  return result;
2529  }
2530 
2531  /// Returns a new deep copy of the complex conjugate of the input tensor (complex types only)
2532 
2533  /// \ingroup tensor
2534  template <class T>
2536  Tensor<T> result(t.ndim(),t.dims(),false);
2537  BINARY_OPTIMIZED_ITERATOR(T,result,const T,t,*_p0 = conditional_conj(*_p1));
2538  return result;
2539  }
2540 }
2541 
2542 #undef TENSOR_SHARED_PTR
2543 
2544 #endif // MADNESS_TENSOR_TENSOR_H__INCLUDED
double q(double t)
Definition: DKops.h:18
Provides routines for internal use optimized for aligned data.
Interface templates for the archives (serialization).
Declares BaseTensor.
std::complex< double > double_complex
Definition: cfft.h:14
The base class for tensors defines generic capabilities.
Definition: basetensor.h:85
bool conforms(const BaseTensor *t) const
Returns true if this and *t are the same shape and size.
Definition: basetensor.h:159
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
bool iscontiguous() const
Returns true if the tensor refers to contiguous memory locations.
Definition: basetensor.h:168
void mapdim_inplace(const std::vector< long > &map)
General permutation of dimensions.
Definition: basetensor.cc:156
long _stride[TENSOR_MAXDIM]
Increment between elements in each dimension.
Definition: basetensor.h:97
long stride(int i) const
Returns the stride associated with dimension i.
Definition: basetensor.h:150
long _size
Number of elements in the tensor.
Definition: basetensor.h:93
long id() const
Returns the typeid of the tensor (c.f., TensorTypeData<T> )
Definition: basetensor.h:141
void set_dims_and_size(long nd, const long d[])
Definition: basetensor.h:99
long _id
Id from TensorTypeData<T> in type_data.h.
Definition: basetensor.h:95
const long * dims() const
Returns the array of tensor dimensions.
Definition: basetensor.h:153
void splitdim_inplace(long i, long dimi0, long dimi1)
Splits dimension i.
Definition: basetensor.cc:88
void swapdim_inplace(long i, long j)
Swaps the dimensions.
Definition: basetensor.cc:124
void fusedim_inplace(long i)
Fuses dimensions i and i+1.
Definition: basetensor.cc:107
long _dim[TENSOR_MAXDIM]
Size of each dimension.
Definition: basetensor.h:96
long ndim() const
Returns the number of dimensions in the tensor.
Definition: basetensor.h:144
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
void cycledim_inplace(long shift, long start, long end)
Cyclic shift of dimensions.
Definition: basetensor.cc:134
long _ndim
Number of dimensions (-1=invalid; 0=no supported; >0=tensor)
Definition: basetensor.h:94
void reshape_inplace(const std::vector< long > &d)
Reshapes the tensor inplace.
Definition: basetensor.cc:76
Indexing a non-constant tensor with slices returns a SliceTensor.
Definition: tensor.h:2041
virtual ~SliceTensor()
Definition: tensor.h:2126
SliceTensor(const Tensor< T > &t, const Slice s[])
Definition: tensor.h:2052
SliceTensor< T > & operator=(const SliceTensor< Q > &t)
Definition: tensor.h:2105
SliceTensor< T > & operator=(const Tensor< T > &t)
Definition: tensor.h:2110
SliceTensor< T > & operator=(const T &t)
Definition: tensor.h:2121
SliceTensor(const Tensor< T > &t, const std::array< Slice, TENSOR_MAXDIM > s)
Definition: tensor.h:2048
SliceTensor< T > & operator=(const Tensor< Q > &t)
Definition: tensor.h:2116
SliceTensor< T > & operator=(const SliceTensor< T > &t)
Definition: tensor.h:2099
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
Definition: tensoriter.h:61
long dimj
Definition: tensoriter.h:70
long _s0
Definition: tensoriter.h:71
void reset()
Reset the iterator back to the start ...
Definition: tensoriter.h:354
T * _p0
Definition: tensoriter.h:66
Traits class to specify support of numeric types.
Definition: type_data.h:56
A tensor is a multidimension array.
Definition: tensor.h:317
scalar_type absmax(long *ind=0) const
Return the absolute maximum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1754
Tensor< T > reshape(const std::vector< long > &d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition: tensor.h:1405
const T & operator()(long i, long j, long k, long l) const
4-d indexing operation without bounds checking.
Definition: tensor.h:901
T & operator()(const std::vector< long > ind)
General indexing operation with bounds checking.
Definition: tensor.h:1032
void deallocate()
Definition: tensor.h:397
const Tensor< T > operator()(const Slice &s0) const
Return a 1d SliceTensor that views the specified range of the 1d Tensor.
Definition: tensor.h:1105
T *MADNESS_RESTRICT _p
Definition: tensor.h:321
const Tensor< T > reshape(long dim0, long dim1, long dim2) const
Returns new view/tensor rehaping to conforming 3-d tensor with given dimensions.
Definition: tensor.h:1471
Tensor< T > reshape(long dim0, long dim1, long dim2)
Returns new view/tensor rehaping to conforming 3-d tensor with given dimensions.
Definition: tensor.h:1460
Tensor(const std::vector< long > &d, bool dozero=true)
Create and optionally zero new n-d tensor. This is the most general constructor.
Definition: tensor.h:536
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
IsSupported< TensorTypeData< Q >, Tensor< TENSOR_RESULT_TYPE(T, Q)> >::type operator*(const Q &x) const
Multiplication of tensor by a scalar of a supported type to produce a new tensor.
Definition: tensor.h:617
Tensor(long d0, long d1)
Create and zero new 2-d tensor.
Definition: tensor.h:481
SliceTensor< T > operator()(long i, const Slice &s1, const Slice &s2)
Return a 2d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1195
Tensor< T > fusedim(long i)
Returns new view/tensor fusing contiguous dimensions i and i+1.
Definition: tensor.h:1587
Tensor< TENSOR_RESULT_TYPE(T, Q) > operator+(const Tensor< Q > &t) const
Addition of two tensors to produce a new tensor.
Definition: tensor.h:592
const T & operator()(const long ind[]) const
Politically incorrect general indexing operation without bounds checking.
Definition: tensor.h:1016
IsSupported< TensorTypeData< Q >, Tensor< TENSOR_RESULT_TYPE(T, Q)> >::type operator/(const Q &x) const
Divide tensor by a scalar of a supported type to produce a new tensor.
Definition: tensor.h:630
scalar_type absmin(long *ind=0) const
Return the absolute minimum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1733
const Tensor< T > fusedim(long i) const
Returns new view/tensor fusing contiguous dimensions i and i+1.
Definition: tensor.h:1596
SliceTensor< T > operator()(const std::array< Slice, TENSOR_MAXDIM > &s)
General slicing operation.
Definition: tensor.h:1080
const Tensor< T > swapdim(long idim, long jdim) const
Returns new view/tensor swaping dimensions i and j.
Definition: tensor.h:1614
const Tensor< T > mapdim(const std::vector< long > &map) const
Returns new view/tensor permuting the dimensions.
Definition: tensor.h:1633
Tensor< T > mapdim(const std::vector< long > &map)
Returns new view/tensor permuting the dimensions.
Definition: tensor.h:1624
Tensor< T > splitdim(long i, long dimi0, long dimi1)
Returns new view/tensor splitting dimension i as dimi0*dimi1 to produce conforming d+1 dimension tens...
Definition: tensor.h:1569
const Tensor< T > operator()(long i, const Slice &s1, long k) const
Return a 1d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1285
T sum() const
Returns the sum of all elements of the tensor.
Definition: tensor.h:1662
const Tensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2) const
Return a 3d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1185
SliceTensor< T > operator()(const Slice &s0, long j, long k)
Return a 1d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1295
BaseTensor * base()
Returns a pointer to the base class.
Definition: tensor.h:1834
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition: tensor.h:409
Tensor(long d0, long d1, long d2, long d3, long d4, long d5)
Create and zero new 6-d tensor.
Definition: tensor.h:527
T type
C++ typename of this tensor.
Definition: tensor.h:406
SliceTensor< T > operator()(long i, const Slice &s1)
Return a 1d SliceTensor that views the specified range of the 2d Tensor.
Definition: tensor.h:1115
SliceTensor< T > operator()(const std::vector< Slice > &s)
General slicing operation.
Definition: tensor.h:1060
SliceTensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2, const Slice &s3, const Slice &s4, const Slice &s5)
Return a 1-6d SliceTensor that views the specified range of the 6d Tensor.
Definition: tensor.h:1359
SliceTensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2, const Slice &s3)
Return a 1-4d SliceTensor that views the specified range of the 4d Tensor.
Definition: tensor.h:1315
Tensor< T > reshape(long dim0, long dim1, long dim2, long dim3)
Returns new view/tensor rehaping to conforming 4-d tensor with given dimensions.
Definition: tensor.h:1483
bool conforms(const Tensor< Q > &t) const
Test if *this and t conform.
Definition: tensor.h:1657
const BaseTensor * base() const
Returns a pointer to the base class.
Definition: tensor.h:1839
Tensor< T > & fill(T x)
Inplace fill with a scalar (legacy name)
Definition: tensor.h:562
T & operator()(long i, long j)
2-d indexing operation without bounds checking.
Definition: tensor.h:825
TensorIterator< T, Q > binary_iterator(const Tensor< Q > &q, long iterlevel=0, bool optimize=true, bool fusedim=true, long jdim=default_jdim) const
Return iterator over two tensors.
Definition: tensor.h:1854
const Tensor< T > reshape(long dim0, long dim1, long dim2, long dim3, long dim4, long dim5) const
Returns new view/tensor rehaping to conforming 6-d tensor with given dimensions.
Definition: tensor.h:1549
IsSupported< TensorTypeData< Q >, Tensor< TENSOR_RESULT_TYPE(T, Q)> >::type operator+(const Q &x) const
Add a scalar of the same type to all elements of a tensor producing a new tensor.
Definition: tensor.h:643
const Tensor< T > operator()(const Slice &s0, long j, const Slice &s2) const
Return a 2d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1225
const Tensor< T > operator()(const Slice &s0, long j) const
Return a 1d constant Tensor that views the specified range of the 2d Tensor.
Definition: tensor.h:1145
const Tensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2, const Slice &s3, const Slice &s4) const
Return a 1-5d constant Tensor that views the specified range of the 5d Tensor.
Definition: tensor.h:1348
Tensor< T > & operator=(const Tensor< T > &t)
Assignment is shallow (same as copy constructor)
Definition: tensor.h:443
const Tensor< T > operator()(const std::vector< Slice > &s) const
General slicing operation (const)
Definition: tensor.h:1070
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition: tensor.h:1384
SliceTensor< T > operator()(const Slice &s0, const Slice &s1, long k)
Return a 2d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1235
Tensor< T > & gaxpy(T alpha, const Tensor< T > &t, T beta)
Inplace generalized saxpy ... this = this*alpha + other*beta.
Definition: tensor.h:1804
const T & operator()(long i, long j, long k, long l, long m) const
5-d indexing operation without bounds checking.
Definition: tensor.h:940
Tensor< T > flat()
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition: tensor.h:1555
TensorIterator< T > unary_iterator(long iterlevel=0, bool optimize=true, bool fusedim=true, long jdim=default_jdim) const
Return iterator over single tensor.
Definition: tensor.h:1844
TensorTypeData< T >::float_scalar_type float_scalar_type
C++ typename of the floating point type associated with scalar real type.
Definition: tensor.h:412
Tensor< T > & conj()
Inplace complex conjugate.
Definition: tensor.h:716
Tensor< T > & operator=(T x)
Inplace fill tensor with scalar.
Definition: tensor.h:553
Tensor(const Tensor< T > &t)
Copy constructor is shallow (same as assignment)
Definition: tensor.h:428
const Tensor< T > cycledim(long nshift, long start, long end) const
Returns new view/tensor cycling the sub-dimensions (start,...,end) with shift steps.
Definition: tensor.h:1649
SliceTensor< T > operator()(const Slice &s0, long j, const Slice &s2)
Return a 2d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1215
const Tensor< T > operator()(const Slice &s0, const Slice &s1) const
Return a 2d constant Tensor that views the specified range of the 2d Tensor.
Definition: tensor.h:1165
Tensor< T > & operator+=(const Tensor< Q > &t)
Inplace addition of two tensors.
Definition: tensor.h:572
SliceTensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2)
Return a 3d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1175
T trace(const Tensor< T > &t) const
Return the trace of two tensors (no complex conjugate invoked)
Definition: tensor.h:1776
T & operator()(const long ind[])
Politically incorrect general indexing operation without bounds checking.
Definition: tensor.h:1000
T & operator[](long i)
1-d indexing operation using [] without bounds checking.
Definition: tensor.h:780
Tensor< T > cycledim(long nshift, long start, long end)
Returns new view/tensor cycling the sub-dimensions (start,...,end) with shift steps.
Definition: tensor.h:1641
Tensor< T > & fillindex()
Inplace fill with the index of each element.
Definition: tensor.h:748
const Tensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2, const Slice &s3, const Slice &s4, const Slice &s5) const
Return a 1-6d constant Tensor that views the specified range of the 6d Tensor.
Definition: tensor.h:1371
const Tensor< T > operator()(const Slice &s0, const Slice &s1, long k) const
Return a 2d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1245
T max(long *ind=0) const
Return the maximum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1703
Tensor(long d0)
Create and zero new 1-d tensor.
Definition: tensor.h:472
T product() const
Return the product of all elements of the tensor.
Definition: tensor.h:1676
SliceTensor< T > operator()(long i, const Slice &s1, long k)
Return a 1d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1275
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type operator*=(const Q &x)
Inplace multiplication by scalar of supported type.
Definition: tensor.h:675
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type operator+=(const Q &x)
Inplace increment by scalar of supported type.
Definition: tensor.h:696
const Tensor< T > reshape(long dim0, long dim1, long dim2, long dim3) const
Returns new view/tensor rehaping to conforming 4-d tensor with given dimensions.
Definition: tensor.h:1495
Tensor< T > & unaryop(opT &op)
Inplace apply a unary function to each element of the tensor.
Definition: tensor.h:1792
const Tensor< T > operator()(long i, const Slice &s1) const
Return a 1d SliceTensor that views the specified range of the 2d Tensor.
Definition: tensor.h:1125
T min(long *ind=0) const
Return the minimum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1683
Tensor< T > operator-() const
Unary negation producing a new tensor.
Definition: tensor.h:663
Tensor< T > reshape(long dim0, long dim1)
Returns new view/tensor rehaping to conforming 2-d tensor with given dimensions.
Definition: tensor.h:1439
virtual ~Tensor()
Definition: tensor.h:1881
const Tensor< T > reshape(const std::vector< long > &d) const
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition: tensor.h:1413
void allocate(long nd, const long d[], bool dozero)
Definition: tensor.h:324
const Tensor< T > reshape(long dim0, long dim1, long dim2, long dim3, long dim4) const
Returns new view/tensor rehaping to conforming 5-d tensor with given dimensions.
Definition: tensor.h:1521
T sumsq() const
Returns the sum of the squares of the elements.
Definition: tensor.h:1669
const Tensor< T > flat() const
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition: tensor.h:1561
Tensor< T > & emul(const Tensor< T > &t)
Inplace multiply by corresponding elements of argument Tensor.
Definition: tensor.h:1798
const Tensor< T > reshape(long dim0, long dim1) const
Returns new view/tensor rehaping to conforming 2-d tensor with given dimensions.
Definition: tensor.h:1449
Tensor(long d0, long d1, long d2, long d3)
Create and zero new 4-d tensor.
Definition: tensor.h:502
Tensor(long d0, long d1, long d2)
Create and zero new 3-d tensor.
Definition: tensor.h:491
const T & operator[](long i) const
1-d indexing operation using [] without bounds checking.
Definition: tensor.h:791
Tensor< T > & operator-=(const Tensor< Q > &t)
Inplace subtraction of two tensors.
Definition: tensor.h:582
TensorIterator< T, Q, R > ternary_iterator(const Tensor< Q > &q, const Tensor< R > &r, long iterlevel=0, bool optimize=true, bool fusedim=true, long jdim=default_jdim) const
Return iterator over three tensors.
Definition: tensor.h:1865
const T & operator()(long i) const
1-d indexing operation without bounds checking.
Definition: tensor.h:813
const Tensor< T > operator()(const Slice &s0, long j, long k) const
Return a 1d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1305
const T & operator()(long i, long j) const
2-d indexing operation without bounds checking.
Definition: tensor.h:838
bool has_data() const
Definition: tensor.h:1886
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type operator-=(const Q &x)
Inplace decrement by scalar of supported type.
Definition: tensor.h:707
T & operator()(long i, long j, long k)
3-d indexing operation without bounds checking.
Definition: tensor.h:852
const T & operator()(long i, long j, long k) const
3-d indexing operation without bounds checking.
Definition: tensor.h:867
Tensor< T > swapdim(long idim, long jdim)
Returns new view/tensor swaping dimensions i and j.
Definition: tensor.h:1605
IsSupported< TensorTypeData< Q >, Tensor< TENSOR_RESULT_TYPE(T, Q)> >::type operator-(const Q &x) const
Subtract a scalar of the same type from all elements producing a new tensor.
Definition: tensor.h:656
TENSOR_RESULT_TYPE(T, Q) trace_conj(const Tensor< Q > &t) const
Return the trace of two tensors with complex conjugate of the leftmost (i.e., this)
Definition: tensor.h:1784
Tensor< T > reshape(long dim0, long dim1, long dim2, long dim3, long dim4)
Returns new view/tensor rehaping to conforming 5-d tensor with given dimensions.
Definition: tensor.h:1508
T & operator()(long i, long j, long k, long l, long m, long n)
6-d indexing operation without bounds checking.
Definition: tensor.h:961
SliceTensor< T > operator()(const Slice &s0, long j)
Return a 1d SliceTensor that views the specified range of the 2d Tensor.
Definition: tensor.h:1135
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
const Tensor< T > splitdim(long i, long dimi0, long dimi1) const
Returns new view/tensor splitting dimension i as dimi0*dimi1 to produce conforming d+1 dimension tens...
Definition: tensor.h:1578
Tensor(long d0, long d1, long d2, long d3, long d4)
Create and zero new 5-d tensor.
Definition: tensor.h:514
const Tensor< T > operator()(long i, const Slice &s1, const Slice &s2) const
Return a 2d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1205
Tensor()
Default constructor does not allocate any data and sets ndim=-1, size=0, _p=0, and id.
Definition: tensor.h:415
const TensorIterator< T > & end() const
End point for forward iteration.
Definition: tensor.h:1876
const T & operator()(const std::vector< long > ind) const
General indexing operation with bounds checking.
Definition: tensor.h:1046
T & operator()(long i, long j, long k, long l, long m)
5-d indexing operation without bounds checking.
Definition: tensor.h:920
const Tensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2, const Slice &s3) const
Return a 1-4d constant Tensor that views the specified range of the 4d Tensor.
Definition: tensor.h:1326
SliceTensor< T > operator()(const Slice &s0)
Return a 1d SliceTensor that views the specified range of the 1d Tensor.
Definition: tensor.h:1095
static bool bounds_checking()
Return true if bounds checking was enabled at compile time.
Definition: tensor.h:768
T & operator()(long i)
1-d indexing operation without bounds checking.
Definition: tensor.h:802
void clear()
Frees all memory and resests to state of default constructor.
Definition: tensor.h:1884
T & operator()(long i, long j, long k, long l)
4-d indexing operation without bounds checking.
Definition: tensor.h:883
const T * ptr() const
Returns a pointer to the internal data.
Definition: tensor.h:1829
TENSOR_SHARED_PTR< T > _shptr
Definition: tensor.h:322
const Tensor< T > operator()(const std::array< Slice, TENSOR_MAXDIM > &s) const
General slicing operation (const)
Definition: tensor.h:1088
SliceTensor< T > operator()(long i, long j, const Slice &s2)
Return a 1d SliceTensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1255
const Tensor< T > operator()(long i, long j, const Slice &s2) const
Return a 1d constant Tensor that views the specified range of the 3d Tensor.
Definition: tensor.h:1265
Tensor(long nd, const long d[], bool dozero=true)
Politically incorrect general constructor.
Definition: tensor.h:545
Tensor< T > reshape(long dim0, long dim1, long dim2, long dim3, long dim4, long dim5)
Returns new view/tensor rehaping to conforming 6-d tensor with given dimensions.
Definition: tensor.h:1535
const Tensor< T > reshape(long dim0) const
Returns new view/tensor rehapings to conforming 1-d tensor with given dimension.
Definition: tensor.h:1429
Tensor< T > & fillrandom()
Inplace fill with random values ( [0,1] for floats, [0,MAXSIZE] for integers)
Definition: tensor.h:724
const T & operator()(long i, long j, long k, long l, long m, long n) const
6-d indexing operation without bounds checking.
Definition: tensor.h:983
SliceTensor< T > operator()(const Slice &s0, const Slice &s1, const Slice &s2, const Slice &s3, const Slice &s4)
Return a 1-5d SliceTensor that views the specified range of the 5d Tensor.
Definition: tensor.h:1337
Tensor< T > & screen(double x)
Inplace set elements of *this less than x in absolute magnitude to zero.
Definition: tensor.h:758
Tensor< T > reshape(long dim0)
Returns new view/tensor rehapings to conforming 1-d tensor with given dimension.
Definition: tensor.h:1421
Tensor< TENSOR_RESULT_TYPE(T, Q) > operator-(const Tensor< Q > &t) const
Subtraction of two tensors to produce a new tensor.
Definition: tensor.h:604
SliceTensor< T > operator()(const Slice &s0, const Slice &s1)
Return a 2d SliceTensor that views the specified range of the 2d Tensor.
Definition: tensor.h:1155
const Tensor< T > reshape(int ndimnew, const long *d) const
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition: tensor.h:1395
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition: derivatives.cc:72
const double m
Definition: gfit.cc:199
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
archive_array< T > wrap(const T *, unsigned int)
Factory function to wrap a dynamically allocated pointer as a typed archive_array.
Definition: archive.h:913
Tensor< T > conj_transpose(const Tensor< T > &t)
Returns a new deep copy of the complex conjugate transpose of the input tensor.
Definition: tensor.h:2027
Tensor< typename Tensor< T >::scalar_type > real(const Tensor< T > &t)
Return a new tensor holding the real part of each element of t (complex types only)
Definition: tensor.h:2513
Tensor< TENSOR_RESULT_TYPE(T, Q) > & fast_transform(const Tensor< T > &t, const Tensor< Q > &c, Tensor< TENSOR_RESULT_TYPE(T, Q) > &result, Tensor< TENSOR_RESULT_TYPE(T, Q) > &workspace)
Restricted but heavily optimized form of transform()
Definition: tensor.h:2443
void inner_result(const Tensor< T > &left, const Tensor< Q > &right, long k0, long k1, Tensor< TENSOR_RESULT_TYPE(T, Q) > &result)
Accumulate inner product into user provided, contiguous, correctly sized result tensor.
Definition: tensor.h:2295
Tensor< typename Tensor< T >::scalar_type > arg(const Tensor< T > &t)
Return a new tensor holding the argument of each element of t (complex types only)
Definition: tensor.h:2502
Tensor< typename Tensor< T >::scalar_type > imag(const Tensor< T > &t)
Return a new tensor holding the imaginary part of each element of t (complex types only)
Definition: tensor.h:2524
const double beta
Definition: gygi_soltion.cc:62
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
Macros and tools pertaining to the configuration of MADNESS.
#define MADNESS_PRAGMA_GCC(x)
Definition: madness_config.h:205
Internal use only.
double norm(const T &t)
Definition: adquad.h:42
T mynorm(T t)
Definition: tensor.h:260
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
std::vector< CCPairFunction< T, NDIM > > operator+(const std::vector< CCPairFunction< T, NDIM >> c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition: ccpairfunction.h:1047
double abs(double x)
Definition: complexfun.h:48
void mTxmT(long dimi, long dimj, long dimk, T *MADNESS_RESTRICT c, const T *a, const T *b)
Matrix += Matrix transpose * matrix transpose ... MKL interface version.
Definition: mxm.h:238
void outer_result(const Tensor< T > &left, const Tensor< T > &right, Tensor< T > &result)
Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...)
Definition: tensor.h:2222
Q conditional_conj(const Q &coeff)
For real types return value, for complex return conjugate.
Definition: tensor.h:255
GenTensor< TENSOR_RESULT_TYPE(R, Q)> general_transform(const GenTensor< R > &t, const Tensor< Q > c[])
Definition: gentensor.h:274
std::vector< CCPairFunction< T, NDIM > > operator-(const std::vector< CCPairFunction< T, NDIM >> c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition: ccpairfunction.h:1055
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
void mxm(long dimi, long dimj, long dimk, T *MADNESS_RESTRICT c, const T *a, const T *b)
Matrix += Matrix * matrix ... BLAS/MKL interface version.
Definition: mxm.h:199
GenTensor< TENSOR_RESULT_TYPE(R, Q)> transform_dir(const GenTensor< R > &t, const Tensor< Q > &c, const int axis)
Definition: lowranktensor.h:1099
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
void mTxm(long dimi, long dimj, long dimk, T *MADNESS_RESTRICT c, const T *a, const T *b)
Matrix += Matrix transpose * matrix ... MKL interface version.
Definition: mxm.h:212
Function< Q, NDIM > convert(const Function< T, NDIM > &f, bool fence=true)
Type conversion implies a deep copy. No communication except for optional fence.
Definition: mra.h:2032
response_space transpose(response_space &f)
Definition: basic_operators.cc:10
void mTxmq_padding(long dimi, long dimj, long dimk, long ext_b, cT *c, const aT *a, const bT *b)
Definition: mtxmq.h:96
std::ostream & operator<<(std::ostream &os, const particle< PDIM > &p)
Definition: lowrankfunction.h:397
TENSOR_RESULT_TYPE(T, R) inner(const Function< T
Computes the scalar/inner product between two functions.
Definition: vmra.h:1059
static double pop(std::vector< double > &v)
Definition: SCF.cc:113
Function< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
static void aligned_zero(long n, T *a)
Definition: aligned.h:55
std::vector< CCPairFunction< T, NDIM > > operator*(const double fac, const std::vector< CCPairFunction< T, NDIM > > &arg)
Definition: ccpairfunction.h:1084
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
double imag(double x)
Definition: complexfun.h:56
std::string type(const PairType &n)
Definition: PNOParameters.h:18
static const long default_jdim
Definition: tensoriter.h:57
double real(double x)
Definition: complexfun.h:52
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition: vmra.h:689
void mxmT(long dimi, long dimj, long dimk, T *MADNESS_RESTRICT c, const T *a, const T *b)
Matrix += Matrix * matrix transpose ... MKL interface version.
Definition: mxm.h:225
std::enable_if< std::is_base_of< ProjectorBase, projT >::value, OuterProjector< projT, projQ > >::type outer(const projT &p0, const projQ &p1)
Definition: projector.h:457
void swap(Vector< T, N > &l, Vector< T, N > &r)
Swap the contents of two Vectors.
Definition: vector.h:497
void mTxmq(long dimi, long dimj, long dimk, T *MADNESS_RESTRICT c, const T *a, const T *b, long ldb=-1)
Matrix = Matrix transpose * matrix ... MKL interface version.
Definition: mxm.h:257
Definition: mraimpl.h:50
static long abs(long a)
Definition: tensor.h:218
static const double b
Definition: nonlinschro.cc:119
static const double a
Definition: nonlinschro.cc:118
Implement dummy posix_memalign if it is missing on the system.
int posix_memalign(void **memptr, std::size_t alignment, std::size_t size)
Definition: posixmem.h:44
std::complex< float > float_complex
Definition: ran.h:39
double Q(double a)
Definition: relops.cc:20
static const double c
Definition: relops.cc:10
static const long k
Definition: rk.cc:44
Definition: test_ccpairfunction.cc:22
Definition: type_data.h:146
static void load(const Archive &s, Tensor< T > &t)
Definition: tensor.h:1913
Default load of an object via serialize(ar, t).
Definition: archive.h:666
static void store(const Archive &s, const Tensor< T > &t)
Definition: tensor.h:1898
Default store of an object via serialize(ar, t).
Definition: archive.h:611
static Q op(const Q &coeff)
Definition: tensor.h:248
For real types return value, for complex return conjugate.
Definition: tensor.h:239
static Q op(const Q &coeff)
Definition: tensor.h:240
static const double s0
Definition: tdse4.cc:83
#define TENSOR_ALIGNMENT
#define IS_UNALIGNED(p)
Definition: tensor.h:234
#define UNARY_UNOPTIMIZED_ITERATOR(X, x, exp)
Definition: tensor_macros.h:678
#define BINARY_OPTIMIZED_ITERATOR(X, x, Y, y, exp)
Definition: tensor_macros.h:701
#define UNARY_OPTIMIZED_ITERATOR(X, x, exp)
Definition: tensor_macros.h:658
#define TERNARY_OPTIMIZED_ITERATOR(X, x, Y, y, Z, z, exp)
Definition: tensor_macros.h:719
#define TENSOR_MAXDIM
Definition: tensor_macros.h:194
Declares and implements TensorException.
#define TENSOR_ASSERT(condition, msg, value, t)
Definition: tensorexcept.h:130
#define TENSOR_EXCEPTION(msg, value, t)
Definition: tensorexcept.h:126
Declares TensorIterator.
const double alpha
Definition: test_coulomb.cc:54
void d()
Definition: test_sig.cc:79
const double offset
Definition: testfuns.cc:143
std::size_t axis
Definition: testpdiff.cc:59
double k0
Definition: testperiodic.cc:66
double k1
Definition: testperiodic.cc:67