I'm trying to improve efficiency of my code in which main task is a solution to huge (~$10^4\times 10^4$) but sparse linear system
(In fact my aim is to solve nonlinear equations $F(x)=0$, but the problem can be reduced with Newton's method to iterative procedure, with $A=J=F'(x)$ as Jacobian matrix of a system, for which I have an explicit form. Since I have good guess for a starting point it takes only 1-2 steps for
FindRoot to converge, but it consumes great ammount of RAM ~9GB, which increases with number of steps taken). The problem is both the time it takes for Mathematica to get the solution and the ammount of memory it uses e.g.
LinearSolve[A, b, Method-> "Multifrontal"]
takes around 140s (from the beginning it works on two cores) while the memory consumption is enormous, with a peak of 6GB (where it comes from?). Thanks to the question (82645) I've found the
"Pardiso" method which is better in terms of memory (at the beginning it consumes only 1GB, at the end ~3.1GB) but is much slower 415s. Interesting thing about
"Pardiso" is that for most of the time (6*60s) it uses only one core, then for the next 60s it uses two (see CPU load graph below; 2 physical cores with HT) and consumes additional 2.1GB of RAM. What it does at these phases (unpacking sparse array, data transfer to MKLs Pardiso solver)?
Is there any chance to improve efficiency for this problem? Both in terms of execution time and memory consumption? Matrix $A$ is quite sparse (density $\leq 1\%$), at the end of the day I would like to increase size of the the problem (dimension of $A$) wich with current timings will make it impossible.
I think that this is independent problem of any details of the system (its sparsity pattern). Below I've uploaded files with compressed (use
Uncompress) $A$ matrix (large file 80MB) and $b$ vector (the $32769\times 32769$ test case).
I'm suspecting that
"Multifrontal" method internally in this case (consistent timings, but for 2 intermediate steps I see that it consumes ~1.5 times memory as
LinearSolve). Iterative methods are not an option here (as far as I know) since these take forever (both for
LinearSolve[..., Method -> "Krylov"] as these seem to be equivalent).