Tensors or multidimension arrays

Collaboration diagram for Tensors or multidimension arrays:

Files
file	slice.h
	Declares and implements Slice.
file	tensor.cc
	Completes the implementation of Tensor and instantiates all specializations for fast compiles.
file	tensor.h
	Defines and implements most of Tensor.
file	tensor_macros.h
	Macros for easy and efficient iteration over tensors.
file	tensortrain.h
	Defines and implements the tensor train decomposition as described in I.V. Oseledets, Siam J. Sci. Comput. 33, 2295 (2011).

Namespaces
	madness
	Holds machinery to set up Functions/FuncImpls using various Factories and Interfaces.

Classes
class	madness::BaseTensor
	The base class for tensors defines generic capabilities. More...
class	madness::Slice
	A slice defines a sub-range or patch of a dimension. More...
class	madness::TensorException
	Tensor is intended to throw only TensorExceptions. More...
class	madness::TensorIterator< T, Q, R >

Functions
template<class T >
GenTensor< T >	madness::outer (const GenTensor< T > &left, const GenTensor< T > &right)
	Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...) More...
template<class T >
GenTensor< T >	madness::outer_low_rank (const Tensor< T > &left, const Tensor< T > &right)
	Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...) More...
template<typename T , typename Q >
IsSupported< TensorTypeData< Q > , Tensor< T > >::type	madness::operator+ (Q x, const Tensor< T > &t)
	The class defines tensor op scalar ... here define scalar op tensor. More...
template<typename T , typename Q >
IsSupported< TensorTypeData< Q > , Tensor< T > >::type	madness::operator* (const Q &x, const Tensor< T > &t)
	The class defines tensor op scalar ... here define scalar op tensor. More...
template<typename T , typename Q >
IsSupported< TensorTypeData< Q > , Tensor< T > >::type	madness::operator- (Q x, const Tensor< T > &t)
	The class defines tensor op scalar ... here define scalar op tensor. More...
template<class T >
Tensor< T >	madness::copy (const Tensor< T > &t)
	Returns a new contiguous tensor that is a deep copy of the input. More...
template<class T , class Q >
Tensor< TENSOR_RESULT_TYPE(T, Q)>	madness::transform_dir (const Tensor< T > &t, const Tensor< Q > &c, int axis)
	Transforms one dimension of the tensor t by the matrix c, returns new contiguous tensor. More...
template<class T >
Tensor< T >	madness::transpose (const Tensor< T > &t)
	Returns a new deep copy of the transpose of the input tensor. More...
template<class T >
Tensor< T >	madness::conj_transpose (const Tensor< T > &t)
	Returns a new deep copy of the complex conjugate transpose of the input tensor. More...
template<class T >
std::ostream &	madness::operator<< (std::ostream &s, const Tensor< T > &t)
	Print (for human consumption) a tensor to the stream. More...
template<class T >
Tensor< T >	madness::outer (const Tensor< T > &left, const Tensor< T > &right)
	Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...) More...
template<class T , class Q >
Tensor< TENSOR_RESULT_TYPE(T, Q)>	madness::inner (const Tensor< T > &left, const Tensor< Q > &right, long k0=-1, long k1=0)
	Inner product ... result(i,j,...,p,q,...) = sum(z) left(i,j,...,z)*right(z,p,q,...) More...
template<class T , class Q >
void	madness::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. More...
template<class T , class Q >
Tensor< TENSOR_RESULT_TYPE(T, Q)>	madness::transform (const Tensor< T > &t, const Tensor< Q > &c)
	Transform all dimensions of the tensor t by the matrix c. More...
template<class T , class Q >
Tensor< TENSOR_RESULT_TYPE(T, Q)>	madness::general_transform (const Tensor< T > &t, const Tensor< Q > c[])
	Transform all dimensions of the tensor t by distinct matrices c. More...
template<class T , class Q >
Tensor< TENSOR_RESULT_TYPE(T, Q) > &	madness::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() More...
template<class T >
Tensor< typename Tensor< T > ::scalar_type >	madness::abs (const Tensor< T > &t)
	Return a new tensor holding the absolute value of each element of t. More...
template<class T >
Tensor< typename Tensor< T > ::scalar_type >	madness::arg (const Tensor< T > &t)
	Return a new tensor holding the argument of each element of t (complex types only) More...
template<class T >
Tensor< typename Tensor< T > ::scalar_type >	madness::real (const Tensor< T > &t)
	Return a new tensor holding the real part of each element of t (complex types only) More...
template<class T >
Tensor< typename Tensor< T > ::scalar_type >	madness::imag (const Tensor< T > &t)
	Return a new tensor holding the imaginary part of each element of t (complex types only) More...
template<class T >
Tensor< T >	madness::conj (const Tensor< T > &t)
	Returns a new deep copy of the complex conjugate of the input tensor (complex types only) More...
	madness::TensorIterator< T, Q, R >::TensorIterator (const Tensor< T > *t0, const Tensor< Q > *t1=0, const Tensor< R > *t2=0, long iterlevel=0, bool optimize=true, bool fusedim=true, long jdim=default_jdim)
	Constructor for general iterator to compose operations over up to three tensors. More...

Detailed Description

Introduction

A tensor is a multi-dimensional array and does not incorporate any concepts of covariance and contravariance.

When a new tensor is created, the underlying data is also allocated. E.g.,

Tensor<double> a(3,4,5)

creates a new 3-dimensional tensor and allocates a contiguous block of 60 doubles which are initialized to zero. The dimensions (numbered from the left starting at 0) are in C or row-major order. Thus, for the tensor a , the stride between successive elements of the right-most dimension is 1. For the middle dimension it is 5. For the left-most dimension it is 20. Thus, the loops

for (i=0; i<3; ++i)
    for (j=0; j<4; ++j)
        for (k=0; k<5; ++k)
            a(i,j,k) = ...

will go sequentially (and thus efficiently) through memory. If the dimensions have been reordered (e.g., with swapdim() or map() ), or if the tensor is actually a slice of another tensor, then the layout in memory may be more complex and may not reflect a contiguous block of memory.

Multiple tensors may be used to provide multiple identical or distinct views of the same data. E.g., in the following

Tensor<double> a(2,3);  // A new tensor initialized to zero
Tensor<double> b = a;

a and b provide identical views of the same data, thus

b(1,2) = 99;
cout << a(1,2) << endl;  // Outputs 99
cout << b(1,2) << endl;  // Outputs 99

Shallow copy and assignment

It is important to appreciate that the views and the data are quite independent. In particular, the default copy constructor and assignment operations only copy the tensor (the view) and not the data — i.e., the copy constructor and assigment operations only take shallow copies. This is for both consistency and efficiency. Thus, assigning one tensor to another generates another view of the same data, replacing any previous view and not moving or copying any of the data. E.g.,

Tensor<double> a(2,3);   // A new tensor initialized to zero
Tensor<double> c(3,3,3); // Another new tensor
Tensor<double> b = a;    // b is a view of the same data as a
a = c;                   // a is now a view of c's data
b = c                    // b is now also a view of c's data and the
// data allocated originally for a is freed

The above example also illustrates how reference counting is used to keep track of the underlying data. Once there are no views of the data, it is automatically freed.

There are only two ways to actually copy the underlying data. A new, complete, and contigous copy of a tensor and its data may be generated with the copy() function. Or, to copy data from one tensor into the data viewed by another tensor, you must use a Slice.

Indexing

One dimensional tensors (i.e., vectors) may be indexed using either square brackets (e.g., v[i] ) or round brackets (e.g., v(i) ). All higher-dimensional tensors must use only round brackets (e.g., t(i,j,k) ). This is due to C++'s restriction that the indexing operator ([] ) can only have one argument. The indexing operation should generate efficient code.

For the sake of efficiency, no bounds checking is performed by default by most single element indexing operations. Checking can be enabled at compile time by defining -DTENSOR_BOUNDS_CHECKING for application files including tensor.h. The MADNESS configure script has the option –enable-tensor-bound-checking to define the macro in madness_config.h . The general indexing operation that takes a std::vector<long> index and all slicing operations always perform bounds checking. To make indexing with checking a bit easier, a factory function has been provided for vectors ... but note you need to explicitly use longs as the index.

Tensor<long> a(7,7,7);
a(3,4,5) += 1;                    // OK ... adds 1 to element (3,4,5)
a(3,4,9) += 1;                    // BAD ... undetected out-of-bounds access
a(vector_factory(3L,4L,9L)) += 1; // OK ... out-bounds access will
// be detected at runtime.

Slicing

Slices generate sub-tensors — i.e., views of patches of the data. E.g., to refer to all but the first and last elements in each dimension of a matrix use

a(Slice(1,-2),Slice(1,-2))

Or to view odd elements in each dimension

a(Slice(0,-1,2),Slice(0,-1,2))

A slice or patch of a tensor behaves exactly like a tensor except for assignment. When a slice is assigned to, the data is copied with the requirement that the source and destinations agree in size and shape (i.e., they conform). Thus, to copy the all of the data from a to b,

Tensor<double> a(3,4), b(3,4);
a = 1;                              // Set all elements of a to 1
b = 2;                              // Set all elements of b to 2
a(Slice(0,-1,1),Slice(0,-1,1)) = b; // Copy all data from b to a
a(_,_) = b(_,_);                    // Copy all data from b to a
a(___) = b(___);                    // Copy all data from b to a
a(Slice(1,2),Slice(1,2) = b;        // Error, do not conform

Special slice values _ ,_reverse, and ___ have been defined to refer to all elements in a dimension, all elements in a dimension but reversed, and all elements in all dimensions, respectively.

Iteration and algorithms

See tensor_macros.h for documentation on the easiest mechanisms for iteration over elements of tensors and tips for optimization. See TensorIterator for the most general form of iteration.

Function Documentation

template<class T >

Tensor< typename Tensor<T>::scalar_type > madness::abs

(

const Tensor< T > &

t

)

Return a new tensor holding the absolute value of each element of t.

References BINARY_OPTIMIZED_ITERATOR.

template<class T >

Tensor< typename Tensor<T>::scalar_type > madness::arg

(

const Tensor< T > &

t

)

Return a new tensor holding the argument of each element of t (complex types only)

References BINARY_OPTIMIZED_ITERATOR.

Referenced by madness::alloc_am_arg(), madness::World::args(), madness::FunctionImpl< Q, NDIM >::do_binary_op(), madness::FunctionImpl< Q, NDIM >::do_mul(), madness::Solver< T, NDIM >::do_rhs(), madness::Solver< T, NDIM >::do_rhs_simple(), madness::free_am_arg(), ScatteringWF::gamma(), main(), madness::new_am_arg(), testing::internal::ParseGoogleTestFlagsOnlyImpl(), psi_exact(), madness::startup(), and madness::FunctionImpl< Q, NDIM >::traverse_tree().

template<class T >

Tensor<T> madness::conj

(

const Tensor< T > &

t

)

Returns a new deep copy of the complex conjugate of the input tensor (complex types only)

References BINARY_OPTIMIZED_ITERATOR.

template<class T >

Tensor<T> madness::conj_transpose

(

const Tensor< T > &

t

)

Returns a new deep copy of the complex conjugate transpose of the input tensor.

References madness::conj(), and TENSOR_ASSERT.

Referenced by madness::Solver< T, NDIM >::build_fock_matrix(), madness::Solver< T, NDIM >::do_rhs(), madness::Solver< T, NDIM >::do_rhs_simple(), madness::Solver< T, NDIM >::initial_guess(), madness::my_conj_transpose(), and madness::Solver< T, NDIM >::print_fock_matrix_eigs().

template<class T >

Tensor<T> madness::copy

(

const Tensor< T > &

t

)

Returns a new contiguous tensor that is a deep copy of the input.

Returns

Returns a new contiguous tensor that is a deep copy of the input

References BINARY_OPTIMIZED_ITERATOR.

template<class T , class Q >

Tensor< TENSOR_RESULT_TYPE(T,Q) >& madness::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()

Both dimensions of c must be the same and match all dimensions of the input tensor t. All tensors must be contiguous.

Performs the same operation as transform but it requires that the caller pass in workspace and a preallocated result, hoping that that both can be reused. If the result and workspace are reused between calls, then no tensor constructors need be called and cache locality should be improved. By passing in the workspace, this routine is kept thread safe.

The input, result and workspace tensors must be distinct.

All input tensors must be contiguous and fastest execution will result if all dimensions are even and data is aligned on 16-byte boundaries. The workspace and the result must be of the same size as the input t . The result tensor need not be initialized before calling fast_transform.

result(i,j,k,...) <-- sum(i',j', k',...) t(i',j',k',...)  c(i',i) c(j',j) c(k',k) ...

The input dimensions of t must all be the same .

References IS_ODD, IS_UNALIGNED, mTxm(), madness::mTxmq(), madness::mTxmq_padding(), std::swap(), and madness::TENSOR_RESULT_TYPE().

Referenced by madness::FunctionImpl< Q, NDIM >::err_box(), madness::FunctionImpl< T, NDIM >::filter(), main(), madness::FunctionImpl< T, NDIM >::project(), Test7(), madness::transform(), and madness::FunctionImpl< T, NDIM >::unfilter().

template<class T , class Q >

Tensor<TENSOR_RESULT_TYPE(T,Q)> madness::general_transform	(	const Tensor< T > &	t,
		const Tensor< Q >	c[]
	)

Transform all dimensions of the tensor t by distinct matrices c.

Similar to transform but each dimension is transformed with a distinct matrix.

result(i,j,k...) <-- sum(i',j', k',...) t(i',j',k',...) c[0](i',i) c[1](j',j) c[2](k',k) ...

The first dimension of the matrices c must match the corresponding dimension of t.

References madness::inner(), and madness::TENSOR_RESULT_TYPE().

template<class T >

Tensor< typename Tensor<T>::scalar_type > madness::imag

(

const Tensor< T > &

t

)

Return a new tensor holding the imaginary part of each element of t (complex types only)

References BINARY_OPTIMIZED_ITERATOR.

template<class T , class Q >

Tensor<TENSOR_RESULT_TYPE(T,Q)> madness::inner	(	const Tensor< T > &	left,
		const Tensor< Q > &	right,
		long	k0 = -1,
		long	k1 = 0
	)

Inner product ... result(i,j,...,p,q,...) = sum(z) left(i,j,...,z)*right(z,p,q,...)

By default it contracts the last dimension of the left tensor and the first dimension of the right tensor. These defaults can be changed by specifying k0 and k1 , the index to contract in the left and right side tensors, respectively. The defaults correspond to (k0=-1 and k1=0 ).

References madness::inner_result(), TENSOR_ASSERT, and TENSOR_MAXDIM.

template<class T , class Q >

void madness::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.

This routine may be used to optimize away the tensor constructor of the result tensor in inner loops when the result tensor may be reused or accumulated into. If the user calls this routine directly very little checking is done since it is intended as an optimization for small tensors. As far as the result goes, the caller is completely responsible for providing a contiguous tensor that has the correct dimensions and is appropriately initialized. The inner product is accumulated into result.

References madness::TensorIterator< T, Q, R >::_p0, madness::TensorIterator< T, Q, R >::_s0, mTxm(), mTxmT(), mxm(), mxmT(), madness::TensorIterator< T, Q, R >::reset(), restrict, mpfr::sum(), and madness::TENSOR_RESULT_TYPE().

Referenced by madness::inner(), and madness::ortho5().

template<typename T , typename Q >

IsSupported< TensorTypeData<Q>, Tensor<T> >::type madness::operator*	(	const Q &	x,
		const Tensor< T > &	t
	)

The class defines tensor op scalar ... here define scalar op tensor.

template<typename T , typename Q >

IsSupported< TensorTypeData<Q>, Tensor<T> >::type madness::operator+	(	Q	x,
		const Tensor< T > &	t
	)

The class defines tensor op scalar ... here define scalar op tensor.

template<typename T , typename Q >

IsSupported< TensorTypeData<Q>, Tensor<T> >::type madness::operator-	(	Q	x,
		const Tensor< T > &	t
	)

The class defines tensor op scalar ... here define scalar op tensor.

template<class T >

std::ostream & madness::operator<<	(	std::ostream &	s,
		const Tensor< T > &	t
	)

Print (for human consumption) a tensor to the stream.

template<class T >

GenTensor<T> madness::outer	(	const GenTensor< T > &	left,
		const GenTensor< T > &	right
	)

Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...)

Referenced by madness::FunctionImpl< T, NDIM >::Vphi_op_NS< opT, LDIM >::make_childrens_sum_coeffs(), madness::TwoElectronInterface< double, 6 >::make_coeff(), madness::FunctionImpl< T, NDIM >::Vphi_op_NS< opT, LDIM >::make_sum_coeffs(), madness::FunctionImpl< T, NDIM >::hartree_op< LDIM, leaf_opT >::operator()(), madness::ortho5(), madness::outer_low_rank(), and Test7().

template<class T >

Tensor<T> madness::outer	(	const Tensor< T > &	left,
		const Tensor< T > &	right
	)

Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...)

References madness::TensorIterator< T, Q, R >::_p0, madness::TensorIterator< T, Q, R >::_s0, madness::TensorIterator< T, Q, R >::dimj, madness::TensorIterator< T, Q, R >::reset(), TENSOR_ASSERT, and TENSOR_MAXDIM.

template<class T >

GenTensor<T> madness::outer_low_rank	(	const Tensor< T > &	left,
		const Tensor< T > &	right
	)

Outer product ... result(i,j,...,p,q,...) = left(i,k,...)*right(p,q,...)

References madness::outer().

Referenced by madness::FunctionImpl< T, NDIM >::recursive_apply_op< opT, LDIM >::operator()().

template<class T >

Tensor< typename Tensor<T>::scalar_type > madness::real

(

const Tensor< T > &

t

)

Return a new tensor holding the real part of each element of t (complex types only)

References BINARY_OPTIMIZED_ITERATOR.

template<class T , class Q , class R >

madness::TensorIterator< T, Q, R >::TensorIterator	(	const Tensor< T > *	t0,
		const Tensor< Q > *	t1 = 0,
		const Tensor< R > *	t2 = 0,
		long	iterlevel = 0,
		bool	optimize = true,
		bool	fusedim = true,
		long	jdim = default_jdim
	)

Constructor for general iterator to compose operations over up to three tensors.

Macros have been defined in tensor_macros.h to take the pain out of using iterators. The options have the following effects

• optimize reorders dimensions for optimal strides (only applies if iterlevel=1). If jdim==default_jdim, all dimensions are reordered for optimal stride. If jdim is not the default value, then dimension jdim is excluded from the set of dimensions being optimized.
• fusedim concatenates contiguous dimensions into the inner loop (only applies if iterlevel=1 and jdim=default_jdim).
• iterlevel can have two values 0=elementwise, 1=vectorwise. Elementwise implies that the iterator returns successive elements, and explicitly in the expected order for the tensor (i.e., with the last index varying fastest). Vectorwise implies that the user is responsible for iterating over dimension jdim.
• jdim — if iterlevel=1, then jdim determines which dimension is left for iteration with an explicit for loop. Negative values implies jdim+=ndim following the convention of Slice (and Python). If (optimize) the default is the fastest varying dimension. If (!optimize) the default is the last dimension (jdim=-1). If (fusedim) contiguous dimensions are fused into this inner loop. Specifying a non-default value for jdim disables fusedim and restricts optimization to reordering only the exterior loops (so that the loop the user is iterating over corresponds exactly to those in dimension jdim).

During iteration:

• ind[] will contain the current iteration indices .. BUT if optimize=true , they will not necessarily be in the order of those of the tensor. if fusedim is true, then there may be fewer dimensions than the input tensor.
• _p0, _p1, _p2 will point to the current elements of t0,t1,t2 noting that if iterlevel>0, the user must provide additional for loops to iterate over the additional dimensions
• stride0[], stride1[], stride2[] will contain the strides for each (possibly reordered) dimension. If iterlevel=1, _s0,_s1,s2 contain the stride info for the dimension that the user is responsible for iterating over.
• dimj -> the size of the j'th dimension

References mpfr::dim(), R, std::swap(), std::tr1::T(), TENSOR_ASSERT, and TENSOR_MAXDIM.

template<class T , class Q >

Tensor<TENSOR_RESULT_TYPE(T,Q)> madness::transform	(	const Tensor< T > &	t,
		const Tensor< Q > &	c
	)

Transform all dimensions of the tensor t by the matrix c.

Often used to transform all dimensions from one basis to another

result(i,j,k...) <-- sum(i',j', k',...) t(i',j',k',...) c(i',i) c(j',j) c(k',k) ...

The input dimensions of t must all be the same and agree with the first dimension of c . The dimensions of c may differ in size. If the dimensions of c are the same, and the operation is being performed repeatedly, then you might consider calling fast_transform instead which enables additional optimizations and can eliminate all constructor overhead and improve cache locality.

References c, madness::fast_transform(), madness::inner(), TENSOR_ASSERT, and madness::TENSOR_RESULT_TYPE().

template<class T , class Q >

Tensor<TENSOR_RESULT_TYPE(T,Q)> madness::transform_dir	(	const Tensor< T > &	t,
		const Tensor< Q > &	c,
		int	axis
	)

Transforms one dimension of the tensor t by the matrix c, returns new contiguous tensor.

transform_dir(t,c,1) = r(i,j,k,...) = sum(j') t(i,j',k,...) * c(j',j)

Parameters

[in]	t	Tensor to transform (size of dimension to be transformed must match size of first dimension of c )
[in]	c	Matrix used for the transformation
[in]	axis	Dimension (or axis) to be transformed

Returns

Returns a new, contiguous tensor

References madness::copy(), and madness::inner().

Referenced by madness::FunctionImpl< Q, NDIM >::apply_1d_realspace_push_op().

template<class T >

Tensor<T> madness::transpose

(

const Tensor< T > &

t

)

Returns a new deep copy of the transpose of the input tensor.

References madness::copy(), and TENSOR_ASSERT.

Referenced by madness::FunctionCommonData< T, NDIM >::_init_quadrature(), madness::SCF::analyze_vectors(), madness::Convolution1D< Q >::Convolution1D(), SCF::diag_fock_matrix(), distributed_localize_PM(), madness::gelss(), madness::geqp3(), madness::gesv(), madness::SCF::get_fock_transformation(), madness::SCF::initial_guess(), SCF::localize_boys(), main(), SCF::make_fock_matrix(), madness::SCF::make_fock_matrix(), madness::my_conj_transpose(), madness::orgqr(), madness::ortho5(), madness::qr(), madness::Convolution1D< Q >::rnlij(), madness::syev(), madness::sygv(), madness::test_syev(), and madness::transform().