I need to calculate $\log(\det (\mathbf M_i))$ where the $\mathbf M_i$'s are large sparse matrices, which are real, symmetric and positive semi-definite. I hope to have between $10$ and $100$ of those. Physically, each $\mathbf M_i$ is the stiffness matrix of a large system, and if it helps, I know how to write it as $\mathbf M_i=\mathbf B_i^T\mathbf B_i$ (where $\mathbf B_i$ is not a square matrix).
I am not a computational scientist, so I would prefer using written packages rather than implementing my own algorithms. I currently work in MATLAB and/or Mathematica, but I am also fluent in C++, and am willing to work in any other enviroment if this could help.
What I've tried so far is calculating the eigenvalues (using MATLAB's
eig) and summing up their logs, but I suspect that there are other, more efficient ways to do so. For example, I stumbled upon this paper, but I am not sure whether its algorithm is suitable in my case, and I could not find an implemented code.
Another complication is that it might happen that $\mathbf M_i$ is not invertible, i.e. $\det\mathbf M_i=0$, and in this case I'm interested in the product of all the non-zero eigenvalues (in these cases I know exactly how many zero eigenvalues are there).
Any help would be much appreciated.