Download this article:

The template below offers a simple and efficient solution for handling matrices and tensors in C++. The idea is to store the matrix (or tensor) in a standard vector by translating the multidimensional index to a one dimensional index.

E.g. we can store this matrix

	c0	c1	c2	c3
r0	(0,0)	(0,1)	(0,2)	(0,3)
r1	(1,0)	(1,1)	(1,2)	(1,3)
r2	(2,0)	(2,1)	(2,2)	(2,3)

in a vector like this

0	1	2	3	4	...	11
(0,0)	(0,1)	(0,2)	(0,3)	(1,0)	...	(2,3)

The only thing we need to do is to convert two dimensional indices (r,c) into a one dimensional index. Using template we can do this very efficiently compile time, minimizing the runtime overhead.

Contents 1 Template 2 Example Usage 3 Efficiency and speed 4 Alternative Template: 1-based indices 5 Add a comment

[edit] Template

This template version support up to 4 dimensional tensors.

template <int d1,int d2=1,int d3=1,int d4=1>
class TensorIndex {
  public:
	enum {SIZE = d1*d2*d3*d4 };
	enum {LEN1 = d1 };
	enum {LEN2 = d2 };
	enum {LEN3 = d3 };
	enum {LEN4 = d4 };

    static int indexOf(const int i) {
      return i;
    }
    static int indexOf(const int i,const int j) {
      return j*d1 + i;
    }
    static int indexOf(const int i,const int j, const int k) {
      return (k*d2 + j)*d1 + i;
    }
    static int indexOf(const int i,const int j, const int k,const int l) {
      return ((l*d3 + k)*d2 + j)*d1 + i;
    }
};

[edit] Example Usage

The examples below demonstrate how

• to work with a matrix, using a std::vector as storage
• work with a four dimensional tensor using a simple array as storage

#include <vector>

int main(int argc, char* argv[])
{
    // First example
    //=================================

    // A matrix of 3 rows and 4 columns.
    typedef TensorIndex<3,4> M;

    // Store the matrix in a std::vector
    std::vector<double> M_Storage(M::SIZE);

    // Fill it with ones
    for (int row=0; row<M::LEN1; ++row)
        for (int col=0; col<M::LEN2; ++col)
            M_Storage[M::indexOf(row,col)] = 1.0;


    // Second example
    //=================================
    // A 4 x 5 x 6 x 7 cube
    typedef TensorIndex<4,5,6,7> T;

    // Use an array for storage
    double T_Storage[T::SIZE];

    // Fill it with ones
    for (int i=0; i<T::LEN1; ++i)
        for (int j=0; j<T::LEN2; ++j)
            for (int k=0; k<T::LEN3; ++k)
                for (int l=0; l<T::LEN4; ++l)
                    T_Storage[T::indexOf(i,j,k,l)] = 1.0;

	return 0;
}

[edit] Efficiency and speed

The translation of the indices is done compile time. The code snippet from the above example

typedef TensorIndex<3,4> M;

std::vector<double> M_Storage(M::SIZE);

for (int row=0; row<M::LEN1; ++row)
    for (int col=0; col<M::LEN2; ++col)
        M_Storage[M::indexOf(row,col)] = 1.0;

gets translated to

std::vector<double> M_Storage(12);

for (int row=0; row<3; ++row)
    for (int col=0; col<4; ++col)
        M_Storage[4*row+col] = 1.0;

[edit] Alternative Template: 1-based indices

many numerical algorithms and mathematical articles use a notation where indices start couting from 1 (instead of 0 as in C++). The following template can be used to make the transition from that notation to the zero based C++ standard easy.

template <int d1,int d2=1,int d3=1,int d4=1>
class TensorIndex_base1 {
  public:
	enum {SIZE = d1*d2*d3*d4 };
	enum {LEN1 = d1 };
	enum {LEN2 = d2 };
	enum {LEN3 = d3 };
	enum {LEN4 = d4 };

    static int indexOf(const int i) {
      return i -1;
    }
    static int indexOf(const int i,const int j) {
      return j*d1 + i -1 -d1;
    }
    static int indexOf(const int i,const int j, const int k) {
      return (k*d2 + j)*d1 + i -1 -d1 -d1*d2;
    }
    static int indexOf(const int i,const int j, const int k,const int l) {
      return ((l*d3 + k)*d2 + j)*d1 + i -1 -d1 -d1*d2 - d1*d2*d3;
    }
};

...


// Example Usage


    typedef TensorIndex_base1<3,4> M;

    double M_Storage[M::SIZE];

    for (int i=1; i<=M::LEN1; ++i)
        for (int j=1; j<=M::LEN2; ++j)
            M_Storage[M::indexOf(i,j)] = 1.0;

[edit] Add a comment

Talk:A Simple and Extremely Fast C++ Template for Matrices and Tensors