# Solving shifted linear system

**URL:** <https://discourse.julialang.org/t/solving-shifted-linear-system/111265>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [March 6, 2024, 10:11pm UTC](https://discourse.julialang.org/t/solving-shifted-linear-system/111265 "2024-03-06T22:11:48Z")\
**Posts on this page:** 1\
**Showing post:** 10

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [March 10, 2024, 7:16pm UTC](https://discourse.julialang.org/t/solving-shifted-linear-system/111265/10 "2024-03-10T19:16:16Z")

</div>

> [@stevengj](#):
>
> _Update:_ [checking the source code](https://github.com/JuliaLang/julia/blob/5c7d24493ebab184d4517ca556314524f4fcb47f/stdlib/LinearAlgebra/src/tridiag.jl#L250) 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:

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

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

```

gives

```julia
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](https://www.sciencedirect.com/science/article/pii/S0024379598100745). 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.

---

_[View the full topic](https://discourse.julialang.org/t/solving-shifted-linear-system/111265)._