The fact is, a 3kx3k system is not that slow to solve. So unless you are doing this solve many thousands of times, keeping it sparse may well slow you down.
It is not even really that sparse, since for a 3100x3100 matrix with 250000 nonzeros, there are on average
So roughly 80 non-zeros per row. That is not very sparse.
Have you tried a direct solve in sparse form, using a column ordering to reduce fill-in?
Have you compared the solve time if the matrix were a full matrix?
That test was run on my computer. I was watching the actuivity monitor on my computer as it ran, and the test itself was so fast that MATLAB never even fired up multiple cores in the solve, at least that the activity monitor showed in the time slice it was watching.
Even for for a full matrix of that size, it only took a small fraction of a second to solve as a full matrix.
Are you truly needing to solve this probem many thousands of times? If not, then this entire question seems moot.
If so, then is your matrix fixed, so just different right hand sides? If so, then a matrix factorization would be right, or if you have many right hands sides, known in advance, then a direct call to backslash is absolutely correct.
Or is your real problem orders of magnitude larger?
The point of my answer is this really is a small probem, and sparse form may not even be a good choice. The amount of time tken by the incomplete LU there is often a significant fraction of the time it would take to perform the matrix factorization itself.
So if you want help, then it does help if you explain all of these things, because your question as it is does not completely make sense. Why do you think you need to solve this problem as you wish to solve it?
I see that you have provided the matrix L. Must I point out that L is singular? You claim that \ is slow here. The probem is, that backslash, when applied to a FULL matrix of that size, is actually roughly 6 times faster for the full matrix than it is for the sparse matrix.
To remove the issue of your matrix also being singular, I'll perturb the diagonal by a small amount.
Lhat = L + speye(3139)*1e-7;
So if you just store the matrix in a full form, \ is significantly faster for the full matrix! More than 6x faster.