Making a function that can accept different types of templated class inputs

I have an application where I need to reduce the memory usage of my matrix classes. Many of my matrices are symmetric, and so should only require N(N+1)/2 memory instead of N2. I can’t simply have a MatrixSym class that derives from my Matrix class because a child class always uses more memory than the base class. I’ve defined a MatrixAbstract class from which both classes inherit, however I am having trouble getting a matrix multiplication function to work on both of them.

My question

How can I define a C = multiply(A,B) function that accepts either a Matrix or MatrixSym object for any input/output? My minimal working example is not so minimal, because it contains four definitions of multiply that all use the same code and should be combined into one.

Requirements

  • Static allocation only
  • Must use std::array for matrix elements
  • Minimal repeated code: I have many matrix operations and the only thing different about MatrixSym is the way it stores and accesses elements from the underlying matVals array.
#include <array> 
#include <iostream> 

template<unsigned int ROWS, unsigned int COLS, unsigned int NUMEL>
class MatrixAbstract
{
private:
    static constexpr unsigned int numel = NUMEL; 
    std::array<double, NUMEL> matVals;   // The array of matrix elements
public:
    MatrixSuperclass(){}
    virtual unsigned int index_from_rc(const unsigned int& row, const unsigned int& col) const = 0;
    // get the value at a given row and column
    double get_value(const int& row, const int& col) const
    {
        return this->matVals[this->index_from_rc(row,col)];
    }
    // set the value, given the row and column
    void set_value(const int& row, const int& col, double value)
    {
        this->matVals[this->index_from_rc(row,col)] = value;
    }

};

template<unsigned int ROWS, unsigned int COLS, unsigned int NUMEL = ROWS*COLS>
class Matrix : public MatrixSuperclass<ROWS, COLS, NUMEL>
{
public:
    Matrix(){}
    // get the linear index in matVals corresponding to a row,column input
    unsigned int index_from_rc(const unsigned int& row, const unsigned int& col) const {
        return row*COLS + col; 
    }
};

template<unsigned int ROWS, unsigned int COLS = ROWS, unsigned int NUMEL = ROWS*(ROWS+1)/2> 
class MatrixSym : public MatrixSuperclass<ROWS, COLS, NUMEL>
{
public:
    MatrixSym(){}
    // get the linear index in matVals corresponding to a row,column input (Symmetric matrix)
    unsigned int index_from_rc(const unsigned int& row, const unsigned int& col) const {
        unsigned int z;
        return ( ( z = ( row < col ? col : row ) ) * ( z + 1 ) >> 1 ) + ( col < row ? col : row ) ;
    }
};

// THE FOLLOWING FOUR FUNCTIONS ALL USE THE EXACT SAME CODE, ONLY INPUT/OUTPUT TYPES CHANGE

// Multiply a Matrix and Matrix and output a Matrix
template<unsigned int ROWS, unsigned int COLS, unsigned int INNER>
Matrix<ROWS,COLS> multiply (Matrix<ROWS,INNER>& inMatrix1, Matrix<INNER,COLS>& inMatrix2) {
    Matrix<ROWS,COLS> outMatrix;
    for (unsigned int r = 0; r < ROWS; r++) {
        for (unsigned int c = 0; c < COLS; c++) {
            double val = 0.0;
            for (unsigned int rc = 0; rc < INNER; rc++) {
                val += inMatrix1.get_value(r,rc)*inMatrix2.get_value(rc,c);
            }
            outMatrix.set_value(r,c,val);
        }
    }
    return outMatrix;
}

// Multiply a Matrix and MatrixSym and output a Matrix
template<unsigned int ROWS, unsigned int COLS, unsigned int INNER>
Matrix<ROWS,COLS> multiply (Matrix<ROWS,INNER>& inMatrix1, MatrixSym<INNER,COLS>& inMatrix2) {
    Matrix<ROWS,COLS> outMatrix;
    for (unsigned int r = 0; r < ROWS; r++) {
        for (unsigned int c = 0; c < COLS; c++) {
            double val = 0.0;
            for (unsigned int rc = 0; rc < INNER; rc++) {
                val += inMatrix1.get_value(r,rc)*inMatrix2.get_value(rc,c);
            }
            outMatrix.set_value(r,c,val);
        }
    }
    return outMatrix;
}

// Multiply a MatrixSym and Matrix and output a Matrix
template<unsigned int ROWS, unsigned int COLS, unsigned int INNER>
Matrix<ROWS,COLS> multiply (MatrixSym<ROWS,INNER>& inMatrix1, Matrix<INNER,COLS>& inMatrix2) {
    //MatrixSym<ROWS,COLS> outMatrixSym;
    Matrix<ROWS,COLS> outMatrix;
    for (unsigned int r = 0; r < ROWS; r++) {
        for (unsigned int c = 0; c < COLS; c++) {
            double val = 0.0;
            for (unsigned int rc = 0; rc < INNER; rc++) {
                val += inMatrix1.get_value(r,rc)*inMatrix2.get_value(rc,c);
            }
            outMatrix.set_value(r,c,val);
        }
    }
    return outMatrix;
}

// Multiply a MatrixSym and MatrixSym and output a MatrixSym
template<unsigned int ROWS, unsigned int COLS, unsigned int INNER>
MatrixSym<ROWS,COLS> multiply (MatrixSym<ROWS,INNER>& inMatrix1, MatrixSym<INNER,COLS>& inMatrix2) {
    //MatrixSym<ROWS,COLS> outMatrixSym;
    MatrixSym<ROWS,COLS> outMatrix;
    for (unsigned int r = 0; r < ROWS; r++) {
        for (unsigned int c = 0; c < COLS; c++) {
            double val = 0.0;
            for (unsigned int rc = 0; rc < INNER; rc++) {
                val += inMatrix1.get_value(r,rc)*inMatrix2.get_value(rc,c);
            }
            outMatrix.set_value(r,c,val);
        }
    }
    return outMatrix;
}

int main()
{
    Matrix<3,3> A;
    MatrixSym<3> S;
    Matrix<3,3> AtimesA = multiply(A,A);
    Matrix<3,3> AtimesS = multiply(A,S);
    Matrix<3,3> StimesA = multiply(S,A);
    MatrixSym<3> StimesS = multiply(S,S);
    
    // Make sure that symmetric matrix S is indeed smaller than A
    std::cout << "sizeof(A)/sizeof(double) = " << sizeof(A)/sizeof(double) << std::endl;
    std::cout << "sizeof(S)/sizeof(double) = " << sizeof(S)/sizeof(double) << std::endl;
    return 0;
}

Outputs:

sizeof(A)/sizeof(double) = 9
sizeof(S)/sizeof(double) = 15

What I’ve tried

  • If I try to make the function use MatrixAbstract as arguments I need to play around with the template parameters, since I can’t have NUMEL as a template parameter. Further, I can’t instantiate a MatrixAbstract for the return value.
  • I can’t figure out how to template the function for either Matrix or MatrixSym. My hunch is that this is the way to solve it, but I don’t understand how to template the function for a templated class input, but with still being able to use ROWS and COLS as template arguments.
template<class InMatrix1Type, class InMatrix2Type, class OutMatrixType, unsigned int ROWS, unsigned int COLS, unsigned int INNER>
OutMatrixType<ROWS,COLS> multiply (InMatrix1Type<ROWS,INNER>& inMatrix1, InMatrix2Type<INNER,COLS>& inMatrix2)

gives me a compiler error starting with

error: ‘OutMatrixType’ is not a template
 OutMatrixType<ROWS,COLS> multiply (InMatrix1Type<ROWS,INNER>& inMatrix1, InMatrix2Type<INNER,COLS>& inMatrix2)
 ^~~~~~~~~~~~~

Answer

The function multiply can be a template that takes two (possibly different) argument types, provided they support the appropriate interface. Like this:

template <unsigned ROWS, unsigned COLS>
struct Matrix {
    double get_value(int row, int col) const;
    void set_value(int row, int col, double value);
};

template <unsigned ROWS, unsigned COLS>
struct Symmetric_Matrix {
    double get_value(int row, int col) const;
    void set_value(int row, int col, double value);
};

template <unsigned ROWS, unsigned COLS, unsigned INNER,
    template <unsigned, unsigned> class M1,
    template <unsigned, unsigned> class M2>
Matrix<ROWS, COLS> multiply(const M1<ROWS, INNER>& m1,
    const M2<INNER, COLS>& m2) {
    Matrix<ROWS, COLS> result;
    for (unsigned r = 0; r < ROWS; ++r)
        for (unsigned c = 0; c < COLS; ++c) {
            double val = 0.0;
            for (unsigned rc = 0; rc < INNER; ++rc)
                val += m1.get_value(r, rc) * m2.get_value(rc, c);
        result.set_value(r, c, val);
        }
    return result;
}

M1 and M2 are template template parameters; that is, they’re templates that are used as template arguments to multiply; the compiler will figure out their types, and the corresponding values of ROW, COL, and INNER when the function is called.

int main() {
    Matrix<3, 5> res;

    Matrix<3, 4> x1;
    Matrix<4, 5> x2;
    res = multiply(x1, x2);

    Symmetric_Matrix<3, 4> sx1;
    Symmetric_Matrix<4, 5> sx2;
    res = multiply(sx1, sx2);

    res = multiply(x1, sx2);
    res = multiply(sx1, x2);

    return 0;
}

Of course, real code would provide implementations for get_value and set_value, but the multiply code doesn’t care how those are implemented, so having two different types will work just fine.