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
36#include <madness/misc/ran.h>
38
39#include <memory>
40#include <complex>
41#include <vector>
42#include <array>
43#include <cmath>
44#include <cstdlib>
45#include <cstddef>
46
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>
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
217namespace std {
218 static long abs(long a) {
219 return a>=0 ? a : -a;
220 }
221}
222#else
223namespace std {
224 static long abs(long a) {
225 return std::labs(a);
226 }
227}
228#endif
229#endif
230
231
232namespace 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>
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>
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
371MADNESS_PRAGMA_GCC(diagnostic push)
372MADNESS_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]);
376MADNESS_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;
449MADNESS_PRAGMA_GCC(diagnostic push)
450MADNESS_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 }
455MADNESS_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
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
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)
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;
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
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
const long * dims() const
Returns the array of tensor dimensions.
Definition basetensor.h:153
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
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< Q > &t)
Definition tensor.h:2116
SliceTensor< T > & operator=(const SliceTensor< T > &t)
Definition tensor.h:2099
SliceTensor(const Tensor< T > &t, const std::array< Slice, TENSOR_MAXDIM > s)
Definition tensor.h:2048
SliceTensor< T > & operator=(const Tensor< T > &t)
Definition tensor.h:2110
SliceTensor< T > & operator=(const T &t)
Definition tensor.h:2121
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
void deallocate()
Definition tensor.h:397
T *MADNESS_RESTRICT _p
Definition tensor.h:321
Tensor< T > & operator=(T x)
Inplace fill tensor with scalar.
Definition tensor.h:553
const Tensor< T > swapdim(long idim, long jdim) const
Returns new view/tensor swaping dimensions i and j.
Definition tensor.h:1614
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
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
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
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition tensor.h:1726
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
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
Tensor(long d0, long d1)
Create and zero new 2-d tensor.
Definition tensor.h:481
const T * ptr() const
Returns a pointer to the internal data.
Definition tensor.h:1829
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
Tensor< T > & unaryop(opT &op)
Inplace apply a unary function to each element of the tensor.
Definition tensor.h:1792
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
Tensor< T > & fill(T x)
Inplace fill with a scalar (legacy name)
Definition tensor.h:562
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< T > & operator+=(const Tensor< Q > &t)
Inplace addition of two tensors.
Definition tensor.h:572
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
Tensor< T > swapdim(long idim, long jdim)
Returns new view/tensor swaping dimensions i and j.
Definition tensor.h:1605
const Tensor< T > flat() const
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition tensor.h:1561
T sum() const
Returns the sum of all elements of the tensor.
Definition tensor.h:1662
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition tensor.h:1384
const T & operator[](long i) const
1-d indexing operation using [] without bounds checking.
Definition tensor.h:791
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
Tensor< T > reshape(long dim0, long dim1)
Returns new view/tensor rehaping to conforming 2-d tensor with given dimensions.
Definition tensor.h:1439
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
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
T & operator()(long i)
1-d indexing operation without bounds checking.
Definition tensor.h:802
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
T & operator()(long i, long j, long k, long l, long m)
5-d indexing operation without bounds checking.
Definition tensor.h:920
T * ptr()
Returns a pointer to the internal data.
Definition tensor.h:1824
bool conforms(const Tensor< Q > &t) const
Test if *this and t conform.
Definition tensor.h:1657
const Tensor< T > mapdim(const std::vector< long > &map) const
Returns new view/tensor permuting the dimensions.
Definition tensor.h:1633
T & operator()(long i, long j)
2-d indexing operation without bounds checking.
Definition tensor.h:825
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
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
const Tensor< T > operator()(const std::vector< Slice > &s) const
General slicing operation (const)
Definition tensor.h:1070
Tensor< T > & operator-=(const Tensor< Q > &t)
Inplace subtraction of two tensors.
Definition tensor.h:582
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
const T & operator()(long i, long j, long k) const
3-d indexing operation without bounds checking.
Definition tensor.h:867
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
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
Tensor< T > mapdim(const std::vector< long > &map)
Returns new view/tensor permuting the dimensions.
Definition tensor.h:1624
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
const Tensor< T > operator()(const std::array< Slice, TENSOR_MAXDIM > &s) const
General slicing operation (const)
Definition tensor.h:1088
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
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
const Tensor< T > fusedim(long i) const
Returns new view/tensor fusing contiguous dimensions i and i+1.
Definition tensor.h:1596
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition tensor.h:686
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type operator-=(const Q &x)
Inplace decrement by scalar of supported type.
Definition tensor.h:707
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
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
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(const Tensor< T > &t)
Copy constructor is shallow (same as assignment)
Definition tensor.h:428
T & operator()(long i, long j, long k)
3-d indexing operation without bounds checking.
Definition tensor.h:852
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
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 > operator-() const
Unary negation producing a new tensor.
Definition tensor.h:663
const T & operator()(long i, long j) const
2-d indexing operation without bounds checking.
Definition tensor.h:838
T trace(const Tensor< T > &t) const
Return the trace of two tensors (no complex conjugate invoked)
Definition tensor.h:1776
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
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, 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)
Return a 1d SliceTensor that views the specified range of the 1d Tensor.
Definition tensor.h:1095
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
const T & operator()(const long ind[]) const
Politically incorrect general indexing operation without bounds checking.
Definition tensor.h:1016
Tensor< T > & emul(const Tensor< T > &t)
Inplace multiply by corresponding elements of argument Tensor.
Definition tensor.h:1798
Tensor< T > & operator=(const Tensor< T > &t)
Assignment is shallow (same as copy constructor)
Definition tensor.h:443
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
T & operator()(long i, long j, long k, long l)
4-d indexing operation without bounds checking.
Definition tensor.h:883
SliceTensor< T > operator()(const std::array< Slice, TENSOR_MAXDIM > &s)
General slicing operation.
Definition tensor.h:1080
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
const T & operator()(long i, long j, long k, long l) const
4-d indexing operation without bounds checking.
Definition tensor.h:901
T product() const
Return the product of all elements of the tensor.
Definition tensor.h:1676
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
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< TENSOR_RESULT_TYPE(T, Q) > operator-(const Tensor< Q > &t) const
Subtraction of two tensors to produce a new tensor.
Definition tensor.h:604
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
T & operator()(const std::vector< long > ind)
General indexing operation with bounds checking.
Definition tensor.h:1032
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
virtual ~Tensor()
Definition tensor.h:1881
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
void allocate(long nd, const long d[], bool dozero)
Definition tensor.h:324
const BaseTensor * base() const
Returns a pointer to the base class.
Definition tensor.h:1839
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
T sumsq() const
Returns the sum of the squares of the elements.
Definition tensor.h:1669
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type operator*=(const Q &x)
Inplace multiplication by scalar of supported type.
Definition tensor.h:675
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
T & operator[](long i)
1-d indexing operation using [] without bounds checking.
Definition tensor.h:780
Tensor< T > & screen(double x)
Inplace set elements of *this less than x in absolute magnitude to zero.
Definition tensor.h:758
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
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)
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 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
const Tensor< T > reshape(long dim0) const
Returns new view/tensor rehapings to conforming 1-d tensor with given dimension.
Definition tensor.h:1429
bool has_data() const
Definition tensor.h:1886
Tensor< T > reshape(const std::vector< long > &d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition tensor.h:1405
BaseTensor * base()
Returns a pointer to the base class.
Definition tensor.h:1834
const TensorIterator< T > & end() const
End point for forward iteration.
Definition tensor.h:1876
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< 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
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
Tensor< T > reshape(long dim0)
Returns new view/tensor rehapings to conforming 1-d tensor with given dimension.
Definition tensor.h:1421
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
Tensor< T > fusedim(long i)
Returns new view/tensor fusing contiguous dimensions i and i+1.
Definition tensor.h:1587
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
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
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
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
const T & operator()(const std::vector< long > ind) const
General indexing operation with bounds checking.
Definition tensor.h:1046
T & operator()(const long ind[])
Politically incorrect general indexing operation without bounds checking.
Definition tensor.h:1000
Tensor< T > flat()
Returns new view/tensor rehshaping to flat (1-d) tensor.
Definition tensor.h:1555
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, long d4)
Create and zero new 5-d tensor.
Definition tensor.h:514
Tensor()
Default constructor does not allocate any data and sets ndim=-1, size=0, _p=0, and id.
Definition tensor.h:415
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
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 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
static bool bounds_checking()
Return true if bounds checking was enabled at compile time.
Definition tensor.h:768
const T & operator()(long i) const
1-d indexing operation without bounds checking.
Definition tensor.h:813
void clear()
Frees all memory and resests to state of default constructor.
Definition tensor.h:1884
TENSOR_SHARED_PTR< T > _shptr
Definition tensor.h:322
Tensor< T > & fillrandom()
Inplace fill with random values ( [0,1] for floats, [0,MAXSIZE] for integers)
Definition tensor.h:724
Tensor(long nd, const long d[], bool dozero=true)
Politically incorrect general constructor.
Definition tensor.h:545
Tensor< T > & gaxpy(T alpha, const Tensor< T > &t, T beta)
Inplace generalized saxpy ... this = this*alpha + other*beta.
Definition tensor.h:1804
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
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
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 > & fillindex()
Inplace fill with the index of each element.
Definition tensor.h:748
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 > 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
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type operator+=(const Q &x)
Inplace increment by scalar of supported type.
Definition tensor.h:696
Tensor< T > & conj()
Inplace complex conjugate.
Definition tensor.h:716
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
SliceTensor< T > operator()(const std::vector< Slice > &s)
General slicing operation.
Definition tensor.h:1060
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
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 > 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< 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
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.
T mynorm(T t)
Definition tensor.h:260
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
std::ostream & operator<<(std::ostream &os, const particle< PDIM > &p)
Definition lowrankfunction.h:397
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
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
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
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
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 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
response_space transpose(response_space &f)
Definition basic_operators.cc:10
std::vector< CCPairFunction< T, NDIM > > operator*(const double fac, const std::vector< CCPairFunction< T, NDIM > > &arg)
Definition ccpairfunction.h:1084
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::enable_if< std::is_base_of< ProjectorBase, projT >::value, OuterProjector< projT, projQ > >::type outer(const projT &p0, const projQ &p1)
Definition projector.h:457
std::vector< CCPairFunction< T, NDIM > > operator-(const std::vector< CCPairFunction< T, NDIM > > c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition ccpairfunction.h:1055
static double pop(std::vector< double > &v)
Definition SCF.cc:113
static void aligned_zero(long n, T *a)
Definition aligned.h:55
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
double imag(double x)
Definition complexfun.h:56
GenTensor< TENSOR_RESULT_TYPE(R, Q)> transform_dir(const GenTensor< R > &t, const Tensor< Q > &c, const int axis)
Definition lowranktensor.h:1099
std::string type(const PairType &n)
Definition PNOParameters.h:18
static const long default_jdim
Definition tensoriter.h:57
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 real(double x)
Definition complexfun.h:52
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
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 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 d
Definition nonlinschro.cc:121
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 double m
Definition relops.cc:9
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.
AtomicInt sum
Definition test_atomicint.cc:46
static const double alpha
Definition testcosine.cc:10
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
#define TENSOR_RESULT_TYPE(L, R)
This macro simplifies access to TensorResultType.
Definition type_data.h:205