Skip to content

Fix whitespace #118

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Oct 14, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 2 additions & 3 deletions examples/blocktridiagonal.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
using InfiniteLinearAlgebra, BlockBandedMatrices, BandedMatrices, BlockArrays, InfiniteArrays, FillArrays, LazyArrays, Test

A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)),
Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)),
A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)),
Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)),
Vcat([fill(1.0,1,2),Matrix(1.0I,2,2)], Fill(Matrix(1.0I,2,2), ∞)))



2 changes: 1 addition & 1 deletion examples/intervalarithmetic.jl
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
using InfiniteLinearAlgebra
import IntervalArithmetic
import IntervalArithmetic
6 changes: 3 additions & 3 deletions examples/jacobi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ ql(Δ - λ*I).L[1,1]
end

λ2 = Interval(m,λ.hi)
ql(Δ - λ1.hi*I).L[1,1]
ql(Δ - λ1.hi*I).L[1,1]


ql(Δ - Interval(4.2)I).L[1,1]
Expand All @@ -41,7 +41,7 @@ Z,A,B = 0.5,-4.0,0.5

A = 4
rig_qltail(-2,10,B)

d,e = d2,e2

rig_qltail(Z,A,B,d,e)
Expand Down Expand Up @@ -106,4 +106,4 @@ X[2,:] .= (zero(T), X[1,1], X[1,2]);
X[1,:] .= (Z,A,B);
QL = _qlfactUnblocked!(X)

X, QL.τ[end]
X, QL.τ[end]
6 changes: 3 additions & 3 deletions examples/periodicschrodinger.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@ using InfiniteLinearAlgebra, Plots
# Consider a rank-1 perturbation of periodic semi-infinite Schrödinger (Jacobi) operator.
# We can construct it as a block matrix as follows:

A = BlockTridiagonal(Vcat([[0. 1.; 0. 0.]],Fill([0. 1.; 0. 0.], ∞)),
Vcat([[-1. 1.; 1. 1.]], Fill([-1. 1.; 1. 1.], ∞)),
A = BlockTridiagonal(Vcat([[0. 1.; 0. 0.]],Fill([0. 1.; 0. 0.], ∞)),
Vcat([[-1. 1.; 1. 1.]], Fill([-1. 1.; 1. 1.], ∞)),
Vcat([[0. 0.; 1. 0.]], Fill([0. 0.; 1. 0.], ∞)))

A[1,1] = 2 # perturbation
Expand All @@ -16,7 +16,7 @@ n = 100; scatter(eigvals(A[1:n,1:n]), zeros(n); label="finite-section")
# In any case, we can also calculate the spectrum using ∞-QL:
xx = [-0.95:0.05:0.95; 2.25:0.125/4:4.0]
plot!(xx, (x -> ql(A-x*I).L[1,1]).(xx); label="L[1,1]")
xx =
xx =
plot!(xx, (x -> ql(A-x*I).L[1,1]).(xx))


Expand Down
4 changes: 2 additions & 2 deletions examples/perttoeplitz.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,11 +23,11 @@ A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)), 0 => Vcat(a, Fill(0,∞)), 1 =>
x,y = range(-2,2; length=200),range(-2,2;length=200)
z = abs.(ℓ.(x' .+ y.*im))
contour(x,y,z; nlevels=50, title="z^2 + 1/(4z)")
scatter!(eigvals(Matrix(A[1:100,1:100]')); label="A' finite section")
scatter!(eigvals(Matrix(A[1:100,1:100]')); label="A' finite section")
scatter!(eigvals(Matrix(A[1:100,1:100])); label="A finite section")

x = range(-2,2,length=1_000)
plot(x, abs.(ℓ.(x.+eps()im)); label="abs(L[1,1])")
plot(x, abs.(ℓ.(x.+eps()im)); label="abs(L[1,1])")

a = [-0.1,0.2,0.3]
A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)), 0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
Expand Down
10 changes: 5 additions & 5 deletions examples/toeplitz.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,9 @@ using InfiniteLinearAlgebra, BandedMatrices, PyPlot
###

function ℓ11(A,λ; kwds...)
try
abs(ql(A-λ*I; kwds...).L[1,1])
catch DomainError
try
abs(ql(A-λ*I; kwds...).L[1,1])
catch DomainError
-1.0
end
end
Expand Down Expand Up @@ -64,7 +64,7 @@ clf(); qlplot(A; branch=findsecond, x=range(-2,3; length=100), y=range(-2.5,2.5;
clf(); qlplot(transpose(A); x=range(-2,3; length=100), y=range(-2.5,2.5;length=100)); symbolplot(A; color=:black); title("Limacon, transpose")


###
###
# bull-head
###

Expand Down Expand Up @@ -95,7 +95,7 @@ clf(); qlplot(transpose(A); branch=branch(2), x=range(-2,2; length=100), y=range
# Whale
###

A = BandedMatrix(-4 => Fill(im,∞), -3 => Fill(4,∞), -2 => Fill(3+im,∞), -1 => Fill(10,∞),
A = BandedMatrix(-4 => Fill(im,∞), -3 => Fill(4,∞), -2 => Fill(3+im,∞), -1 => Fill(10,∞),
1 => Fill(1,∞), 2 => Fill(im,∞), 3 => Fill(-(3+2im),∞), 4=>Fill(-1,∞))

clf(); qlplot(A; x=range(-15,20; length=100), y=range(-20,20;length=100)); symbolplot(A; color=:black); title("Whale")
Expand Down
4 changes: 2 additions & 2 deletions src/InfiniteLinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, si
show, getproperty, copy, copyto!, map, require_one_based_indexing, similar, inv
import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted

import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout,
import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout,
sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
_bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul, CNoPivot
import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
Expand Down Expand Up @@ -46,7 +46,7 @@ else
import Base: oneto, unitrange
end

if VERSION ≥ v"1.7-"
if VERSION ≥ v"1.7-"
LinearAlgebra._cut_B(x::AbstractVector, r::InfUnitRange) = x
LinearAlgebra._cut_B(X::AbstractMatrix, r::InfUnitRange) = X
end
Expand Down
6 changes: 3 additions & 3 deletions src/banded/hessenbergq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, F::QLHessenberg)
show(io, mime, F.L)
end

@inline function getL(F::QLHessenberg, _)
@inline function getL(F::QLHessenberg, _)
m, n = size(F)
tril!(getfield(F, :factors)[end-min(m,n)+1:end, 1:n], max(n-m,0))
end
Expand Down Expand Up @@ -78,7 +78,7 @@ for Typ in (:UpperHessenbergQ, :LowerHessenbergQ)
end
end

$Typ(q::AbstractVector{<:AbstractMatrix{T}}) where T =
$Typ(q::AbstractVector{<:AbstractMatrix{T}}) where T =
$Typ{T,typeof(q)}(q)
end
end
Expand Down Expand Up @@ -111,7 +111,7 @@ function materialize!(L::MatLmulVec{<:HessenbergQLayout{'L'}})
v = view(x, n:n+1)
mul!(t, Q.q[n], v)
v .= t
n > nz && norm(t) ≤ 10floatmin(real(T)) && return x
n > nz && norm(t) ≤ 10floatmin(real(T)) && return x
end
x
end
Expand Down
16 changes: 8 additions & 8 deletions src/banded/infbanded.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ function _inf_banded_sub_materialize(::BandedColumns, V)
b = parentindices(V)[1].b
data = bandeddata(A)
l,u = bandwidths(A)
if -l ≤ b ≤ u
if -l ≤ b ≤ u
data[u+1-b, max(1,b+1):end]
else
Zeros{eltype(V)}(∞) # Not type stable
Expand Down Expand Up @@ -86,7 +86,7 @@ function BandedMatrix{T}(kv::Tuple{Vararg{Pair{<:Integer,<:Vcat{<:Any,1,<:Tuple{
m,n = mn
@assert isinf(n)
l,u = lu
M = mapreduce(x -> length(x.second.args[1]) + max(0,x.first), max, kv) # number of data rows
M = mapreduce(x -> length(x.second.args[1]) + max(0,x.first), max, kv) # number of data rows
data = zeros(T, u+l+1, M)
t = zeros(T, u+l+1)
for (k,v) in kv
Expand All @@ -110,30 +110,30 @@ function BandedMatrix(Ac::Adjoint{T,<:InfToeplitz}) where T
A = parent(Ac)
l,u = bandwidths(A)
a = A.data.args[1]
_BandedMatrix(reverse(conj(a)) * Ones{T}(1,∞), ℵ₀, u, l)
_BandedMatrix(reverse(conj(a)) * Ones{T}(1,∞), ℵ₀, u, l)
end

function BandedMatrix(Ac::Transpose{T,<:InfToeplitz}) where T
A = parent(Ac)
l,u = bandwidths(A)
a = A.data.args[1]
_BandedMatrix(reverse(a) * Ones{T}(1,∞), ℵ₀, u, l)
_BandedMatrix(reverse(a) * Ones{T}(1,∞), ℵ₀, u, l)
end

function BandedMatrix(Ac::Adjoint{T,<:PertToeplitz}) where T
A = parent(Ac)
l,u = bandwidths(A)
a,b = A.data.args
Ac_fd = BandedMatrix(_BandedMatrix(Hcat(a, b[:,1:l+1]), size(a,2)+l, l, u)')
_BandedMatrix(Hcat(Ac_fd.data, reverse(conj(b.args[1])) * Ones{T}(1,∞)), ℵ₀, u, l)
_BandedMatrix(Hcat(Ac_fd.data, reverse(conj(b.args[1])) * Ones{T}(1,∞)), ℵ₀, u, l)
end

function BandedMatrix(Ac::Transpose{T,<:PertToeplitz}) where T
A = parent(Ac)
l,u = bandwidths(A)
a,b = A.data.args
Ac_fd = BandedMatrix(transpose(_BandedMatrix(Hcat(a, b[:,1:l+1]), size(a,2)+l, l, u)))
_BandedMatrix(Hcat(Ac_fd.data, reverse(b.args[1]) * Ones{T}(1,∞)), ℵ₀, u, l)
_BandedMatrix(Hcat(Ac_fd.data, reverse(b.args[1]) * Ones{T}(1,∞)), ℵ₀, u, l)
end


Expand Down Expand Up @@ -231,7 +231,7 @@ end
function BandedMatrix(A::PertToeplitz{T}, (l,u)::Tuple{Int,Int}) where T
@assert A.u == u # Not implemented
a, b = A.data.args
t = b.args[1] # topelitz part
t = b.args[1] # topelitz part
t_pad = vcat(t,Zeros(l-A.l))
data = Hcat([vcat(a,Zeros{T}(l-A.l,size(a,2))) repeat(t_pad,1,l)], t_pad * Ones{T}(1,∞))
_BandedMatrix(data, ℵ₀, l, u)
Expand Down Expand Up @@ -509,7 +509,7 @@ copyto!(dest::AbstractArray, L::Ldiv{BidiagonalToeplitzLayout,Lay}) where Lay =

# copy for AdjOrTrans
copy(A::Adjoint{T,<:BandedMatrix{T,<:Any,OneToInf{Int}}}) where T = copy(parent(A))'
copy(A::Transpose{T,<:BandedMatrix{T,<:Any,OneToInf{Int}}}) where T = transpose(copy(parent(A)))
copy(A::Transpose{T,<:BandedMatrix{T,<:Any,OneToInf{Int}}}) where T = transpose(copy(parent(A)))


##
Expand Down
2 changes: 1 addition & 1 deletion src/banded/infqltoeplitz.jl
Original file line number Diff line number Diff line change
Expand Up @@ -191,4 +191,4 @@ ql(A::InfToeplitz; kwds...) = _inf_ql(A; kwds...)
ql(A::PertToeplitz; kwds...) = _inf_ql(A; kwds...)

ql(A::Adjoint{<:Any,<:InfToeplitz}) = ql(BandedMatrix(A))
ql(A::Adjoint{<:Any,<:PertToeplitz}) = ql(BandedMatrix(A))
ql(A::Adjoint{<:Any,<:PertToeplitz}) = ql(BandedMatrix(A))
2 changes: 1 addition & 1 deletion src/blockbanded/infblocktridiagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -86,4 +86,4 @@ function BlockBandedMatrix(A::BlockTriPertToeplitz{T}, (l,u)::NTuple{2,Int}) whe
B = _BlockSkylineMatrix(Vcat(data,tl), BlockSkylineSizes(A, (l,u)))
end

BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A))
BlockBandedMatrix(A::BlockTriPertToeplitz) = BlockBandedMatrix(A, blockbandwidths(A))
4 changes: 2 additions & 2 deletions src/infcholesky.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ MemoryLayout(::Type{AdaptiveCholeskyFactors{T,DM,M}}) where {T,DM,M} = AdaptiveL


function partialcholesky!(F::AdaptiveCholeskyFactors{T,<:BandedMatrix}, n::Int) where T
if n > F.ncols
if n > F.ncols
_,u = bandwidths(F.data.array)
resizedata!(F.data,n+u,n+u);
ncols = F.ncols
Expand Down Expand Up @@ -106,4 +106,4 @@ function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveCholeskyFactors}, B::CachedVec
partialcholesky!(parent(R), n)
ldiv!(UpperTriangular(view(parent(R).data.data,oneto(n),oneto(n))), view(B.data,oneto(n)))
B
end
end
4 changes: 2 additions & 2 deletions src/infconv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,11 @@ end


function conv(r::InfRanges, x::AbstractVector)
length(x) ≠ 1 && throw(ArgumentError("conv(::$(typeof(r)), ::$(typeof(x))) not implemented"))
length(x) ≠ 1 && throw(ArgumentError("conv(::$(typeof(r)), ::$(typeof(x))) not implemented"))
first(x)*r
end
function conv(x::AbstractVector, r::InfRanges)
length(x) ≠ 1 && throw(ArgumentError("conv(::$(typeof(r)), ::$(typeof(x))) not implemented"))
length(x) ≠ 1 && throw(ArgumentError("conv(::$(typeof(r)), ::$(typeof(x))) not implemented"))
first(x)*r
end

Expand Down
30 changes: 15 additions & 15 deletions src/infql.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@ function tailiterate!(X::AbstractMatrix{T}) where T
c,a,b = X[1,:]
h = zero(T)
for _=1:10_000_000
QL = ql!(X)
if h == X[1,3]
QL = ql!(X)
if h == X[1,3]
return QL
end
h = X[1,3]
Expand All @@ -52,17 +52,17 @@ function qltail(Z::Number, A::Number, B::Number)

e = sqrt(n^2 - abs2(B))
d = σ*e*Z/n

ql!([Z A B; 0 d e])
end


ql_hessenberg(A::InfBandedMatrix{T}; kwds...) where T = ql_hessenberg!(BandedMatrix(A, (bandwidth(A,1)+bandwidth(A,2),bandwidth(A,2))); kwds...)

toeptail(B::BandedMatrix{T}) where T =
toeptail(B::BandedMatrix{T}) where T =
_BandedMatrix(B.data.args[end].args[1][1:end-B.u]*Ones{T}(1,∞), size(B,1), B.l-B.u, B.u)

# asymptotics of A[:,j:end] as j -> ∞
# asymptotics of A[:,j:end] as j -> ∞
rightasymptotics(d::Hcat) = last(d.args)
rightasymptotics(d::Vcat) = Vcat(rightasymptotics.(d.args)...)
rightasymptotics(d) = d
Expand All @@ -81,12 +81,12 @@ function ql_hessenberg!(B::InfBandedMatrix{TT}; kwds...) where TT
B̃[end,end] = L∞[1,1]
B̃[end:end,end-l+1:end-1] = adjoint(Q∞)[1:1,1:l-1]*T[l:2(l-1),1:l-1]
F = ql!(B̃)
# fill in data with L∞

# fill in data with L∞
B̃ = _BandedMatrix(B̃.data, size(data,2)+l, l,u)
B̃[size(data,2)+1:end,end-l+1:end] = adjoint(Q∞)[2:l+1,1:l+1] * T[l:2l,1:l]


# combine finite and infinite data
H = Hcat(B̃.data, rightasymptotics(F∞.factors.data))
QLHessenberg(_BandedMatrix(H, ℵ₀, l, 1), Vcat( LowerHessenbergQ(F.Q).q, F∞.q))
Expand All @@ -103,7 +103,7 @@ nzzeros(A::AbstractArray, k) = size(A,k)
nzzeros(::Zeros, k) = 0
nzzeros(B::Vcat, k) = sum(size.(B.args[1:end-1],k))
nzzeros(B::CachedArray, k) = max(B.datasize[k], nzzeros(B.array,k))
function nzzeros(B::AbstractMatrix, k)
function nzzeros(B::AbstractMatrix, k)
l,u = bandwidths(B)
k == 1 ? size(B,2) + l : size(B,1) + u
end
Expand Down Expand Up @@ -176,7 +176,7 @@ end
function _lmul_cache(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}
TS = promote_op(matprod, T, S)
lmul!(A, cache(convert(AbstractVector{TS},x)))
end
end

(*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::AbstractVector) where {T} = _lmul_cache(A, x)
(*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector) where {T} = _lmul_cache(A, x)
Expand All @@ -190,7 +190,7 @@ function blocktailiterate(c,a,b, d=c, e=a)
X = [c a b; z d e]
F = ql!(X)
d̃,ẽ = F.L[1:2,1:2], F.L[1:2,3:4]

d̃,ẽ = QLPackedQ(F.factors[1:2,3:4],F.τ[1:2])*d̃,QLPackedQ(F.factors[1:2,3:4],F.τ[1:2])*ẽ # undo last rotation
if ≈(d̃, d; atol=1E-10) && ≈(ẽ, e; atol=1E-10)
X[1:2,1:2] = d̃; X[1:2,3:4] = ẽ
Expand All @@ -211,7 +211,7 @@ function _blocktripert_ql(A, d, e)
c,a,b = A[Block(N+1,N)],A[Block(N,N)],A[Block(N-1,N)]
P,τ = blocktailiterate(c,a,b,d,e)
B = BlockBandedMatrix(A,(2,1))


BB = _BlockBandedMatrix(B.data.args[1], (fill(2,N+2), fill(2,N)), (2,1))
BB[Block(N),Block.(N-1:N)] .= P[Block(1), Block.(1:2)]
Expand Down Expand Up @@ -290,7 +290,7 @@ function materialize!(M::MatLdivVec{<:TriangularLayout{'L','N',BandedColumns{Per
@inbounds for j in 1:n
iszero(data[j,j]) && throw(SingularException(j))
bj = b[j] = data[j,j] \ b[j]
allzero = j > nz && iszero(bj)
allzero = j > nz && iszero(bj)
for i in (j+1:n) ∩ colsupport(data,j)
b[i] -= data[i,j] * bj
allzero = allzero && iszero(b[i])
Expand All @@ -316,7 +316,7 @@ function _ql(::SymTridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds
T = eltype(A)
d,d∞ = _data_tail(A.dv)
ev,ev∞ = _data_tail(A.ev)

m = max(length(d), length(ev)+1)
dat = zeros(T, 3, m)
dat[1,2:1+length(ev)] .= ev
Expand All @@ -338,7 +338,7 @@ function _ql(::TridiagonalLayout, ::NTuple{2,OneToInf{Int}}, A, args...; kwds...
d,d∞ = _data_tail(A.d)
dl,dl∞ = _data_tail(A.dl)
du,du∞ = _data_tail(A.du)

m = max(length(d), length(du)+1, length(dl))
dat = zeros(T, 3, m)
dat[1,2:1+length(du)] .= du
Expand Down
Loading