# Efficient solution of huge sparse linear system

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 $$Ax=b$$ (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 FindRoot uses "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 SparseArrayKrylovLinearSolve[] and LinearSolve[..., Method -> "Krylov"] as these seem to be equivalent).

• Is A constructed as an explicit SparseArray? If not, try making it one to see if that improves on memory use and/or LinearSolve speed. – Daniel Lichtblau Nov 15 '15 at 22:36
• Have you seen this? – J. M. is away Nov 16 '15 at 3:11
• @DanielLichtblau Yes, $A$ is constructed as SparseArray object. @J.M. Yes I've seen that, I'm calling FindRoot with specified Jacobian function FindRoot[F[X], {X, X0}, Method -> {"Newton", "StepControl" -> "LineSearch", "UpdateJacobian" -> 1}, Jacobian :> J[X], AccuracyGoal -> 10, PrecisionGoal -> 10] with F and J a 'black box' functions. – mmal Nov 16 '15 at 9:02
• I suppose that one solution would be to explicitly call Pardiso solver from Fortran and use LibraryLink`. I'll try to experiment with that later today. – mmal Nov 16 '15 at 10:46