I wrote a custom function computing the rank of some matrices (I need an exact computation, and the entries are mod p integers for some small prime p which is why I’m not using the built-in rank function). This is nothing fancy, just a fairly straightforward Gaussian elimination, the code is below. I’ve been using it on larger and larger matrices and it’s working pretty well and giving me an answer after a couple minutes.
Yet, when the matrix is big enough (roughly 40.000x70.000), the computation starts, do something like 70% of the job fairly quickly and then just stop working, ie it’s not crashing or throwing any error, it’s just not doing anything at all.
Weirdly enough the following ugly hack works: I modified my function in order to periodically save the current state of the matrix on the disk. Every time this happens, I kill Julia, start again, load the matrix I just save and start computing the rank again. If I do this one or two times I get my result again in a couple of minutes.
It kinda seems like this is a memory issue, although my matrices are not that big… Does anyone knowns what’s going on ?
Here’s the fairly naive code for the rank. A sparse matrix is represented as a Vector of SparseVector’s. The first function swaps the first column with a short(er) one. The second is just basic Gaussian elimination, deleting the pivot after it’s been used to save some memory.
function sparseFindShort!(n,ll) i=0 j=length(ll)+1 @inbounds for k in 1:n m=nnz(ll[k]) if m==1 i=k break end if m<j && m>0 i=k j=nnz(ll[k]) end end if i!=0 ll,ll[i]=ll[i],ll end return nothing end function sparseRank!(n,l) if n>=1 sparseFindShort!(n,l) if nnz(l)==0 deleteat!(l,1) return sparseRank!(n-1,l) else ind=l.nzind l/=(l.nzval) @inbounds Threads.@threads for j in 2:n if l[j][ind]!=0 l[j]-= l[j][ind]*l end end deleteat!(l,1) return 1+sparseRank!(n-1,l) end else return 0 end end