Created
May 20, 2022 15:15
-
-
Save YingboMa/b77de655fa60f4ea132236e16b94c516 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
=# |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment