This question borders between a programming and math question (more math). I have a little matrix knowledge but this is past my ability, so any help is very much appreciated.


Question
--------


I have a very large sparse matrix (say one million rows square). I'd like to solve this for many different sets of constant values (to clarify, by constant values, I mean the 1xn matrix 'b' in Ax = b). Performance is absolutely crucial - so it would make sense to do whatever simplification I can in advance of introducing each set of constants.


Goal
----


I'd like to get the time complexity for solving each constant set down to O(n)... I may be dreaming however.


Other Info
----------


* The matrix really is a million (maybe even 100 million) rows
* It's *very* sparse
* It has a bandwidth ranging between 2 and as much as a few hundred
* Even for larger bandwidths, each row has max. 7 or 8 non zero entries
* It takes an (extremely rough) diagonal form
* Each matrix 'b' in Ax = b is actually very sparse itself
* Due to the nature of the problem, a solution should always be available
* All nonzero values are real, floating point numbers. They can be negative however.



I have a very large sparse matrix, 'A', in Ax = b. What work in advance of getting 'b' can be done to reduce solving time?
Instead of reinventing the wheel, look at the major linear algebra packages. Many have procedures for large sparse matrices. One method is "sparse $LU$ decomposition" (google that phrase for tutorials on it). Once the matrix is factored, solving multiple equations with different right sides becomes fast. The difficult part is factoring into sparse factors, which is tricky, and that's why it would be best to use existing code rather than write your own. It is theoretically possible for a sparse matrix to be impossible to factor into sparse factors (the additional non-zero entries are called "fill-in") but there are heuristics that attempt to keep the factors sparse and in practice they are often successful. The only sure way to know if it will work is to try it.


Note that if your matrix is symmetric there will be special versions for that case.


