C++ Matrix Multiplication Using Classes

Question

I'm trying to write a C++ multiplication program using classes (for the first time - I usually try to avoid OOP) and am having trouble freeing my memory. What I'd ideally like is a matrix class with a function for allocating, initializing, multiplying, and freeing my matrices. What I currently have is

#include <iostream>
#include <cstdio>
#include <ctime>
#include <stdlib.h>

using namespace std;

class matrix 
{
   static const int nrows = 40000;
   static const int ncols = 40000;

   // Allocate matrices for multiplication
   double (*matrix_a)[ncols] = (double (*)[ncols])malloc(sizeof(double) * (nrows) * (ncols));
   double (*matrix_b)[ncols] = (double (*)[ncols])malloc(sizeof(double) * (nrows) * (ncols));
   double (*matrix_c)[ncols] = (double (*)[ncols])malloc(sizeof(double) * (nrows) * (ncols));

   public:
      matrix() 
      {
         // Initialize matrices
         for (int i = 0; i < nrows; i++)
         {
            for (int j = 0; j < nrows; j++)
            {
               matrix_a[i][j] = (double)rand() / (double)RAND_MAX;
               matrix_b[i][j] = (double)rand() / (double)RAND_MAX;
            }
         }
      }

      void multiplication() 
      {
         // Multiply matrices
         for (int i = 0; i < nrows; i++) 
         {
            for (int j = 0; j < ncols; j++) 
            {
               matrix_c[i][j] = 0.0;
               for (int k = 0; k < ncols; k++) 
               {
                  matrix_c[i][j] += matrix_a[i][k] * matrix_b[k][j];
               }
            }
         }
      }

      void free()
      {
         free(matrix_a);
         free(matrix_b);
         free(matrix_c);
      }
};

int main() 
{
   matrix matMul = matrix();
   matMul.multiplication();
   matMul.free();
   return 0;
}

This currently works if you cut out the free() public function (and call to it within main), but I'd like to have this more like my intended description. If anyone could help me mold this into a more close format to my desired program, it would be much appreciated! Thanks in advance.

Answer 1

I will be assuming that you use at least C++11. It also seems that you are coming from C, is that correct?

1. double (*matrix_a)[ncols] = (double (*)[ncols])malloc(sizeof(double) * (nrows) * (ncols));:

   There is no reason to use dynamic allocation, ncols and nrows are constant expressions!. Just use double matrix_a[ncols][nrows]; or better std::array<std::array<double,ncols>,nrows> matrix_a;. That said you usually want to have dynamic sizes or your matrix may be larger than the stack limit. Then you should use std::vector<std::vector<double>> matrix_a{nrows, {ncols}}. C-style arrays and dynamic allocation should always be limited as much as possible. Using proper class objects wrapping them is much safer. std::vector will take care of all memory management for you. If you are concerned that std::vector will layout rows non-continuously, then use std::vector<double> matrix_a{nrows*ncols} instead and properly access the elements by ncols*row+col.

   Also, in C++, you do not use malloc and free. Instead you use new and delete. The difference is that new not only allocates memory but also constructs the object in that memory. In most cases you would need to manually construct the object in the allocated space (although for double it is technically ok here).

2. You should not declare a free method, but rather everything freeing memory belongs in the destructor, which is automatically called as soon as the object itself goes out of scope or is destroyed:

   ~matrix()
   {
        free(matrix_a); // should be "delete[] matrix_a;"
        free(matrix_b); // ^^^
        free(matrix_c); // ^^^
   }
   

   Note that you would not need a destructor if you were using std::array or std::vector.

   If you stick to dynamic memory allocation, you probably also need to implement a custom copy and/or move constructor and assignment operator.

3. The whole idea of a class is that it represents a certain closed/separate entity. Your class does not follow this, it is merely a wrapper to call the multilication method. If you were following OOP principles, you would create a class matrix holding one of the matrices with an overloaded operator+(const matrix&), operator*(const matrix&) and so on, then create two matrix instances in main and call * on them. Like so:

   class Matrix {
   private:
       std::size_t ncol;
       std::site_t nrow;
       std::vector<double> storage;
   public:
       Matrix() : storage() {}
       Matrix(std::size_t ncol, std::size_t nrow) : ncol(ncol), nrow(nrow), storage(ncol*nrow) {}
   
       std::size_t getColCount() const { return ncol; }
       std::size_t getRowCount() const { return nrow; }
   
       double getAt(std::size_t col, std::size_t row) { [implement this] }
       Matrix operator+(const Matrix& other) const { [implement this] }
       Matrix operator*(const Matrix& other) const { [implement this] }
       [other operators...]
   
       static Matrix randomMatrix(std::size_t ncol, std::size_t nrow) { [ implement this ] }
   };
   
   int main() {
       Matrix matrixA = Matrix::randomMatrix();
       Matrix matrixB = Matrix::randomMatrix();
       Matrix matrixC = matrixA*matrixB;
   }
   

Answer 2

Problems

function shadowing

Your code doesn't compile here:

142807.cpp: In member function ‘void matrix::free()’:
142807.cpp:48:23: error: no matching function for call to ‘matrix::free(double (*&)[40000])’
          free(matrix_a);
                       ^
142807.cpp:46:12: note: candidate: void matrix::free()
       void free()
            ^~~~
142807.cpp:46:12: note:   candidate expects 0 arguments, 1 provided

The reason being that your class's free() hides the ::free defined in <stdlib.h>. You could re-write that as

  void free()
  {
     ::free(matrix_a);
     ::free(matrix_b);
     ::free(matrix_c);
  }

but, as mentioned in other answer, the correct tool is the destructor.

---

includes

You can remove most of the header includes - your code uses only <stdlib.h> and none of the others.

---

fixed constants

It's not clear why you've chosen the magic number 40000 or the type double. This might be the moment to introduce templates; see my code below.

---

other issues

Additionally, all the points in Eichhörnchen's answer are all useful; I won't repeat them here.

---

Example code

I ended up with a template class, and chose to use C++17 Concepts Lite as implemented in g++. I also chose to compile with OpenMP to speed up the computations.

For the element type, almost anything that can implicitly convert from the integer 1 is acceptable - char, double, std::complex<float>, etc.

#include <array>
#include <memory>

template<typename T, size_t X, size_t Y>
    // ensure that we can create zero and unity elements
    requires(T{} != T(1))
class matrix
{
    static constexpr size_t WIDTH = X;
    static constexpr size_t HEIGHT = Y;

instance data

The elements are stored in an array. We want to allocate it from heap, as it's likely too big for stack, so we use a smart pointer to release it when the object is destroyed.

    using array_type = std::array<T,X*Y>;
    std::unique_ptr<array_type> p = std::make_unique<array_type>();

A couple of convenience methods to access specific elements (you might choose to make value() public if you want):

    static constexpr size_t index(size_t x, size_t y) { return y*X + x; }
    T& value(size_t x, size_t y) { return (*p)[index(x,y)]; };
    const T& value(size_t x, size_t y) const { return (*p)[index(x,y)]; };

constructors, destructor and assignment operators

public:

    // default constructor - creates a zero matrix (identity for addition)
    matrix()
    {
        p->fill({});
    }
    // copy constructor
    matrix(const matrix<T,X,Y>& other)
    {
        *p = *other.p;
    }
    // move constructor
    matrix(matrix<T,X,Y>&& other)
        : p(std::move(other.p))
    {
    }

    ~matrix() = default;

    // copy assignment
    matrix& operator=(const matrix& other)
    {
        *p = *other.p;
        return *this;
    }
    // move assignment
    matrix& operator=(matrix&& other)
    {
        std::swap(p, other.p);
        return *this;
    }

factory methods

A utility to produce an identity matrix; you might think of others to add to this:

    static matrix identity()
        requires(X == Y)
    {
        matrix m;
        for (size_t i = 0;  i < X;  ++i)
            m.value(i,i) = 1;
        return m;
    }

equality and inequality operators

    bool operator==(const matrix& other) const
    {
        return *p == *other.p;
    }

    bool operator!=(const matrix& other) const
    {
        return *p != *other.p;
    }

simple operators

Element-wise multiplication (division is left as an exercise for the reader):

    matrix& operator*=(T t) {
        for (auto& v: *p)
            v *= t;
        return *this;
    }

You'll probably want an element-wise addition operator for a matrix of the same element type and dimensions; I'll leave that as an exercise, too.

matrix multiplication

We use the template arguments to verify at compile time that the dimensions and element types are compatible:

    template<typename T2, size_t Z>
    auto operator*(const matrix<T2,Y,Z>& other) const
    {
        matrix<decltype(T{}*T2{}),X,Z> product;
#pragma omp parallel
        for (size_t i = 0;  i < product.WIDTH;  ++i) {
            for (size_t j = 0;  j < product.HEIGHT;  ++j) {
                auto& val = product.value(i,j);
                for (size_t k = 0; k < WIDTH; k++)
                    val += value(i,k) * other.value(k,j);
            }
        }
        return product;
    }
};

non-member operators

Here are some operators that don't need to be members:

template<typename T1, typename T2, size_t X, size_t Y>
auto operator*(const matrix<T1,X,Y>& m, T2 t)
{
    auto result = m;
    return result *= t;
}
template<typename T1, typename T2, size_t X, size_t Y>
auto operator*(T2 t, const matrix<T1,X,Y>& m)
{
    return m * t;
}

main()

Finally, let's exercise it, by testing that identity multiplication and value multiplication are both commutative:

int main()
{
    static constexpr int SIZE = 4000;
    auto a = matrix<double, SIZE, SIZE>();
    auto b = matrix<double, SIZE, SIZE>::identity();
    return b*a != a*b
        || 2*b != b*2;
}