Thomas Russell: Programmer & Designer
Have a browse of my thoughts!

Writing a Tensor class in Modern C++

As promised in my article about the mathematics of solving the multidimensional Poisson equation, I am beginning writing about the process of writing a parallelized solution to the Poisson Equation by successive over-relaxation. This article will be about writing a tensor class in modern C++, as we will need a multidimensional grid for our Poisson solver.

Mathematically speaking, a Tensor is a multidimensional object over some appropriate field. In this case we will be considering simple fields, real numbers (to floating-point precision), integers (modulo some appropriate power of 2) and complex numbers. Tensors are mathematically-rich objects, with their own notation and algebras embedded in Ricci calculus, but here we will only be writing a basic Tensor class equipped with a few useful operators, however, it is an easily extendable object, so can easily be modified for most general purposes.

Defining our Tensor class

We begin by defining our barebones Tensor class as follows (which we will call CTensor):

// Our Tensor class:
template <typename T, std::size_t... Sizes>
class CTensor {
        
	// Check that T is trivivally copyable:
	static_assert(std::is_trivially_copyable<T>::value, "We cannot create a Tensor of a non-trivially copyable type");

public:
	CTensor() {
		m_ptValues = new T[total_extent<Sizes...>::value]();
	}

	CTensor( const CTensor& tenRHS )
	{
		m_ptValues = new T[total_extent<Sizes...>::value]();
		std::memcpy( m_ptValues, tenRHS.m_ptValues, total_extent<Sizes...>::value * sizeof(T) );
	}

	template <typename U>
	CTensor( const CTensor<U, Sizes...>& tenRHS )
	{
		static_assert(std::is_convertible<U, T>::value, "Cannot construct a matrix from a matrix of non-convertible type");
		m_ptValues = new T[total_extent<Sizes...>::value]();
		
		for( std::size_t s = 0U; s < total_extent<Sizes...>::value; ++s )
		{
			m_ptValues[s] = static_cast<T>(tenRHS.m_ptValues[s]);
		}
	}

	CTensor( CTensor&& tenRHS )
	{
		m_ptValues = tenRHS.m_ptValues;
		tenRHS.m_ptValues = nullptr;
	}

	~CTensor()
	{
		if( m_ptValues )
			delete[] m_ptValues;
	}

public:
	CTensor&	operator=(CTensor&& tenRHS)
	{
		std::swap( m_ptValues, tenRHS.m_ptValues );

		return *this;
	}

	CTensor&	operator=(const CTensor& tenRHS)
	{
		std::memcpy( m_ptValues, tenRHS.m_ptValues, total_extent<Sizes...>::value * sizeof(T) );

		return *this;
	}

	template <typename U>
	CTensor&	operator=(const CTensor<U, Sizes...>& tenRHS)
	{
		static_assert(std::is_convertible<U, T>::value, "Cannot copy-assign a matrix from a matrix of non-convertible type");

		for( std::size_t s = 0U; s < total_extent<Sizes...>::value; ++s )
		{
			m_ptValues[s] = static_cast<T>(tenRHS.m_ptValues[s]);
		}

		return *this;
	}

public:
	// Our operators will go here!!

private:
	T*			m_ptValues;
	std::size_t m_sTotalSize = total_extent<Sizes...>::value;
};

We note two important things here: Firstly, the Tensor is of a fixed size; this is because in most numeric methods the dimensions of the problem are pre-specified and can be determined at compile time; Secondly, the memory required is allocated on the heap, this is done to prevent a Stack Overflow occurring if we specify a particularly large Tensor (if we had a $256 \times 256 \times 256$ tensor of double then this would require over 134MB, too much for the average stack).

A few useful helper methods

The observant amongst you will have noticed that we make extensive use of a variadic template class total_extent, this is a helper class which allows us to calculate the total extent (i.e. the product of the dimensions) of the Tensor using compile-time arithmetic. We define the helper struct as follows:

// Helper struct to determine total extent:
template <std::size_t Size1, std::size_t... Sizes>
struct total_extent {
	static const std::size_t value = Size1 * total_extent<Sizes...>::value;
};

template <std::size_t Size>
struct total_extent<Size> {
	static const std::size_t value = Size;
};

We can also define another helpful method (which will prove it’s usefulness when we write our indexing operators later), which checks to make sure all of the types of a variadic parameter pack are the same:

// Helper struct to check type-similarity:
template <typename T, typename U, typename... Args>
struct check_same {
	static const bool value = std::is_same<T, U>::value && check_same<U, Args...>::value;
};

template <typename T, typename U>
struct check_same<T, U> {
	static const bool value = std::is_same<T, U>::value;
};

Our last helper method is one which allows us to transform an n-tuple of coordinates (with the most-significant coordinate first) into a one-dimensional memory reference. We do this by considering the following mathematical function $f:\mathbb{N}^{n} \to \mathbb{N}$:

$$f:\begin{pmatrix}x_{n} & \dots & x_{1}\end{pmatrix} \mapsto \sum_{k=1}^{n}(n-k+1)x_{k},$$

Unfortunately, the parameter pack subscripting syntax Coords[k] is not yet valid C++ (although there do exist proposals) so instead we have to develop a recursive function:

template <std::size_t... Coords>
std::size_t coord_to_reference( std::size_t sCurrIndex, std::size_t sPos, Coords... coords )
{
	return sCurrIndex * sPos + coord_to_reference( sCurrIndex - 1U, coords... );
}

std::size_t coord_to_reference( std::size_t sCurrIndex )
{
	return 0U;
}

Writing the indexing operator

We are now ready to write the const and non-const indexing operators for our Tensor class, using our previously specified helper functions. We can write these as:

template <typename... Coords>
T&			operator()( Coords... coords )
{
	static_assert(check_same<std::size_t, Coords...>::value || check_same<int, Coords...>::value, "The coordinate type must be std::size_t or int");
	static_assert(sizeof...(Coords) == sizeof...(Sizes), "There must be the same number of coordinates as the number of dimensions of the Tensor");
	return m_ptValues[coord_to_reference( sizeof...(Sizes), coords... )];
}

template <typename... Coords>
const T&	operator()( Coords... coords ) const
{
	static_assert(check_same<std::size_t, Coords...>::value || check_same<int, Coords...>::value, "The coordinate type must be std::size_t or int");
	static_assert(sizeof...(Coords) == sizeof...(Sizes), "There must be the same number of coordinates as the number of dimensions of the Tensor");
	return m_ptValues[coord_to_reference( sizeof...(Sizes), coords... )];
}

We’ve now completed everything required for the Tensor class to be used in the multidimensional Poisson equation solver.

Extending the Tensor class

Whilst we don’t need anything more in our Tensor class for our Poisson solver, it may be useful for other applications to have a few more operators, for instance, simple Tensor addition/subtraction, equality checking and scalar multiplication:

template <typename U>
auto		operator+( const CTensor<U, Sizes...>& tenRHS ) const -> CTensor<decltype(std::declval<T>() + std::declval<U>()), Sizes...>
{
	CTensor<decltype(std::declval<T>() + std::declval<U>()), Sizes...> tenResult;

	for( std::size_t s = 0U; s < total_extent<Sizes...>::value; ++s )
	{
		tenResult.m_ptValues[s] = m_ptValues[s] + tenRHS.m_ptValues[s];
	}

	return tenResult;
}

template <typename U>
auto		operator-( const CTensor<U, Sizes...>& tenRHS ) const -> CTensor<decltype(std::declval<T>() - std::declval<U>()), Sizes...>
{
	CTensor<decltype(std::declval<T>() + std::declval<U>()), Sizes...> tenResult;

	for( std::size_t s = 0U; s < total_extent<Sizes...>::value; ++s )
	{
		tenResult.m_ptValues[s] = m_ptValues[s] - tenRHS.m_ptValues[s];
	}

	return tenResult;
}

template <typename U>
auto		operator*(U uVal) const -> CTensor<decltype(std::declval<T>() * std::declval<U>()), Sizes...>
{
	CTensor<decltype(std::declval<T>() * std::declval<U>()), Sizes...> tenResult;

	for( std::size_t s = 0U; s < total_extent<Sizes...>::value; ++s )
	{
		tenResult.m_ptValues[s] = m_ptValues[s] * uVal;
	}

	return tenResult;
}

template <typename U>
bool		operator==(const CTensor<U, Sizes...>& tenRHS) const
{
	static_assert(std::is_convertible<U, T>::value, "Cannot compare Tensors of non-convertible values");

	for( std::size_t s = 0U; s < total_extent<Sizes...>::value; ++s )
	{
		if( m_ptValues[s] != static_cast<T>(tenRHS[s]) )
			return false;
	}

	return true;
}

Other useful operators can be defined as-needed, and I will probably define more as I use them in other posts. We can also define a Tensor product function, which is a binary operator which takes two tensors $T \in \mathcal{F}^{n_{1} \times \cdots \times n_{\ell}}$ and $U \in \mathcal{F}^{m_{1} \times \cdots \times m_{p}}$ and returns a third tensor $T \otimes U \in \mathcal{F}^{n_{1} \times \cdots \times n_{\ell} \times m_{1} \times \cdots \times m_{p}}$. This tensor can be defined component-wise by:

$$(T \otimes U)_{i_{1},\dots,i_{\ell},j_{1},\dots,j_{p}} \triangleq T_{i_{1},\dots,i_{\ell}} \times U_{j_{1},\dots,j_{p}}$$

We can define this programmatically as a friend function of CTensor as follows:

template <typename T1, typename T2, std::size_t... Sizes1, std::size_t... Sizes2>
CTensor<decltype(std::declval<T1>() * std::declval<T2>()), Sizes1..., Sizes2...> TensorProduct( const CTensor<T1, Sizes1...>& tensor1, const CTensor<T2, Sizes2...>& tensor2 )
{
	CTensor<decltype(std::declval<T1>() * std::declval<T2>()), Sizes1..., Sizes2...> tenResult;

	for( std::size_t s1 = 0U; s1 < total_extent<Sizes1...>::value; ++s1 )
	{
		for( std::size_t s2 = 0U; s2 < total_extent<Sizes2...>::value; ++s2 )
		{
			tenResult.m_ptValues[total_extent<Sizes2...>::value * s1 + s2] = tensor1.m_ptValues[s1] * tensor2.m_ptValues[s2];
		}
	}

	return tenResult;
}

This is all we’re going to do now with regards to writing a Tensor class; in future blog posts I will show applications involving the CTensor class and various extensions required to do these.

If you found this blog post interesting, you might like the rest of my numerical methods related posts.

ZIP File containing .cpp file: TensorClass

Leave a Reply

Your email address will not be published. Required fields are marked *

Designed and Produced by Thomas Russell © 2014-2017

Log-in | Register