Are you sure you need to do this?
Be aware that determinants for any reasonable sized matrix are very expensive to calculate, so if you're actually doing something else (linear solve, eigen-solve etc) there is probably a better way without calculating the determinant. And there are many external libraries which will solve this stuff for you.
In a real code review, I would be asking why you are implementing this yourself.
Getting an idea of the complexity
The best known algorithms are approx O(n2.4), and the best actually implemented algorithms are O(n3) (as discussed in this SO question). That would put a rank 50 matrix determinant about 4600x slower than a 3x3. So if you are going to need determinants of large matrices, make sure your method will permit that to calculate in an acceptable time frame.
This method, if I understand it correctly, calculates the determinants of n minor matrices each of rank n-1. And that happens recursively. So y(n)=n.y(n-1) -> y(n)/y(n-1) = n. This is a method that will be O(n!)... at 50x50 it will take of order 1063 times longer than the 3x3.
Critique of the code itself
Aside from the questions over the method, we can look at the code itself.
Data structure
An array of pointers to arrays of doubles has some issues in terms of performance; every access of a matrix element requires 2 dereferences; one to the array of pointers, and then from there to an array of doubles.
Although it looks the same, this is less efficient than a multidimensional C array like double[4][4] - being a fixed size makes it possible for the lookup to index directly into the correct member without the additional dereference step; element [2][3] is always the same fixed distance from the beginning, and the whole array can be allocated as a single contiguous block.
The equivalent thing with vectors would not be a vector of vectors (which again requires 2 dereferences) but a single vector of size n2 which you index into via a function. That could be put together as a matrix class:
class SquareMatrix
{
vector<double> _elements;
size_t _rank;
public:
SquareMatrix(size_t rank) : _rank(rank), _elements(rank*rank)
{
// vector empty(rank*rank);
// swap(_elements, empty);
}
double & elem(size_t i, size_t j)
{
return (_elements.begin() + i_rank + j);
}
};
Usage:
SquareMatrix a(10);
a.elem(2,3) = 5;
std::cout << a.elem(2,3);
The improvement from this sort of change should be quite significant, partly from the reduction in dereferencing, but most significantly because the CPU will be able to usually cache the right values - much of the time the next value is adjacent to the prior one. This kind of predictable memory usage is also friendly to branch prediction.
Copying copying
The code is creating a copy of each minor in order to calculate its determinant; so for each minor you copy (n-1)2 doubles in. That small part itself could be improved by reusing the matrix; instead of copying columns b,c,d,e in, then copying a,c,d,e in, you could just replace the relevant values:
a b c d f g h a c d
e f g h -> j k l then i k* l*
i j k l n o p m o* p*
m n o p
See that the marked elements have not changed. So at the very least this part could change to only copy in 2n-1 elements rather than (n-1)2. At 50x50 that means 96% fewer copies.
Small things
I would prefer pre-increment (++value
) over post-increment (value++
) because it is always unambiguous and thus a preferable default even when either would work.
You can do sign = -sign
instead of sign = -1 * sign
; I don't know whether the compiler will optimise it away, but a negation is much faster than a multiplication.
Memory management in C++ is not like in Fortran
Your code leaks almost all the memory it allocates. For a 50x50 matrix it would leak 50 matrices of size 49x49, and recurse to 50 rank-49 determinants which each leak 49 matrices of size 48x48, which is something like O((n+1)!) leaking.
calloc
is a C thing; the C++ equivalent is to call new
. Either way, if you calloc
some memory you must call free
on it to release it, and if you new
something you must delete
it.
Generally in modern C++ you avoid calling either, and instead rely on structures that manage the memory for you, like vector<double>
- this will free the memory it allocates when the vector goes out of scope, which is probably what you're expecting if you're used to Fortran 95/2003.
If you still want to allocate memory via new
, then you can use a smart pointer to ensure that the memory is freed when the array goes out of scope. In C++11 that would likely be a unique_ptr or a shared_ptr. Note that you couldn't just wrap your array of pointers in the smart pointer, or it would only free that array and not all the arrays those pointers point to. It would need to be an array of smart pointers, so that each deletes its array when the top level array deconstructor runs.
Summary
Whilst this is a good initial attempt, and not a bad coding exercise, this code will not achieve what you want to achieve. Fixing the leaks will be an educational exercise, but will not result in code that can handle the sizes of matrices that you state you need to cope with. In general anything that needs to cope with matrices of any significant size requires a different kind of method - generally matrix methods that can be done by hand do not scale well.