Best way to solve specific block-tridiagonal linear system (10000x10000 and larger) To provide more context, this system came from energy balance equation on a mesh with (n,m) nodes in each direction. It's a linear system that looks like this (size of system in blocks n = 4, size of blocks m = 4, 'd' means 'some value'):
$$
\left| \begin{array}{cccc|cccc|cccc|cccc}\hline
1 &   &   &   & d & d & d & d &  &  &  &  &  &  &  & \\
1 &-1 &   &   &   &   &   &   &  &  &  &  &  &  &  & \\
1 &   &-1 &   &   &   &   &   &  &  &  &  &  &  &  & \\
1 &   &   &-1 &   &   &   &   &  &  &  &  &  &  &  & \\
\hline
1 &   &   &   & d & d &   & d & d &  &  &  &  &  &  & \\
  & 1 &   &   & d & d & d &   &  & d &  &  &  &  &  & \\
  &   & 1 &   &   & d & d & d &  &  & d &  &  &  &  & \\
  &   &   & 1 & d &   & d & d &  &  &  & d &  &  &  & \\
\hline
&  &  &  & 1 &   &   &   & d & d &   & d & d &  &  &  \\
&  &  &  &   & 1 &   &   & d & d & d &   &  & d &  &  \\
&  &  &  &   &   & 1 &   &   & d & d & d &  &  & d &  \\
&  &  &  &   &   &   & 1 & d &   & d & d &  &  &  & d \\
\hline
&  &  &  &   &  &  &  & 1 &   &   &   & d & d &   & d \\
&  &  &  &   &  &  &  &   & 1 &   &   & d & d & d &   \\
&  &  &  &   &  &  &  &   &   & 1 &   &   & d & d & d \\
&  &  &  &   &  &  &  &   &   &   & 1 & d &   & d & d \\
\hline
\end{array} \right|
$$
[Part of Update 3: further down you may find a genius idea of an uninformed lunatic. The idea below ruins the diagonal dominance for the matrix, so while this New Revolutionary Elimination is significantly faster, it's useless.]
By moving first m equations to the last line, i get almost-upper-triangular system:
$$
\left| \begin{array}{cccc|cccc|cccc|cccc}\hline
1 &   &   &   & d & d &   & d & d &  &  &  &  &  &  & \\
  & 1 &   &   & d & d & d &   &  & d &  &  &  &  &  & \\
  &   & 1 &   &   & d & d & d &  &  & d &  &  &  &  & \\
  &   &   & 1 & d &   & d & d &  &  &  & d &  &  &  & \\
\hline
&  &  &  & 1 &   &   &   & d & d &   & d & d &  &  &  \\
&  &  &  &   & 1 &   &   & d & d & d &   &  & d &  &  \\
&  &  &  &   &   & 1 &   &   & d & d & d &  &  & d &  \\
&  &  &  &   &   &   & 1 & d &   & d & d &  &  &  & d \\
\hline
&  &  &  &   &  &  &  & 1 &   &   &   & d & d &   & d \\
&  &  &  &   &  &  &  &   & 1 &   &   & d & d & d &   \\
&  &  &  &   &  &  &  &   &   & 1 &   &   & d & d & d \\
&  &  &  &   &  &  &  &   &   &   & 1 & d &   & d & d \\
\hline
1 &   &   &   & d & d & d & d &  &  &  &  &  &  &  & \\
1 &-1 &   &   &   &   &   &   &  &  &  &  &  &  &  & \\
1 &   &-1 &   &   &   &   &   &  &  &  &  &  &  &  & \\
1 &   &   &-1 &   &   &   &   &  &  &  &  &  &  &  & \\
\hline
\end{array} \right|
$$
Using the gaussian eliminiation, i can upper-triangularize it in O(nm^2 + m^3) operations.
Other approach i've tested is to use matrix version of Thomas algorithm, but it costs O(nm^3) actions.
I really need this speed-up of gaussian elimination, as both n and m are going to be somewhat like 100-200-300, but gaussian elimination aggregates double rounding error tremendously fast. With N and M over 15, roots of such system calculated with the use of elimination already differ from "true" roots (calculated with the use of LU decomposition) by relative difference of 0.1. If i scale system to even larger M,N, the only thing that i get right with elimination is an order of magnitude of roots, and even that i wont get anymore around N=M=70.
Thomas algorithm works perfect, but i need this speed-up badly.
So, the question is. Is there a way to solve such system in less time than O(n*m^3) without losing precision? Mathematically such thing is possible (gaussian elimination ftw), but is there a way i can get precise for at least 3-4-5 significant values roots fast? 
Update 1:
forgot to post some more info:
precision fall comes from pure gaussian elimination of last lower-right m-sized block. I could solve that last block as separate system with some precise algorithm, but that would complicate overall method a lot. Is there any other ways?
Update 2:
Sadly, information in Update 1 is incomplete. Further test have shown that, while last block indeed causes most of error, the whole elimination process wrecks precision. If i do not triangularize last block and try to solve the entire system from this point with QR/LU, around N=M=50 i only get order of magnitude of roots correctly.

Update 3:
Huh, wasn't expecting any answers after such a pause :D
But @Armadillo Jim noticing this post means it still comes up on some searches, so i guess i can provide an update for everyone interested. Answering his question on why i don't use iterative solvers:
The matrices are very badly conditioned.
From time to time i have to solve batches of 100+ such systems corresponding to radiative transfer in different parts of spectrum, and noble gases' spectra are not the smoothest things that exist. This results in a ton of systems that only have their shape similar.
That means if i am to use iterative solver, i have to recalculate preconditioner for every matrix each time i want to solve aforementioned batch of systems, and preconditioners are just wasted after that, as i am solving those systems exactly once for any given left part (A in AX=B). Of course i could use no preconditioner at all, but this results in thousands of iterations and outright diverging.
If somebody is interested, possible ways of solving those systems on regular meshes efficiently are the family of O(N^1.5)..(N^2) Nested Dissection methods (Nested Dissection of a Regular Finite Element Mesh, Alan George, Siam J. on Numerical Analysis, vol. 10, #2, Apr. 1973, 345-363) and their somewhat descendant based on specific properties of hierarchicaly-blockseparable matrices (O(N), A.Gillman, Fast direct solvers for elliptic partial differential equations, phd thesis, 2011).
Sadly i may have to move to adaptive meshes and cell-based approach, as now i have to add gas dynamics to the model and single test may run for a week straight on our department's cluster, and if that's the case, i can as well raise mesh resolution and use iterative solvers or what else is used for sparse matrices.
 A: Ah, this question seems to be quite old so you probably found a solution or went on another problem by now, but I just could not resist making some observations.
My first reaction as it is a sparse system would be to use some Krylov subspace method as they use matrix-vector multiplications, and then factor out $d$ as it is recycled many times over should save many scalar multiplications you will then be left with a (larger) sparse matrix full of $\pm 1$ elements so in that pass you can replace multiplication in the scalar products with addition and subtraction altogether.
Also you have "leave one out" structure on some of the blocks, they would benefit computationally from instead factoring out the sum and take difference between sum and individual elements.
$$\left[\begin{array}{rrrr}1&1&0&1\\1&1&1&0\\0&1&1&1\\1&0&1&1\end{array}\right] = \left[\begin{array}{rrrrr}1&0&0&-1&0\\1&0&0&0&-1\\1&-1&0&0&0\\1&0&-1&0&0\end{array}\right] \left[\begin{array}{cccc}1&1&1&1\\1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right]$$
3 additions and 4 subtractions instead of 8 additions.. hmm not as impressive for 4x4 blocks, but it will get much more efficient on larger blocks. For example 8x8 you would have 8 rows with 6 additions each for 40 additions in total but with factoring you would have 7 additions and 8 subtractions.
In general we would go from $n(n-2) \rightarrow 2n-1$. From quadratic to linear complexity.
