Generic matrices implementation

Question

I am implementing a litte generic math library. What I have done is to write my generic matrix and vector class. I'm curious if I have it done right so far (implementation wise not totally mathematical correctness wise). Reason: I'm relatively new to template programming and I am not sure what is right, what could be done better and what should I try to avoid at all costs if I use templates in C++.

So far I managed to multiply different matrices of (obviously) different sizes while maintaining static allocation by using only std::array, so no dynamic memory allocation.

One thing is needed to be mentioned: my Vectors and Matrices start count their elements with #1 not #0 like it is usually done in C++ with arrays. The reason for this is to be near as possible to 'proper' math, but I am open to get convinced otherwise.

First of all, my generic Vector class is quite boring but needed for my matrix class, so I put it inside this post.

#include <array>
#include <cmath>

namespace jslmath
{
    template<size_t Dimension,typename NumberType = double>
    class Vector
    {
    public:
        using Value = NumberType;
        using Storage = std::array<NumberType, Dimension>;

    private:
        Storage mField;
    public:
        Vector(const Vector&) = default;
        Vector(Vector&&)  = default;
        virtual ~Vector() = default;

        template<typename ...Targs>
        Vector(Targs... args): mField({args...}){}

        Vector(const Storage& args) : mField(args){}

        Vector& operator=(const Vector&) = default;
        Vector& operator=(Vector&&) = default;

        Value operator[](size_t index)const
        {
            return at(index);
        }
        Value at(size_t index)const
        {
            if (index <= Dimension && index != 0)
                return mField[index - 1];
            throw;//need improvment
        }

        Vector operator+(Vector vec) const{ return add(vec); }
        Vector operator-(Vector vec) const{ return sub(vec); };
        Vector operator/(Vector vec) const{ return div(vec); };
        Vector operator*(Vector vec) const{ return mul(vec); };

        Vector operator+(Value scalar) const{ return add(scalar); };
        Vector operator-(Value scalar) const{ return sub(scalar); };
        Vector operator/(Value scalar) const{ return div(scalar); };
        Vector operator*(Value scalar) const{ return mul(scalar); };

        Vector add(Vector vec)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] + vec[i + 1];
            return {temp};
        }
        Vector sub(Vector vec)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] - vec[i + 1];
            return { temp };
        }
        Vector div(Vector vec)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] / vec[i + 1];
            return { temp };
        }
        Vector mul(Vector vec)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] * vec[i + 1];
            return { temp };
        }


        Value dot(Vector vec) const
        {
            Value tmp = 0;
            for (auto i = 1; i <= Dimension; i++)
                tmp = vec[i] + at(i);
            return tmp;
        }

        Value magnitude() const
        {
            return std::sqrt(dot(*this));
        }
        Value magnitudeSq() const
        {
            return (dot(*this));
        }
        Value distance(Vector vec) const
        {
            auto tmp = vec - *this;
            return tmp.magnitude();
        }
        void normalize()
        {
            *this = *this * (1.0 / magnitude());
        }
        double Angle(Vector vec)
        {
            return std::acos(dot(vec)/ std::sqrt(magnitudeSq()*vec.magnitudeSq()));
        }

        Vector add(Value scalar)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] + scalar;
            return { temp };
        }
        Vector sub(Value scalar)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] - scalar;
            return { temp };
        }
        Vector div(Value scalar)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] / scalar;
            return { temp };
        }
        Vector mul(Value scalar)const
        {
            Storage temp;
            for (auto i = 0; i < Dimension; i++)
                temp[i] = mField[i] * scalar;
            return { temp };
        }

    };

So from here on now the interesting part starts:

    template<size_t N,size_t M,typename NumberType = double>
    class Matrix
    {
    public:
        using Value = NumberType;
        using Storage = std::array<NumberType, N*M>;
        using RowVec = Vector<M, NumberType>;
        using ColVec = RowVec;

The first part of my class is just some lifetime saver and also makes the code more readable:

private:
        Storage     mGrid;

    public:

        Matrix(const Matrix&) = default;
        Matrix(Matrix&&) = default;


        template<typename ...Targs>
        Matrix(Targs... args) : mGrid({ args... }) { }

        Matrix(const Storage& args) : mGrid(args) {  }

        RowVec operator[](size_t index) const
        {
            return Row(index);
        }
        Value operator()(size_t Row, size_t Col) const 
        {
            if(Row != 0 && Col != 0)
                return mGrid[(Col-1) + M*(Row-1)];
            throw;
        }

        constexpr size_t Height() { return N; };
        constexpr size_t Width() { return M; }

        Storage& data() { return mGrid; }

        RowVec Row(size_t index)const
        {
            if (index <= N && index != 0)
            {
                typename RowVec::Storage temp;
                for (auto i = 1; i <= M; i++)
                    temp[i - 1] = (*this)(index, i);
                return RowVec(temp);
            }
            throw;
        }
        RowVec Col(size_t index)const
        {
            if (index <= N && index != 0)
            {
                typename RowVec::Storage temp;
                for (auto i = 1; i <= N; i++)
                    temp[i - 1] = (*this)(i, index);
                return RowVec(temp);
            }
            throw;
        }

        void transpose(){
            Storage tmp;
            for (auto i = 0; i < N * M; i++) {
                int row = i / N;
                int col = i % M;
                tmp[i] = mGrid[M * col + row];
            }
            mGrid = tmp;
        }

        template< template<size_t, size_t, typename>class B, size_t I, size_t J, typename Type>
        auto operator*(B<I, J, Type>& b)
        {
            return mul(b);
        }

        auto operator*(Value val)
        {
            return mul(val);
        }

The most complicated part of my code so far. I have to calculate the size of my new matrix at compile time so that I am able to keep my goal to avoid dynamic allocation. The syntax of template templates is quite odd to me but it does work.

        template< template<size_t,size_t,typename>class B, size_t I, size_t J, typename Type>
        auto mul(B<I,J,Type>& b) -> decltype(Matrix<N, J, Value>{})
        {
            Matrix<N, J, Value> result;
            for (auto i = 0; i < N; ++i)
                for (auto j = 0; j < J; ++j)
                {
                    for (int k = 0; k < I; ++k)
                    {
                        int _a = M * i + k;
                        int _b = J * k + j;
                        result.data()[J * i + j] += this->mGrid[_a] * b.data()[_b];
                    }
                }
            return result;
        }

        Matrix mul(Value b)
        {
            Storage result;
            for (auto i = 0; i < N; ++i)
                for (auto j = 0; j < M; ++j)
                {
                    result[M*i + j] = mGrid[M*i + j] * b;
                }
            return { result };
        }

};

}

And a small test application:

int main(int argc, char* argv[])
{
    jslmath::Matrix<3, 3> Test3( 1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0 );
    jslmath::Matrix<3, 1> Test4( 1.0, 2.0, 3.0 );


    auto T = Test3 * Test4;


    return T(1,1);
}

With the help of an online compiler I already found out that VS is not particularly great at compiling this code example compared to gcc and clang. What totally surprised me was that clang was so good at optimizing my code that in the and all that was left were a single return statement with value 14!

Clang did all the computation of multiplying two different sized matrices at compile time. I am flabbergasted about that. It is also the reason why I posted this code here. I want to know from you what I did well and what could be improved. I'm not really sure how I did that on the first try without even thinking of active optimization from my side.

Answer 1

Some Points

So far I managed to multiply different matrices of (obviously) different sizes while maintaining static allocation by using only std::array, so no dynamic memory allocation.

Sure. If this is something you really want.

One thing is needed to be mentioned: my Vectors and Matrices start count their elements with #1 not #0 like it is usually done in C++ with arrays. The reason for this is to be near as possible to 'proper' math, but I am open to get convinced otherwise.

That depends entirely on the user base who will be using this class. If they are mathematicians who are used to 1 based array then fine. But think that most people who use C like languages are already used to using 0 based arrays so have this may confuse people.

Some Thoughts

A lot of Matrix libraries actually delay the multiplication (and other operations) until the value inside the matrix is required. That way you don't pay for operations that you don't need.

Also by deferring operations you can potentially eliminate null operations or simplify operations that are cumulative.

Example:

 {
     Matrix<4,5>    x(init);
     Matrix<4,5>    y(init);

     Matrix<4,5>    z = x + y;  // Here the + operator
                                // Can loop over all the elements
                                // and do the operation for each element.

     std::cout << z[1][1] << "\n"; // Here we use only one value
 }                                 // Then z goes out of scope and is destroyed.
                                   // So we just did a bunch of operations
                                   // that are not needed.

// If at the point where we did the operation + we returned an
// object that knows about x and y but did not immediately do the operation
// Then we accesses element [1][1] we see the work has not been done
// and just do the operation for that location. We don't work out all
// the elements just the one we want and only when we need it.
// That is a deferred operation.

Example 2:

 Matrix<4,5>    a(init);
 Matrix<5,3>    b(init);
 Matrix<4,3>    c = a * b;
 Matrix<4,3>    d = c;
 Matrix<4,3>    e = c - d;

In this example we can see that all elements of e will be zero. So calculating the value of c first is a waste of time. By deferring calculations of the elements we can sometimes determine the result without having to do all the expensive operations.

So at run-time you can perform operations that cancel out other operations and thus you do not need to perform expensive operations if there results do not generate a value that effects the result.

Code Review

If your Value (or NumberType) is always a simple POD then this is fine.

        template<typename ...Targs>
        Vector(Targs... args): mField({args...}){}

But maybe you want arrays/matrices of complex types then it may be useful to forward the values.

        template<typename ...Targs>
        Vector(Targs&&... args): mField({std::forward<Targs>(args)...}){}

In C++ normally the operator[] is unchecked to allow for optimal speed (you don't need to check the ranges in the function if you have already checked it outside the function) while the method at() does do a range check (because you have not done the check externally).

Your access access operators change this behavior.

        Value operator[](size_t index)const
        {
            // Changes an unchecked operation into a checked operation.
            return at(index);
        }
        Value at(size_t index)const
        {
            // Uses the unchecked access operator[]
            // But does manually checking the mField.at() will check.
            if (index <= Dimension && index != 0)
                return mField[index - 1];
            throw;
        }

The second thing to think about is returning by value. Usually when you have a container you return accesses to the element by reference. Thus allowing you to modify the value intuitively. This also makes writting operator+= simpler (then operator+ can be written in terms of operator+=.

The third thing is that throw; (without an expression is wrong). This is used to re-throw a currently propogating exception from within a catch clause. There is no catch here so this results in a call to std::terminate

        // Const and non const version return by reference
        // Use correct underlying method for checking and non checking.
        Value const& operator[](size_t index) const
        {
            return mField[internalIndex(index)];
        }
        Value&       operator[](size_t index)
        {
            return mField[internalIndex(index)];
        }
        Value const& at(size_t index) const
        {
            return mField.at(internalIndex(index));
        }
        Value&       at(size_t index)
        {
            return mField.at(internalIndex(index));
        }

Creating Operators.

Normally you implement operator+ in terms of operator+=.

        Vector operator+(Vector vec) const{ return add(vec); }

But you should also pass the parameter by const reference to avoid a copy.

        Vector& operator+=(Vector const& vec)
        {
            for (auto i = 1; i <= Dimension; ++i)
                (*this)[i] += vec[i];
            return *this;
        }
        Vector operator+(Vector const& vec) const
        {
            Vector  temp(*this);
            return temp += vec;
        }

Same applies to all the operators.

Matrix definitions

    using RowVec = Vector<M, NumberType>;
    using ColVec = RowVec;

Should ColVec not have a size N?

    using ColVec = Vector<N, NumberType>;

Answer 2

Consider utilizing the empty base class optimization

This will be the hardest to implement as it requires careful refactoring of your class. A typical container is implemented like this:

• It has an allocator member

• It has begin, end and capacity pointer members

Because this bloats the container, it's possible to utilize the EBCO to make the size of these members just a single pointer by either:

• Encapsulating it in a base class that has these as members, then have the derived class inherit it and use it to access the members; typically called something like "impl"

• Wrap it in a boost::compressed_pair which requires less legwork

I'm unsure how useful this will be since you are using an std::array for storage; however, consider this. The size of Matrix is proportional to its geometry (Matrix<3, 3> is 72 bytes on a 64-bit system) even though it is a non-owning view. So further refactoring leads to these possible routes:

• Give Matrix a pointer to the storage owned by someone else (i.e Vector) and a capacity; have Matrix perform all of its operations in-place (rather than using temporary variables)

• Remove implementation details from Matrix that it doesn't "know" about. You hard-code the storage type, but could use using Storage = Vector<...>::Storage

Consider an STL compatible interface

This one is a no-brainer and just requires refactoring your typedef's and interface to look similar to an STL container. Using value_type, reference etc. Make the Storage typedef private since it is an implementation-detail.

Having a constructor that takes an initializer_list might be useful too.

Throw std::out_of_range

This is a freebie. You throw inside at, might as well make it a logical exception type.