Optimizing matrix multiplication algorithm I have implemented my own version of algorithm to multiply two matrices together and I need someone who know what they are doing to see if there is anything that can be done in a more optimal way. I am also keen on understanding how I can make it robust such that it doesn't crash when non-rectangular matrices are passed as argument i.e. matrix that has uneven number of elements in different Vds it holds.
I am particularly interested in function matrixDot in the code bellow, everything else is only to show how I am using it in my project).
#include &quot;iostream&quot;
#include &lt;vector&gt;

#define LOG(m) std::cout &lt;&lt; m &lt;&lt; std::endl

struct Vd
{
    std::vector&lt;double&gt; v;
};

struct Md
{
    std::vector&lt;Vd&gt; m;

    //fill matrix with num
    void fill(unsigned const int rows, unsigned const int cols, const double num)
    {
        m.clear();
        for (unsigned int i = 0; i &lt; rows; i++)
        {
            Vd tempVec;
            for (unsigned int j = 0; j &lt; cols; j++)
            {
                tempVec.v.push_back(num);
            }
            m.push_back(tempVec);
        }
    }

    friend std::ostream&amp; operator &lt;&lt; (std::ostream&amp; out, const Md&amp; mat)
    {
        out &lt;&lt; &quot;[&quot; &lt;&lt; std::endl &lt;&lt; std::endl;
        for (unsigned int i = 0; i &lt; mat.m.size(); i++)
        {
            out &lt;&lt; &quot;[&quot;;
            for (unsigned int j = 0; j &lt; mat.m[i].v.size(); j++)
            {
                if (j % mat.m[i].v.size() == mat.m[i].v.size() - 1)
                    out &lt;&lt; mat.m[i].v[j] &lt;&lt; &quot;]&quot; &lt;&lt; std::endl &lt;&lt; std::endl;
                else
                    out &lt;&lt; mat.m[i].v[j] &lt;&lt; &quot;, &quot;;
            }
        }
        out &lt;&lt; &quot;]&quot; &lt;&lt; std::endl;
        return out;
    }
};

inline void matrixDot(const Md&amp; m1, const Md&amp; m2, Md&amp; outm)
{
    if (m1.m[0].v.size() &amp;&amp; m2.m.size())
        if (m1.m[0].v.size() != m2.m.size())
        {
            LOG(&quot;Shape mismatch: &quot; &lt;&lt; &quot;matrix1 columns: &quot; &lt;&lt; m1.m[0].v.size() &lt;&lt; &quot;, &quot; &lt;&lt; &quot;matrix2 rows: &quot; &lt;&lt; m2.m.size());
            throw std::exception();
        }

    unsigned int m1x = 0; unsigned int m1y = 0; unsigned int m2y = 0; //m2x = m1y

    while (outm.m.size() &lt; m1.m.size())
    {
        Vd tempv;
        while (tempv.v.size() &lt; m2.m[0].v.size())
        {
            double total = 0.0;
            while (m1x &lt; m1.m[0].v.size())
            {
                total += m1.m[m1y].v[m1x] * m2.m[m1x].v[m2y];
                m1x++;
            }
            tempv.v.push_back(total);
            m1x = 0;
            m2y &lt; m2.m[0].v.size() - 1 ? m2y++ : m2y = 0;
        }
        m1y &lt; m1.m.size() - 1 ? m1y++ : m1y = 0;
        outm.m.push_back(tempv);
    }
}

int main()
{
    Md mat1;
    mat1.fill(5, 2, 1.0);
    Md mat2;
    mat2.fill(2, 6, 2.0);
    Md mat3;
    matrixDot(mat1, mat2, mat3);
    std::cout &lt;&lt; mat3;
}
I have to admit I'm a bit confused, you did the hard parts and have a question on the easy parts? I might be misunderstanding your question.

I am also keen on understanding how I can make it robust such that it doesn't crash when non-rectangular matrices are passed as argument i.e. matrix that has uneven number of elements in different Vds it holds.

You could validate that you have a well formed matrix by
inline bool isValid(const Md&amp; mat)
{
    if (mat.m.size())
    {
        int size = mat.m[0].v.size();
        for (int i = 1; i &lt; mat.m.size(); i++) {
            if (size != mat.m[i].v.size())
            {
                return false;
            }
        }    
    }
    
    return true;
}


and incorporating it into the matrixDot function for validation similar to the shape validation you have right now
    if (m1.m[0].v.size() &amp;&amp; m2.m.size())
        if (m1.m[0].v.size() != m2.m.size())
        {
            LOG(&quot;Shape mismatch: &quot; &lt;&lt; &quot;matrix1 columns: &quot; &lt;&lt; m1.m[0].v.size() &lt;&lt; &quot;, &quot; &lt;&lt; &quot;matrix2 rows: &quot; &lt;&lt; m2.m.size());
            throw std::exception();
        }
        
    if (!isValid(m1))
    {
        LOG(&quot;Invalid matrix :: &quot; &lt;&lt; std::endl);
        std::cout &lt;&lt; m1;
        throw std::exception();
    }
    
    if (!isValid(m2))
    {
        LOG(&quot;Invalid matrix :: &quot; &lt;&lt; std::endl);
        std::cout &lt;&lt; m2;
        throw std::exception();
    }

Any optimizations I can think of are around utilizing std::array instead of std::vector as you are aware of row and column lengths at time of creation.