   Sparse LU decomposition
  I'm implementing a linear algebra package for my game, (and not using an existing package just because), and I have LU decomposition and I want to build a sparse implementation.  Unfortunately I can't seem to find any (free) articles which talk about it.I understand that I need to store my matrix in a format to take advantage of its sparsity.  I can think of two methods: first, store a value in the matrix and its corresponding row and column.  Second, store the values as a dense array, and then store a bitfield to record where the 0s/non 0s are.  The first potentially wastes a lot of space on indices, but it scales linearly with the number of non zero elements.  But it doesn't support random access very well.  The second is space efficient but it still scales with the square of the matrix size.  And it also doesn't support random access very well.  There are probably a few other methods that I'm not even thinking of.  Plus, of course, I'm entirely ignoring cache thrashing issues and scalability to multiple cores (ala z ordering on the matrix).The second problem is performing the decomposition on the sparse matrix in a way to preserve sparseness.  The Numerical Recipes briefly mentions one method: basically choosing as pivot that row with the fewest number of non zero elements (IIRC, I couldn't find the spot where it talks about it).  Plus, of course, the LU decomposition method ties in with the method of sparse storage.  If you decompose in a way that creates more non zero elements but allows for more efficient forward and backward substitution in the future it's a net win.And last, I want to understand the theoretical order of complexity of the decomposition and forward and backward substitutions on the sparse matrix, for both worst case and average case.Obviously these are pretty complex topics, so I'm not looking for a rote solution.  I'm just not finding much relevant literature online or in the few text books I have on hand.  Mostly just a sentence or maybe a section here and there.  And never a comparison of different methods.   
  You can take a look at boost::ublas.  Even if you don't want their code, there's documentation for it.  I think they explain their storage options in the docs.   ;  UMFpack is used in matlab, and the source code is available.  I downloaded the paper free a while ago, but don't know where it is.As for storage methods, the one you didn't mention is operator form.  I store some sparse matrices in operator form allowing me to do A*x fast.  For FEM, you can store the stiffness matrix as A = A1+A2+... and then do a multifrontal factorization without having to explicitly assemble (I think).  I thought of doing this but took a look at the papers and decided to be lazy and use an iterative method. Have you tried GMRES (out of curiosity)   ;  I took a look at the docs for boost::ublas and UMFpack but didn't glean anything useful.  I'd take a look at the source but generally I have a hard time reconstructing the math from the source (source tends to be written for speed instead of clarity of concept.  Plus people don't really play with proper notation in source (so few unicode .cpp files :P )).Quote:Have you tried GMRES (out of curiosity)Only that it exists and is sort of a hybrid direct/iterative method (iterative, but converges to an exact solution after a known number of iterations).I'm not sure if it applies to my problem, though.  Instead of Ax = b, I have something more like (A * B * C)x = b, where (A, B, C) are all the same size and B is a diagonal matrix and A and C are extremely sparse (and presumably not singular).QR decomposition keeps popping up when I search for articles, so there might be something there.   ;  The conjugate gradient method is a candidate for you type of system.  You can implement this for sparse matrices stored using a std::map.  A good resource for this and other matrix algorithms is "Matrix Computations" by Golub and Van Loan.;  Yea, CG is really great if your matrix is symmetric + positive definite.  GMRES is good for a general case, unless your matrix is ill-conditioned (restarted GMRES has a tendency to stall with ill-conditioned matrices).   ;  While the (A*B*C) matrix is symmetric PD, it isn't necessarily sparse.  Likewise the individual matrices are sparse, but not necessarily symmetric (not sure about PD for them).I'm aware of CG (that it exists.  I haven't looked into it), but I'm trying to avoid iterative solutions, at least for the moment.  Though I am aware that at a certain point numerical issues on a large system require an iterative refinement.  And of course an iterative method with a known number of iterations for convergence can be thought of as a direct method...There's really a whole wonderland of interesting linear algebra techniques, but I'm trying to understand sparse LU decomposition first before I move on, because I'm comfortable with the formulation for dense matrices and it isn't too hard to think about the crout/doolittle algorithm.   ;  You might look at the sparse cholesky factorizations b/c your matrix is positive definite.  These are much easier than sparse LU (and half the work).  Although since ABC is not necessarily sparse, I'm not sure you gain anything from a sparse formulation.   