4
\$\begingroup\$

I've been learning C++ on and off, and recently learnt about overloading and templates. I wanted to used these concepts to create classes about matrices using std::array.

Note that this has to be compiled using std=c++17 or newer, following a clarification involving constexpr.

#include <array>
#include <iostream>

template <std::size_t R, std::size_t C>
class Matrix {
   protected:
    std::array<std::array<double, C>, R> matrix;

   public:
    /// @brief sets all elements of the matrix to 0
    Matrix() {
        for (std::size_t row = 0; row < R; row++) {
            for (std::size_t column = 0; column < C; column++) {
                matrix[row][column] = 0;
            }
        }
    }

    /// @brief sets the matrix to the given matrix
    /// @param matrix the input data
    Matrix(std::array<std::array<double, C>, R> matrix) : matrix(matrix) {}

    /// @brief adds the given matrix to the current matrix, element by element
    /// @param matrix the input data
    /// @return the sum of the two matrices
    Matrix<R, C> operator+(Matrix<R, C> matrix) {
        Matrix<R, C> returnMatrix;

        for (std::size_t row = 0; row < R; row++) {
            for (std::size_t column = 0; column < C; column++) {
                returnMatrix[row][column] = this->matrix[row][column] + matrix[row][column];
            }
        }

        return returnMatrix;
    }

    /// @brief subtracts the given matrix from the current matrix, element by element
    /// @param matrix the input data
    /// @return the difference of the two matrices
    Matrix<R, C> operator-(Matrix<R, C> matrix) {
        Matrix<R, C> returnMatrix;

        for (std::size_t row = 0; row < R; row++) {
            for (std::size_t column = 0; column < C; column++) {
                returnMatrix[row][column] = this->matrix[row][column] - matrix[row][column];
            }
        }

        return returnMatrix;
    }

    /// @brief multiplies the given matrix with the current matrix
    /// @param matrix the input data
    /// @return the product of the two matrices
    /// @warning the number of columns of the first matrix must be equal to the number of rows of the second matrix. I'm unable to figure out a way to check this, since the number of rows and columns are template parameters
    Matrix<R, C> operator*(Matrix<R, C> matrix) {
        Matrix<R, C> returnMatrix;

        for (std::size_t k = 0; k < R; k++) {
            for (std::size_t i = 0; i < R; i++) {
                for (std::size_t j = 0; j < C; j++) {
                    returnMatrix[i][j] += this->matrix[i][k] * matrix[k][j];
                }
            }
        }

        return returnMatrix;
    }

    /// @brief allows you to index a row of the matrix
    /// @param row index of the row
    /// @return the row of the matrix, as an array
    std::array<double, C>& operator[](std::size_t row) {
        return matrix[row];
    }

    /// @brief return the matrix
    /// @return the `std::array<std::array>` object
    std::array<std::array<double, C>, R> getMatrix() {
        return matrix;
    }

    /// @brief returns the element at the given row and column
    /// @param row index of the row
    /// @param column index of the column
    int getElement(std::size_t row, std::size_t column) {
        return matrix[row][column];
    }

    /// @brief sets the matrix to the given matrix
    /// @param matrix the input data
    void setMatrix(std::array<std::array<double, C>, R> matrix) {
        this->matrix = matrix;
    }

    /// @brief sets the element at the given row and column to the given value
    /// @param row index of the row
    /// @param column index of the column
    void setElement(std::size_t row, std::size_t column, double value) {
        matrix[row][column] = value;
    }

    /// @brief prints the matrix to the console
    void printMatrix() {
        printf("\n[%dx%d] matrix\n", R, C);

        for (std::size_t row = 0; row < R; row++) {
            for (std::size_t column = 0; column < C; column++) {
                printf("%.0f\t", matrix[row][column]);
            }

            printf("\n");
        }

        printf("\n");
    }
};

/// @brief square matrix, a matrix with the same number of rows and columns
/// @tparam N number of rows and columns
template <std::size_t N>
class SquareMatrix : public Matrix<N, N> {
   public:
    /// @brief reuse the constructor of the base class
    SquareMatrix() : Matrix<N, N>() {}

    /// @brief sets the matrix to the given matrix. reuses the constructor
    /// @param matrix  the input data
    SquareMatrix(std::array<std::array<double, N>, N> matrix) : Matrix<N, N>(matrix) {}

    /// @brief extracts the submatrix of the given matrix, ignoring the given row and column
    /// @param ignoreRow row to ignore
    /// @param ignoreColumn column to ignore
    /// @return `SquareMatrix<N - 1>` object, the submatrix
    SquareMatrix<N - 1> subMatrix(std::size_t ignoreRow, std::size_t ignoreColumn) {
        SquareMatrix<N - 1> returnMatrix;

        int subMatrixRow = 0, subMatrixColumn = 0;

        for (std::size_t matrixRow = 0; matrixRow < N; matrixRow++) {
            for (std::size_t matrixColumn = 0; matrixColumn < N; matrixColumn++) {
                if (matrixRow != ignoreRow && matrixColumn != ignoreColumn) {
                    returnMatrix[subMatrixRow][subMatrixColumn++] = this->matrix[matrixRow][matrixColumn];

                    if (subMatrixColumn == N - 1) {
                        subMatrixColumn = 0;
                        subMatrixRow++;
                    }
                }
            }
        }

        return returnMatrix;
    }

    /// @brief calculates the determinant of the matrix. 
    /// see Laplace expansion. https://en.wikipedia.org/wiki/Determinant#Laplace_expansion
    /// @return the determinant of the matrix
    double getDeterminant() {
        double determinant = 0;

        // base case 1
        if constexpr (N == 1) {
            determinant = this->matrix[0][0];
        }

        // base case 2
        else if constexpr (N == 2) {
            determinant = (this->matrix[0][0] * this->matrix[1][1]) - (this->matrix[0][1] * this->matrix[1][0]);
        }

        // recursive case
        else if constexpr (N > 0) {
            int sign = 1;

            for (std::size_t dimension = 0; dimension < N; dimension++) {

                // calculate cofactor and add to determinant
                determinant += sign * this->matrix[0][dimension] * subMatrix(0, dimension).getDeterminant();
                sign = -sign;
            }
        }

        else {
            throw std::invalid_argument("expected square matrix");
        }

        return determinant;
    }
};

int main() {
    static const std::size_t LENGTH = 4;

    SquareMatrix<LENGTH> matrixOne({{{7, -2, 2, 1},
                                     {3, 1, -5, 2},
                                     {2, 2, -5, 3},
                                     {3, -2, 5, 1}}});

    SquareMatrix<LENGTH> matrixTwo({{{1, 2, 3, 4},
                                     {5, 6, 7, 8},
                                     {9, -8, 7, 6},
                                     {5, 4, -3, -2}}});

    matrixOne.printMatrix();
    matrixTwo.printMatrix();

    printf("\nmatrixOne + matrixTwo\n");
    (matrixOne + matrixTwo).printMatrix();

    printf("\nmatrixOne - matrixTwo\n");
    (matrixOne - matrixTwo).printMatrix();

    printf("\nmatrixOne * matrixTwo\n");
    (matrixOne * matrixTwo).printMatrix();

    printf("determinant of matrixOne: %.0f\n", matrixOne.getDeterminant());
    printf("determinant of matrixTwo: %.0f\n", matrixTwo.getDeterminant());
}

Program output

g++ ./determinant-inheritance.cpp -o determinant-inheritance -std=c++20 && ./determinant-inheritance

[4x4] matrix
7       -2      2       1
3       1       -5      2
2       2       -5      3
3       -2      5       1


[4x4] matrix
1       2       3       4
5       6       7       8
9       -8      7       6
5       4       -3      -2


matrixOne + matrixTwo

[4x4] matrix
8       0       5       5
8       7       2       10
11      -6      2       9
8       2       2       -1


matrixOne - matrixTwo

[4x4] matrix
6       -4      -1      -3
-2      -5      -12     -6
-7      10      -12     -3
-2      -6      8       3


matrixOne * matrixTwo

[4x4] matrix
20      -10     18      22
-27     60      -25     -14
-18     68      -24     -12
43      -42     27      24

determinant of matrixOne: 47
determinant of matrixTwo: -640
\$\endgroup\$

2 Answers 2

4
\$\begingroup\$

Avoid C functions

I see you are using printf() in your code. Consider avoiding this and using the C++ way of printing instead. The main disadvantage of printf() is that it is not type safe. The main advantage of printf() is perhaps the format string, but since C++20 there is std::format(), and in C++23 we will get the best of everything with std::print(). If you cannot use a newer C++ version yet, you can link with the {fmt} library to get the same functionality.

Note that with the "old" way of printing in C++, you could have written:

void printMatrix() {
    std::cout << "\n[" << R << 'x' << C << "] matrix\n";

    for (auto& row: matrix) {
        for (auto& value: row) {
            std::cout << std::setprecision(0) << value << '\t';
        }

        std::cout << '\n';
    }

    std::cout << '\n';
}

Note that if you want to keep using C's printf(), you should #include <cstdio> instead of #include <iostream>.

Use I/O operator overloading to print

Wouldn't it be nice if you could print by write the following?

Matrix<…> matrix(…);

std::cout << matrix;

You can quite easily change your printMatrix() function into an overload of operator<<(std::ostream&). The main advantage is that you no longer are hardcoding where the output goes: you can then print to std::cerr, to a std::stringstream, a std::fstream and whatever else provides a std::ostream interface.

Avoid making unnecessary copies of matrices

While your code might be functionally correct, it makes copies of matrices in lots of places. For example, because operator+() takes the parameter matrix by value, whenever it is called a copy is made of that matrix. To avoid that, pass the input by const reference:

Matrix operator+(const Matrix& matrix) {
    …
}

You should not do the same for the return value; copy elision will take of that.

No need to repeat the template parameters inside the template itself

Inside the declaration of class Matrix you don't need to repeat the template parameters if they are the same as the object's, you can just write Matrix and it will be equivalent to Matrix<R, C>. Prefer the short form, there is less chance of making mistakes this way, and if you really do want different template parameters you can still do that, but it will then clearly stand out.

About matrix-matrix multiplications

Your operator*() only allows two matrices of the same size to be multiplied. However, in mathematics, you can multiply two matrices of different size, as long as the number of columns on the left hand matches the number of rows on the right hand of the multiplication. You can easily change your operator to do this. First, you know that the matrix on the right hand side should have the same number of rows as the left hand side has columns, that's just C. But then the right hand side can have it's own number of columns, so let's name that RHC. Then you can just write:

template<std::size_t RHC>
Matrix<R, RHC> operator*(const Matrix<C, RHC>& matrix) {
    …
}

And you just have to change one of the for-loops to use RHC.

Consider writing out-of-class operators

Instead of writing the operators as class member functions, you can also make them out-of-class friend functions:

template<…>
class Matrix {
    …
    friend Matrix operator+(const Matrix& lhs, const Matrix& rhs);
    …
};

template<std::size_t R, std::size_t C>
Matrix<R, C> operator+(const Matrix<R, C>& lhs, const Matrix<R, C>& rhs) {
    Matrix<R, C> result;
    …
        result[row][column] = lhs[row][column] + rhs[row][column];
    …
    return result;
}

You might not see the point of this, but it has several advantages. In particular, consider that you might want to allow multiplying a matrix by a scalar value. You could add an overload of the member function operator*():

Matrix operator*(double value) {…}

That would allow you to write:

Matrix<…> matrix(…);
auto matrix2 = matrix * 2;

But it would not allow you to write:

auto matrix2 = 2 * matrix;

And there is no way to achieve that with a member function. However, the friend operators get two parameters, so then you can do:

friend Matrix operator*(const Matrix& lhs, double rhs);
friend Matrix operator*(double lhs, const Matrix& rhs);

In this case, because the operation is commutative, you only have to fully implement one of them:

template<std::size_t R, std::size_t C>
friend Matrix<R, C> operator*(const Matrix<R, C>& lhs, double rhs) {
    Matrix result = lhs;

    for (auto& row: result)
        for (auto& value: row)
            value *= rhs;

    return result;
}

Then the other can just call the first one with the arguments swapped:

template<std::size_t R, std::size_t C>
friend Matrix<R, C> operator*(double lhs, const Matrix<R, C>& rhs) {
    return rhs * lhs;
}

Apart from the binary operators themselves, you can also make members like getDeterminant() friends instead:

friend double getDeterminant(const Matrix& matrix);

You would then need to call it like getDeterminant(matrix) instead of matrix.getDeterminant(), but that's similar to how you use std::abs() and other standard math functions.

Element access

You have getElement() and setElement() members, but have you noticed that standard containers like std::array and std::vector don't have such functions at all? They have overloads of operator[] and have a member function at() that allows read/write access to a given element. You can do the same:

double& at(std::size_t row, std::size_t column) {
    return matrix[row][column];
}

Because a reference is returned, the caller can then do:

matrix.at(1, 2) = 3;

Also note that your getElement() returned an int instead of a double.

You can also overload operator[]. Before C++23, it was not possible to pass more than one parameter inside the brackets, but you could have cheated by writing something like:

double& operator[](std::pair<std::size_t, std::size_t> index) {
    return matrix[index.first][index.second];
}

And call it like:

matrix[{1, 2}] = 3;

Make member functions const where appropriate

Member functions themselves can be marked const to indicate that they won't modify any member variables. Consider writing:

const Matrix<…> matrix(…);
auto value = matrix.getElement(1, 2); // compile error!

This wouldn't compile, because getElement() is not marked const, so the compiler will not allow you to call it on a const object. The simple solution is to make getElement() const:

double getElement(std::size_t row, std::size_t column) const {
   …
}

But that works here because this function returns by value, and never has to modify any members. But the above at() has a problem: if you want to be able to use it to write to an element, it needs to return a non-const reference. But that is not allowed if you apply it to a const object. You can however create a second overload for const objects:

double& at(std::size_t row, std::size_t column) {
    return matrix[row][column];
}

const double& at(std::size_t row, std::size_t column) const {
    return matrix[row][column];
}

You'll note that the body of those functions are the same, which is annoying. There are two ways around it; having one call the other but using std::const_cast to cheat. Another way is to use C++23's explicit object parameter.

About inheritance

Inheritance is one way to create a SquareMatrix class. However, it has some drawbacks, and there are alternatives you could use. Consider for example making getDeterminant() an out-of-class function that takes a SquareMatrix as a parameter. Now I can write:

SquareMatrix<…> square;
auto det = getDeterminant(square);

But now consider I wrote:

Matrix<3, 3> square;
auto det = getDeterminant(square); // compile error

This fails to compile because getDeterminant() only works on SquareMatrixes, even though Matrix<3, 3> is also square. You could create an overload for getDeterminant() to take square Matrix objects as well, but that's more work and extra code, with more chance of bugs.

Instead of using inheritance, consider creating a type alias:

template<std::size_t N>
using SquareMatrix = Matrix<N, N>;

Your getDeterminant() function will then look like:

template<std::size_t N>
double getDeterminant(const SquareMatrix<N>& matrix) {…}

But SquareMatrix was just an alias, so I can pass in a Matrix<3, 3> and it will still work. If I pass in a Matrix<4, 5> it will fail to compile.

Another option would be to use static_assert(), SFINAE or C++20's concepts to ensure getDeterminant() only compiles for square matrices. For example:

template<std::size_t R, std::size_t C>
double getDeterminant(const Matrix<R, C>& matrix) {
    static_assert(R == C, "The determinant is only defined for square matrices.");
    …
}

Or with C++20:

template<std::size_t R, std::size_t C>
requires (R == C)
double getDeterminant(const Matrix<R, C>& matrix) {
    …
}

In this simple case you can also just make sure template deduction only works for square matrices:

template<std::size_t N>
double getDeterminant(const Matrix<N, N>& matrix) {…}

But in fact I already showed that, because this is equivalent to writing getDeterminant(const SquareMatrix<N>& matrix).

\$\endgroup\$
3
  • \$\begingroup\$ Thanks a lot for your insights! Could you please clarify a little about the subheading Consider writing out-of-class operators? In particular, I don't understand how you can "define one overload in terms of the other" for the scalar-matrix multiplication. You have to write the detailed routine at least once, right? \$\endgroup\$ Jun 22, 2023 at 13:58
  • \$\begingroup\$ Yes, you have to do that. I amended the example. \$\endgroup\$
    – G. Sliepen
    Jun 22, 2023 at 14:28
  • \$\begingroup\$ Gotcha. Thanks a lot again for the detailed feedback. If you're interested, here's the issue incorporating it all into my code: github.com/eccentricOrange/c-cpp-mirror/issues/1 \$\endgroup\$ Jun 22, 2023 at 18:21
2
\$\begingroup\$

After G. Sliepen’s answer, there is not much left to comment on. I’ll point out one simplification that I think is worth while: instead of

std::array<std::array<double, C>, R> matrix;

use

std::array<double, C * R> matrix;

and then compute your indices to the various matrix elements. So instead of matrix[row][column] write matrix[row * C + column].

std::array is static, and so an array of arrays has, in fact, contiguous memory. When doing the same with std::vector, one would produce a very inefficient container because it’s not contiguous. Still, even though it’s contiguous, there’s more writing involved. For example the initializer that sets all elements to zero needs a double loop. With the flat memory layout, it’d just be

for (auto& value : matrix) { value = 0; }

(Or you could use std::fill). The addition and subtraction functions will be similarly simplified.

You will have to explicitly compute the linear index from the two indices, but you can delegate this to a private function, so the code will not look any more complex. For example you could write a function element(row, column). The cost is the same, with the array of arrays the compiler computes the index in the same way. (With a vector of vectors there’s a double dereference, which is part of why that is so expensive.)

You can now quite easily write an iterator that visits elements in turn. One version of the begin function would return the iterator set up to iterate over one row, another to iterate over one column, another over the diagonal, another over all elements of the matrix, etc. All this iterator needs is the index to the current element, and a step size, how much to increment the index to get to the next element in the row/column/diagonal/etc. Doing the same with an array of arrays is quite challenging (unless you want to cheat and just take the pointer to the first element as the pointer to a contiguous buffer, which I guess is possible but ugly, and defeats the purpose of the array of arrays).

\$\endgroup\$
1
  • 1
    \$\begingroup\$ Indeed. Also, this is a good place to start using std::mdspan if you can use C++23. \$\endgroup\$
    – G. Sliepen
    Jun 22, 2023 at 14:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.