Porting Matlab arrays to C/C++

Question

I am porting a matlab program to C/C++. I have several problems with it but one of the most important ones is here: Matlab treats arrays with any dimension the same. Assume we have a function like this,

function result = f(A, B, C)
result = A + 2 * B + C;

A, B, and C can be arrays of any dimensions/size. I am not a C/C++ pro, but I guess it is not a simple & clean job in C. One idea is to use void pointers to pass the arrays to the function. What should I do with the dimensions and array operations (+/*) then? The other idea is to use C++ classes. I might be able to write a template class with all the required array operations such as (*, +, <<, >>, ...). But I am sure it's going to be an exhausting job. Does anybody have a better idea? Any simple/multidimensional/single header file/opensource array class that supports + and * operators?

Answer 1

You could have a look at the boost::ublas library. It supports vectors, matrices, linear algebra etc.

Answer 2

I'd recommend using Armadillo. The docs even have a syntax conversion table for Matlab users.

Answer 3

I would recommend you to use C++ classes. Of course do not implement those classes yourself. There are people out there doing an excellent job in such a task. @Darrem suggested boos::ublas. I would recommend eigen. It much more functionality that ublas, it is well maintained, it supports a wide amount of compilers and its performance it is close and sometimes better than proprietary libraries like Intel MKL.

Answer 4

Since non of the answers worked for me, I made my own class. The class works well for me. But I am pretty sure it is not designed so C++'y. It also may result in memory leaks. If you found a problem with it, please leave me a comment. I'd upgrade the code if I made a new version to keep the answer useful for others. By the way, the class works for 1, 2, and 3 dimensions now, but it can be easily extended to any number of dimensions. Here is the code,

#ifndef XARRAY_H_INCLUDED
#define XARRAY_H_INCLUDED

#include <string>
#include <sstream>
#include <vector>
#include <assert.h>

using namespace std;

template <class T = double>
class XArray
{
    // Fields (keep data)
    int index_helper[10];
    // cells of the array
    vector<T> table;
    // dimensions of the array
    vector<int> dims;

public:
    XArray(){}

    XArray(unsigned int n, int *d)
    {
        dims.resize(n);
        int size = 1;
        for (unsigned int i = 0; i < n; i++) {
            size *= d[i];
            dims[i] = d[i];
        }
        table.resize(size);
    }

    XArray(unsigned int d1)
    {
        dims.resize(1);
        dims[0] = d1;
        table.resize(d1);
    }

    XArray(unsigned int d1, unsigned int d2)
    {
        dims.resize(2);
        dims[0] = d1;
        dims[1] = d2;
        table.resize(d1 * d2);
    }

    XArray(unsigned int d1, unsigned int d2, unsigned int d3)
    {
        dims.resize(3);
        dims[0] = d1;
        dims[1] = d2;
        dims[2] = d3;
        table.resize(d1 * d2 * d3);
    }

    XArray(const XArray<T>& xa)
    {
        this->table = xa.table;
        this->dims = xa.dims;
    }

    int dim(int i)
    {
        return dims[i];
    }

    int num_dims()
    {
        return dims.size();
    }

    T& operator()(int i)
    {
        index_helper[0] = i;
        return get_helper(1, index_helper);
    }

    T& operator()(int i, int j)
    {
        index_helper[0] = i;
        index_helper[1] = j;
        return get_helper(2, index_helper);
    }

    T& operator()(int i, int j, int k)
    {
        index_helper[0] = i;
        index_helper[1] = j;
        index_helper[2] = k;
        return get_helper(3, index_helper);
    }

    XArray<T> operator*(double m)
    {
        XArray<T> r = *this;
        for (unsigned int i = 0; i < table.size(); i++) {
            r.table[i] *= m;
        }

        return r;
    }

    XArray<T> operator/(double m)
    {
        XArray<T> r = *this;
        for (unsigned int i = 0; i < table.size(); i++) {
            r.table[i] /= m;
        }

        return r;
    }

    XArray<T> operator+(const XArray<T> &that)
    {
        assert(this->dims.size() == that.dims.size());
        for (unsigned int i = 0; i < dims.size(); i++) {
            assert(this->dims[i] == that.dims[i]);
        }

        XArray<T> r = *this;
        for (unsigned int i = 0; i < table.size(); i++) {
            r.table[i] += that.table[i];
        }

        return r;
    }

    XArray<T> operator-(const XArray<T> &that)
    {
        assert(this->dims.size() == that.dims.size());
        for (unsigned int i = 0; i < dims.size(); i++) {
            assert(this->dims[i] == that.dims[i]);
        }

        XArray<T> r = *this;
        for (unsigned int i = 0; i < table.size(); i++) {
            r.table[i] -= that.table[i];
        }

        return r;
    }
private:
    T& get_helper(unsigned int n, int *indices)
    {
        assert(n == dims.size());

        int multiplier = 1;
        int index = 0;

        for (unsigned int i = 0; i < n; i++) {
            //cerr << "index " << i << " out of range. Expected [0, " << dims[i] - 1
            //     << "] found " << indices[i] << endl;
            assert(indices[i] >= 0 && indices[i] < dims[i]);
            index += indices[i] * multiplier;
            multiplier *= dims[i];
        }

        return table[index];
    }
};

template <class T>
ostream &operator<<(ostream &stream, XArray<T> xa)
{
    int d = xa.num_dims();

    if(d == 1)
    {
        int n = xa.dim(0);
        for(int i = 0; i < n; i++)
        {
            stream << xa(i);
            if(i < n - 1)
            {
                stream << ", ";
            }
        }
    }
    else
    {
        stream << "XArray[";
        for(int i = 0; i < d; i++)
        {
            stream << xa.dim(i);
            if(i < d - 1)
            {
                stream << "x";
            }
        }
        stream << "]";
    }

    return stream;
}

#endif // XARRAY_H_INCLUDED