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;
}
