periodic_qr.jl

using LinearAlgebra, BandedMatrices

n = 3000
A = Matrix(SymTridiagonal(Symmetric(rand(n, n)))); A[end, 1] = A[1, end] = 1;
b = rand(size(A, 1))

function periodic_givens_rotation!(A::AbstractMatrix{T}) where T
    m, n = size(A)
    gs = Vector{LinearAlgebra.Givens{T}}(undef, 2n-2)
    jj = 0
    @inbounds for j in 1:n-1
        g1, r = LinearAlgebra.givens(@view(A[:, j]), j, m)
        gs[jj+=1] = g1
        #A = g1 * A
        A[j, j] = r
        A[m, j] = 0
        lmul!(g1, @view(A[:, j+1]))
        if n-1 != j && n-1 != j+1
            lmul!(g1, @view(A[:, n-1]))
        end
        if n != j && n != j+1
            lmul!(g1, @view(A[:, n]))
        end

        g2, r = LinearAlgebra.givens(@view(A[:, j]), j, j+1)
        gs[jj+=1] = g2
        A[j, j] = r
        A[j+1, j] = 0
        lmul!(g2, @view(A[:, j+1]))
        j + 2 <= n && lmul!(g2, @view(A[:, j+2]))
        #A = g2 * A
        if !(j <= n-1 <= j+2)
            lmul!(g2, @view(A[:, n-1]))
        end
        if !(j <= n <= j+2)
            lmul!(g2, @view(A[:, n]))
        end
    end
    gs, UpperTriangular(A)
end

function solveAb!(A, b)
    gs, R = periodic_givens_rotation!(A)

    bb = copy(b)
    for g in gs
        lmul!(g, bb)
    end

    R11 = BandedMatrix(@view(R[1:end-2, 1:end-2]), (0, 2))
    R12 = @view R[1:end-2, end-1:end]
    R22 = @view(R[end-1:end, end-1:end])
    b1, b2 = @view(bb[1:end-2]), @view(bb[end-1:end])
    x2 = R22 \ b2
    x1 = R11 \ (b1 - R12 * x2)
    [x1; x2]
end
norm(A * solveAb!(copy(A), b) - b)

using BenchmarkTools, SparseArrays
@btime (lu!(B) \ $b) setup = (B=copy(A)) evals=1;
@btime ldlt($(sparse(A))) \ $b;
@btime solveAb!(B, $b) setup = (B=copy(A)) evals=1;
#=
julia> @btime (lu!(B) \ $b) setup = (B=copy(A)) evals=1;
  108.699 ms (4 allocations: 46.97 KiB)

julia> @btime ldlt($(sparse(A))) \ $b;
  653.285 μs (69 allocations: 1.31 MiB)

julia> @btime solveAb!(B, $b) setup = (B=copy(A)) evals=1;
  448.327 μs (30 allocations: 516.44 KiB)
=#