How Do I solve large systems given UMFPACK memory limitations?

Question

I am trying to solve a system of equations (A x = b) for 3D heat diffusion (i.e. each equation has at most 7 terms not including the constant "b" term) using UMFPACK with boost numeric bindings to C++.

e.g.:

namespace umf = boost::numeric::bindings::umfpack;

umf::umf_solve(A,x,b);

I can populate the A matrix and the "b" vector for a 200K+ set of equations while using roughly 200MB of memory. Then while in the UMFPACK routines, the memory used by my program jumps up to 2.3GB+.

Furthermore, at some point above ~211K equations, the routine passes back all zeros for "x", without throwing an error. I assume this is related to some kind of memory limit.

This is the same on every machine I've tested on (my MacBook, my Ubuntu VM, and a node on a Red Hat supercomputer).

My questions are:

1. Why is UMFPACK returning zeros without an error?
2. Why does UMFPACK require so much memory?
3. What options do I have for solving the system for larger sets of equations?

Answer 1

I can't give you a good answer for 1, but I can give you decent answers for 2 and 3.

1. I'm not terribly familiar with the Boost interface to UMFPACK. In the C interface, normally, you call UMFPACK routines that will allocate memory. If the memory cannot be allocated because there is not enough free memory available, UMFPACK will return a null pointer. You can then test for routines returning null pointers to handle errors; this approach is standard in C programming. Looking at the documentation, it seems like you'll have to check the return values of each function; it doesn't look like the library throws any exceptions. If you don't check return values, and only examine the solution returned by UMFPACK, that could explain the behavior you observe.

2. UMFPACK is a sparse direct solver; it carries out a sparse LU factorization of your $A$ matrix. As $A$ becomes larger, more of the entries of the LU factorization become structurally nonzero -- they are "filled in", and this phenomenon is called fill-in -- and structurally nonzero entries are represented explicitly in memory because they could be numerically nonzero. Fill-in is why you run out of memory; memory limitations tend to make sparse direct solvers unattractive for problems exceeding around 100K unknowns.

3. Use an iterative method for solving sparse linear equations instead of a sparse direct method. These methods don't necessarily require you to explicitly store $A$; they typically only require a method for calculating the matrix-vector product $Ax$ for any $x$. For best results, you should supply a preconditioner for the iterative method. For a diffusion problem, multigrid is a good preconditioner; multigrid methods are also worth using without a preconditioner for diffusion equations. Other people can probably offer better advice on iterative methods to use for this type of problem.

Answer 2

Geoff mentioned the limitations of sparse direct solvers. To solve a 3D problem of size $N = n^3$, an optimal direct solver will need $N^{4/3}$ memory and $N^2$ time. See this question for more details.

But it sounds like you have a structured grid diffusion problem. If the diffusion coefficient is smooth (or constant), you can solve the problem using about four floating point values worth of memory and 50 flops per grid point, using the celebrated Full Multigrid algorithm (developed by Achi Brandt in the 1970s). If the diffusion coefficient is very rough, it is easiest to use an algebraic multigrid method. Both geometric and algebraic multigrid methods are available in PETSc and other packages, depending on your needs. [Disclaimer: I am a developer of PETSc.]