System of linear equations, resulting from a weighted graph. How to solve this numerically? We have a problem that leads to a system of linear equations which has to be solved numerically. There are thousands of algorithms to solve linear equations, but I haven't found any that fits our special requirements.
I've tried to describe these requirements in the pdf linked below:
http://pdfcast.org/pdf/linalgproblem (Document has three pages)
Since our problem is probably a common one, we believe that there are algorithms specially designed for the kind of matrices we are facing.
Do you know someone who had the same problem?
Do you know a fitting algorithm? Do you know an implementation in c++? That would help us a lot!
PS: The matrices to solve will have about 1000-10000 rows and 5000-20000 lines.
EDIT: I still hope to hear about a more specific recommendation than to use an algorithm for sparse matrices. I believe I'm not the first one having to solve this problem.
 A: OK, I read your PDF.  I'd like to offer some comments which won't fit well into the comment box, hence this pseudo-answer:


*

*As you note, your system is likely to be (at least partially) overdetermined, and the data is likely to have errors.  This makes a big difference in the kind of algorithm you're going to need.

*You mention that consecutive routes are likely to differ by one node, e.g. A, AB, ABC, ABCD, etc.  Are the measurement errors for such routes likely to be correlated?  If so, you'll also need to take this into account.


For the first point, you might try formulating your problem as something like
$$\begin{aligned}
d(A,B) + d(B,E) &= 3.5 + e_1 \\
d(A,B) + d(B,C) + d(C,N) &= 8.6 + e_2
\end{aligned}$$
etc., and try to minimize the sum of the squares of the error terms $e_i$ (possibly with weighting factors to account for varying measurement accuracy) subject to the given constraints.  I'm not familiar enough with these kinds of linear fitting problems to recommend a specific algorithm or implementation, but I believe such do exists.
As for the correlations, in the simple case where $E[e_{i+1}|e_1,\dotsc,e_i] = e_i$ within a given chain of routes, (i.e. where errors within a chain are simply cumulative), there is fortunately an easy way to decorrelate the data: just subtract the previous route from the next within each chain.  This will also conveniently simplify your data, giving you lots of independent single-edge measurements instead of lots of overlapping multi-edge measurements.
A: I think I have solved the problem. (Though, we haven't tested it ye.)
The answer is Singular Value Decomposition. It gives us the possibility to discover column degenracies in our equation. We can remove the unsolvable variables from the equation and then calculate the least square solution for the rest of the unknowns.
When I look at the algorithm, I feel pretty sure it will also be possible optimize it for integer values. 
I've found all this in the book "Numerical Recipies, Second Edition, Cambridge University Press."
