C++ Matrix Class

nota bene.

This answer has 20 upvotes now, but it is not intended as an endorsement of std::valarray.

In my experience, time is better spent installing and learning to use a full-fledged math library such as Eigen. Valarray has fewer features than the competition, but it isn't more efficient or particularly easier to use.

If you only need a little bit of linear algebra, and you are dead-set against adding anything to your toolchain, then maybe valarray would fit. But, being stuck unable to express the mathematically correct solution to your problem is a very bad position to be in. Math is relentless and unforgiving. Use the right tool for the job.

The standard library provides std::valarray<double>. std::vector<>, suggested by a few others here, is intended as a general-purpose container for objects. valarray, lesser known because it is more specialized (not using "specialized" as the C++ term), has several advantages:

It does not allocate extra space. A vector rounds up to the nearest power of two when allocating, so you can resize it without reallocating every time. (You can still resize a valarray; it's just still as expensive as realloc().)
You may slice it to access rows and columns easily.
Arithmetic operators work as you would expect.

Of course, the advantage over using C is that you don't need to manage memory. The dimensions can reside on the stack, or in a slice object.

std::valarray<double> matrix( row * col ); // no more, no less, than a matrix
matrix[ std::slice( 2, col, row ) ] = pi; // set third column to pi
matrix[ std::slice( 3*row, row, 1 ) ] = e; // set fourth row to e

C++ is mostly a superset of C. You can continue doing what you were doing.

That said, in C++, what you ought to do is to define a proper Matrix class that manages its own memory. It could, for example be backed by an internal std::vector, and you could override operator[] or operator() to index into the vector appropriately (for example, see: How do I create a subscript operator for a Matrix class? from the C++ FAQ).

To get you started:

class Matrix
{
public:
    Matrix(size_t rows, size_t cols);
    double& operator()(size_t i, size_t j);
    double operator()(size_t i, size_t j) const;

private:
    size_t mRows;
    size_t mCols;
    std::vector<double> mData;
};

Matrix::Matrix(size_t rows, size_t cols)
: mRows(rows),
  mCols(cols),
  mData(rows * cols)
{
}

double& Matrix::operator()(size_t i, size_t j)
{
    return mData[i * mCols + j];
}

double Matrix::operator()(size_t i, size_t j) const
{
    return mData[i * mCols + j];
}

(Note that the above doesn't do any bounds-checking, and I leave it as an exercise to template it so that it works for things other than double.)

You could do it that way. The only difference is you'd need to cast the result from malloc.

Rather, you would use a vector, either as a 1D array with computed indexing or an embedded vector. (The former matches your code better.)

For example:

template <typename T> // often, they are templates
struct matrix
{
    // should probably be hidden away, and the class would
    // provide `at` and `operator()` for access
    int col, row;
    std::vector<T> data;

    matrix(int columns, int rows) :
    col(columns), row(rows), 
    data(col * row)
    {}

}

matrix m(4, 4);
m.data[1 + 1 * 4] = /* ... */;

Or:

template <typename T>
struct matrix
{
    int col, row;
    std::vector<std::vector<T> > data;

    matrix(int columns, int rows) :
    col(columns), row(rows), 
    data(col, std::vector(row))
    {}
}

matrix m(4, 4);
m.data[1][1] = /* ... */;

But these are only examples. You'd want to make a full-fledged class; if you want more advice on that, edit your question and clarify you'd like to know the canonical way of implementing matrix classes.

There are pre-existing matrix classes. My favorite is that from boost, UBLAS.

There's lots of subtleties in setting up an efficient and high quality matrix class. Thankfully there's several good implementations floating about.

Think hard about whether you want a fixed size matrix class or a variable sized one. i.e. can you do this:

// These tend to be fast and allocated on the stack.
matrix<3,3> M;

or do you need to be able to do this

// These tend to be slower but more flexible and partially allocated on the heap 
matrix M(3,3);

There's good libraries that support either style, and some that support both. They have different allocation patterns and different performances.

If you want to code it yourself, then the template version requires some knowledge of templates (duh). And the dynamic one needs some hacks to get around lots of small allocations if used inside tight loops.

You could use a template like :

#include <iostream>
using std::cerr;
using std::endl;

//qt4type
typedef unsigned int quint32;

template <typename T>
void deletep(T &) {}
template <typename T>
void deletep(T* & ptr) {
    delete ptr;
    ptr = 0;
}
template<typename T>
class Matrix {
    public:
        typedef T value_type;
        Matrix() : _cols(0), _rows(0), _data(new T[0]), auto_delete(true) {};
        Matrix(quint32 rows, quint32 cols, bool auto_del = true);

        bool exists(quint32 row, quint32 col) const;
        T & operator()(quint32 row, quint32 col);
        T operator()(quint32 row, quint32 col) const;
        virtual ~Matrix();

        int size() const { return _rows * _cols; }
        int rows() const { return _rows; }
        int cols() const { return _cols; }
    private:
        Matrix(const Matrix &);
        quint32 _rows, _cols;
        mutable T * _data;
        const bool auto_delete;
};
template<typename T>
Matrix<T>::Matrix(quint32 rows, quint32 cols, bool auto_del) : _rows(rows), _cols(cols), auto_delete(auto_del) {
    _data = new T[rows * cols];
}
template<typename T>
inline T & Matrix<T>::operator()(quint32 row, quint32 col) {
    return _data[_cols * row + col];
}
template<typename T>
inline T Matrix<T>::operator()(quint32 row, quint32 col) const {
    return _data[_cols * row + col];
}

template<typename T>
bool Matrix<T>::exists(quint32 row, quint32 col) const {
    return (row < _rows && col < _cols);
}

template<typename T>
Matrix<T>::~Matrix() {
    if(auto_delete){
        for(int i = 0, c = size(); i < c; ++i){
            //will do nothing if T isn't a pointer
            deletep(_data[i]);
        }
    }
    delete [] _data;
}

int main() {
    Matrix< int > m(10,10);
    quint32 i = 0;
    for(int x = 0; x < 10; ++x) {
        for(int y = 0; y < 10; ++y, ++i) {
            m(x, y) = i;
        }
    }
    for(int x = 0; x < 10; ++x) {
        for(int y = 0; y < 10; ++y) {
            cerr << "@(" << x << ", " << y << ") : " << m(x,y) << endl;
        }
    }
}

*edit, fixed a typo.

C++ Matrix Class

nota bene.

Related

Recent Posts