Allocate data in nested STL vectors contigously

Question

I want to store a 3-dimensional matrix in a way that all the data is contiguous. This is because I need to send the matrix to a master node using MPI, and then concantenate all the smaller matrices into one big matrix. This becomes cumbersome with a nested vector, since you can only send more than one int, double, etc. if they are contiguous in memory, so you have to send each of the innermost vectors separately, using something like

for (int i = 0; i < nx; i++) {
    for (int j = 0; j < ny; j++) {
        int id_of_data = i*ny*nz + j*nz;
        MPI_Send(&local_matrix[i][j][0], nz, MPI_DOUBLE, 0, id_of_data, MPI_COMM_WORLD);
    }
}

and then receive the data on the master node, using id_of_data Right now I'm using std vectors for storing the matrix,

std::vector<std::vector<std::vector<double> > > matrix;

so I'm trying to stick to std::vector for now. I realize it's easy to solve this using arrays/pointers, but I want to try this using std::vectors before resorting to that. The dimensions of the matrix are constant, so I don't have to worry about dynamic allocation.

I can access the elements of the vector using pointers like this

double* a = &matrix[x][y][z];

but I don't know how to do "the opposite", ie. something like

double* vec_ptr = new float(nz);
&matrix[x][y][0] = vec_ptr;

What I would like to do is

double* linear_matrix = new double(nx*ny*nz);

and then somehow make the nested vectors point to this data.

Answer 1

I ended up creating my own classes for this kind of situations. It is very typical for matrix manipulation libraries to take a matrix as a contiguous block of data, but at the same time I would like to be able to fill and read it using the standard m[i][j] accessors. So, here you go. Below is the code to create containers in arbitrary dimensions as a contiguous block of memory. Using boost preprocessor here to be able to do any dimension. If you only need it up to 3d it is quite obvious to de-boostify it :)

memblock.h

#ifndef BOOST_PP_IS_ITERATING

#ifndef MEMBLOCK_INCLUDED
#define MEMBLOCK_INCLUDED
#include <algorithm>
#include <boost/preprocessor.hpp>

#    ifndef MBLOCK_MAX_DIM
#      define MBLOCK_MAX_DIM 4
#    endif

template <typename T> class MemoryBlock1d
{
protected:
    T* m_data;
public:
    MemoryBlock1d() { m_data = 0; }
    MemoryBlock1d(int size) { m_data = 0;  resize(size); }
    MemoryBlock1d(int size, T val) { m_data = 0;  resize(size, val); }
    virtual ~MemoryBlock1d();
    inline T*& data() { return m_data; }
    inline void resize(int size);
    inline void resize(int size, T val);
    inline void clear();
    inline T& operator[](int i) { return m_data[i]; }
};

// Implementation 1d:
template <typename T> inline MemoryBlock1d<T>::~MemoryBlock1d()
{
    if (m_data)
    {
        delete[] m_data;
    }
}

template <typename T> inline void MemoryBlock1d<T>::clear()
{
    if (m_data)
    {
        delete[] m_data;
    }
    m_data = 0;
}

template <typename T> inline void MemoryBlock1d<T>::resize(int size)
{
    if (m_data)
    {
        delete[] m_data;
    }
    m_data = new T[size];
}

template <typename T> void MemoryBlock1d<T>::resize(int size, T val)
{
    resize(size);
    std::fill(m_data, m_data + size, val);
}


// generate multi_dimensional blocks:
#    define BOOST_PP_ITERATION_LIMITS (2, MBLOCK_MAX_DIM)
#    define BOOST_PP_FILENAME_1       "memblock.h"   // this file
#    include BOOST_PP_ITERATE()

#endif // MEMBLOCK_INCLUDED

#else // BOOST_PP_IS_ITERATING
#   define n BOOST_PP_ITERATION()
#   define nprev BOOST_PP_SUB(n, 1)
#   define MemoryBlockNd BOOST_PP_CAT(MemoryBlock, BOOST_PP_CAT(n, d))
#   define MemoryBlockPREVd BOOST_PP_CAT(MemoryBlock, BOOST_PP_CAT(nprev, d))
#   define MB_print(z, m, data) data
#   define MB_timesDATAM(z, m, data) * BOOST_PP_CAT(data, m)

template <typename T> class MemoryBlockNd : public MemoryBlock1d<T>
{
private:
    MemoryBlockPREVd<T*> m_pointers;
public:
    MemoryBlockNd() { MemoryBlock1d<T>::data() = 0; }
    MemoryBlockNd(BOOST_PP_ENUM_PARAMS(n, int size)) { MemoryBlock1d<T>::data() = 0;  resize(BOOST_PP_ENUM_PARAMS(n, size)); }
    MemoryBlockNd(BOOST_PP_ENUM_PARAMS(n, int size), T val) { MemoryBlock1d<T>::data() = 0;  resize(BOOST_PP_ENUM_PARAMS(n, size), val); }
    inline void resize(BOOST_PP_ENUM_PARAMS(n, int size));
    inline void resize(BOOST_PP_ENUM_PARAMS(n, int size), T val);
    inline T BOOST_PP_REPEAT(nprev, MB_print,*) & operator[](int i) { return m_pointers[i]; }
};

// Implementation n-dim:

template <typename T> void MemoryBlockNd<T>::resize(BOOST_PP_ENUM_PARAMS(n, int size))
{
    m_pointers.resize(BOOST_PP_ENUM_PARAMS(nprev, size));
    int sizePointers = size0 BOOST_PP_REPEAT_FROM_TO(1, nprev, MB_timesDATAM, size);
    int sizeData = sizePointers * BOOST_PP_CAT(size, nprev);
    MemoryBlock1d<T>::resize(sizeData);

    T *p = MemoryBlock1d<T>::data();
    T **ptr = m_pointers.data();

    for (int i=0; i < sizePointers; i++, p += BOOST_PP_CAT(size, nprev)) ptr[i] = p;
}

template <typename T> void MemoryBlockNd<T>::resize(BOOST_PP_ENUM_PARAMS(n, int size), T val)
{
    resize(BOOST_PP_ENUM_PARAMS(n, size));
    int sizeData = size0 BOOST_PP_REPEAT_FROM_TO(1, n, MB_timesDATAM, size);
    std::fill(MemoryBlock1d<T>::data(), MemoryBlock1d<T>::data() + sizeData, val);
}

#  undef MB_timesDATAM
#  undef MB_print
#  undef MemoryBlockPREVd
#  undef MemoryBlockNd
#  undef nprev
#  undef n
#endif // BOOST_PP_IS_ITERATING

memblock.cpp

// test the class here
#include <stdio.h>
#include <memblock.h>

int main() {
  MemoryBlock3d<double> cube(2, 2, 2);
  double *p = &cube[0][0][0];

  for (int i=0; i < 8; i++)
    p[i] = 0.5 + i;

  for (int i=0; i < 2; i++)
    for (int j=0; j < 2; j++)
      for (int k=0; k < 2; k++)
        printf("%f\n", cube[i][j][k]);

  return 0;
}