Solving shifted linear system

[mstewart]

@stevengj

Update: checking the source code it seems that it does form an ldlt factorization explicitly, however, though it is still true that it is O(n); in principle this could probably be elided into the solve.

That’s not ideal in terms of stability. It looks like that is currently using LDL^T with no pivoting and no blocking so that D is really diagonal without any 2\times 2 blocks. That can be unstable:

A = SymTridiagonal([
  3e-10 1 0
  1     2 1 
  0     1 2
])

F = ldlt(A)
display(A - F.L * F.D * F.L')

gives

3×3 Matrix{Float64}:
 0.0  0.0         0.0
 0.0  4.76837e-7  0.0
 0.0  0.0         0.0

There’s a blocking strategy due to Bunch that doesn’t do any interchanges, but does use 2\times 2 block pivots and was proven to be stable by N. Higham in this paper. It’s still linear complexity. Edit: L and D do require more storage in the stable algorithm, so ldlt! isn’t something that can be done fully in-place, but it seems worth it.