diff --git a/src/ApproxFunBase.jl b/src/ApproxFunBase.jl index edb7da60..2a78bd7b 100644 --- a/src/ApproxFunBase.jl +++ b/src/ApproxFunBase.jl @@ -70,7 +70,7 @@ import BlockBandedMatrices: blockbandwidth, blockbandwidths, blockcolstop, block subblockbandwidth, subblockbandwidths, _BlockBandedMatrix, _BandedBlockBandedMatrix, BandedBlockBandedMatrix, BlockBandedMatrix, isblockbanded, isbandedblockbanded, bb_numentries, BlockBandedSizes, - BandedBlockBandedSizes + BandedBlockBandedSizes, DefaultBandedBlockBandedMatrix import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_chebyshevtransform, plan_chebyshevtransform!, plan_ichebyshevtransform, plan_ichebyshevtransform! diff --git a/src/Caching/almostbanded.jl b/src/Caching/almostbanded.jl index 5205e80e..35c81afa 100644 --- a/src/Caching/almostbanded.jl +++ b/src/Caching/almostbanded.jl @@ -1,17 +1,19 @@ ## Caches -function CachedOperator(io::InterlaceOperator{T,1};padding::Bool=false) where T - ds=domainspace(io) - rs=rangespace(io) +CachedOperator(io::InterlaceOperator{T,1};padding::Bool=false) where T = + _interlace_CachedOperator(io, padding) - ind=findall(op->isinf(size(op,1)), io.ops) - if length(ind) ≠ 1 || !isbanded(io.ops[ind[1]]) # is almost banded +function _interlace_CachedOperator(io::InterlaceOperator{T,1}, padding) where T + ds = domainspace(io) + rs = rangespace(io) + + ind = findall(op->isinf(size(op,1)), io.ops) + if length(ind) ≠ 1 || !(isbanded(io.ops[ind[1]])::Bool) # is almost banded return default_CachedOperator(io;padding=padding) end i=ind[1] bo=io.ops[i] - lin,uin=bandwidths(bo) - + lin,uin=bandwidths(bo)::NTuple{2,Int} # calculate number of rows interlaced @@ -20,13 +22,12 @@ function CachedOperator(io::InterlaceOperator{T,1};padding::Bool=false) where T md=0 for k=1:length(io.ops) if k ≠ i - d=dimension(rs[k]) + d=dimension(rs[k])::Int nds+=d md=max(md,d) end end - isend=true for k=i+1:length(io.ops) if dimension(rs[k]) == md @@ -34,10 +35,10 @@ function CachedOperator(io::InterlaceOperator{T,1};padding::Bool=false) where T end end - numoprows=isend ? md-1 : md - n=nds+numoprows + numoprows=(isend ? md-1 : md) ::Int + n=(nds+numoprows) ::Int - (l,u) = (max(lin+nds,n-1),max(0,uin+1-ind[1])) + (l,u) = (max(lin+nds,n-1),max(0,uin+1-ind[1]))::Tuple{Int,Int} # add extra rows for QR if padding @@ -49,13 +50,13 @@ function CachedOperator(io::InterlaceOperator{T,1};padding::Bool=false) where T # populate the finite rows jr=1:n+u - ioM=io[1:n,jr] + ioM=io[1:n,jr]::RaggedMatrix{T} bcrow=1 oprow=0 for k=1:n - K,J=io.rangeinterlacer[k] + K,J= io.rangeinterlacer[k] if K ≠ i # fill the fill matrix @@ -173,12 +174,15 @@ function CachedOperator(io::InterlaceOperator{T,2};padding::Bool=false) where T end +resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}}, row::Integer, col::Integer) where {T} = + _almostbanded_resizedata!(co, co.op, row, col) + +resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}}, row::Integer, ::Colon) where {T} = + _almostbanded_resizedata!(co, co.op, row, :) # Grow cached interlace operator -function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, - InterlaceOperator{T,1,DS,RS,DI,RI,BI}}, - n::Integer,::Colon) where {T<:Number,DS,RS,DI,RI,BI} +function _almostbanded_resizedata!(co, op::InterlaceOperator{T,1}, n::Integer,::Colon) where {T<:Number} if n ≤ co.datasize[1] return co end @@ -187,12 +191,12 @@ function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, pad!(co.data,n,n+u) r=rank(co.data.fill) - ind=findfirst(op->isinf(size(op,1)),co.op.ops) + ind=findfirst(op->isinf(size(op,1)),op.ops) k=1 - for (K,J) in co.op.rangeinterlacer + for (K,J) in op.rangeinterlacer if K ≠ ind - co.data.fill.V[co.datasize[2]:end,k] = co.op.ops[K][J,co.datasize[2]:n+u] + co.data.fill.V[co.datasize[2]:end,k] = op.ops[K][J,co.datasize[2]:n+u] k += 1 if k > r break @@ -202,7 +206,7 @@ function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, kr=co.datasize[1]+1:n jr=max(1,kr[1]-l):n+u - BLAS.axpy!(1.0,view(co.op.ops[ind],kr .- r,jr), + BLAS.axpy!(1.0,view(op.ops[ind],kr .- r,jr), view(co.data.bands,kr,jr)) co.datasize=(n,n+u) @@ -211,14 +215,11 @@ end -function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, - InterlaceOperator{T,2,DS,RS,DI,RI,BI}}, - n::Integer,::Colon) where {T<:Number,DS,RS,DI,RI,BI} +function _almostbanded_resizedata!(co, io::InterlaceOperator{T,2}, n::Integer,::Colon) where T<:Number if n ≤ co.datasize[1] return co end - io=co.op ds=domainspace(io) rs=rangespace(io) di=io.domaininterlacer @@ -242,7 +243,7 @@ function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, K=k=1 while k ≤ r if isfinite(dimension(rs[ri[K][1]])) - co.data.fill.V[co.datasize[2]:end,k] = co.op[K,co.datasize[2]:n+u] + co.data.fill.V[co.datasize[2]:end,k] = io[K,co.datasize[2]:n+u] k += 1 end K += 1 @@ -259,14 +260,12 @@ function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, end -resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, - InterlaceOperator{T,1,DS,RS,DI,RI,BI}}, -n::Integer,m::Integer) where {T<:Number,DS,RS,DI,RI,BI} = resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:) +_almostbanded_resizedata!(co, ::InterlaceOperator{T,1}, n::Integer,m::Integer) where {T<:Number} = + resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:) -resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T}, - InterlaceOperator{T,2,DS,RS,DI,RI,BI}}, -n::Integer,m::Integer) where {T<:Number,DS,RS,DI,RI,BI} = resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:) +_almostbanded_resizedata!(co, ::InterlaceOperator{T,2}, n::Integer,m::Integer) where {T<:Number} = + resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:) @@ -285,8 +284,8 @@ end function resizedata!(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T}, - MM,DS,RS,BI}}, - ::Colon,col) where {T,MM,DS,RS,BI} + DS,RS,BI}}, + ::Colon,col) where {T,DS,RS,BI} if col ≤ QR.ncols return QR end @@ -351,8 +350,8 @@ end # BLAS versions, requires BlasFloat function resizedata!(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T}, - MM,DS,RS,BI}}, -::Colon,col) where {T<:BlasFloat,MM,DS,RS,BI} + DS,RS,BI}}, +::Colon,col) where {T<:BlasFloat,DS,RS,BI} if col ≤ QR.ncols return QR end diff --git a/src/Caching/banded.jl b/src/Caching/banded.jl index 9a22e6be..fb659a58 100644 --- a/src/Caching/banded.jl +++ b/src/Caching/banded.jl @@ -41,8 +41,8 @@ end function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T}, - MM,DS,RS,BI}}, - ::Colon,col) where {T,MM,DS,RS,BI} + DS,RS,BI}}, + ::Colon,col) where {T,DS,RS,BI} if col ≤ QR.ncols return QR end @@ -94,8 +94,8 @@ end function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T}, - MM,DS,RS,BI}}, - ::Colon,col) where {T<:BlasFloat,MM,DS,RS,BI} + DS,RS,BI}}, + ::Colon,col) where {T<:BlasFloat,DS,RS,BI} if col ≤ QR.ncols return QR end diff --git a/src/Caching/blockbanded.jl b/src/Caching/blockbanded.jl index cf284054..9409edb4 100644 --- a/src/Caching/blockbanded.jl +++ b/src/Caching/blockbanded.jl @@ -195,13 +195,13 @@ QROperator(R::CachedOperator{T,BlockBandedMatrix{T}}) where {T} = # always resize by column resizedata!(QR::QROperator{CachedOperator{T,BlockBandedMatrix{T}, - MM,DS,RS,BI}}, - ::Colon, col::Int) where {T,MM,DS,RS,BI} = + DS,RS,BI}}, + ::Colon, col::Int) where {T,DS,RS,BI} = resizedata!(QR, :, block(domainspace(QR.R_cache),col)) function resizedata!(QR::QROperator{CachedOperator{T,BlockBandedMatrix{T}, - MM,DS,RS,BI}}, - ::Colon, COL::Block) where {T<:BlasFloat,MM,DS,RS,BI} + DS,RS,BI}}, + ::Colon, COL::Block) where {T<:BlasFloat,DS,RS,BI} MO = QR.R_cache W = QR.H R = MO.data diff --git a/src/Caching/matrix.jl b/src/Caching/matrix.jl index 18f0636e..0eacbfae 100644 --- a/src/Caching/matrix.jl +++ b/src/Caching/matrix.jl @@ -1,5 +1,5 @@ CachedOperator(::Type{Matrix},op::Operator;padding::Bool=false) = - CachedOperator(op,Array{eltype(op)}(0,0),padding) + CachedOperator(op,Array{eltype(op)}(undef,0,0),padding) # Grow cached operator diff --git a/src/Caching/ragged.jl b/src/Caching/ragged.jl index 2825dcc9..4339974b 100644 --- a/src/Caching/ragged.jl +++ b/src/Caching/ragged.jl @@ -64,8 +64,8 @@ QROperator(R::CachedOperator{T,RaggedMatrix{T}}) where {T} = QROperator(R,RaggedMatrix{T}(undef,0,Int[]),0) function resizedata!(QR::QROperator{CachedOperator{T,RaggedMatrix{T}, - MM,DS,RS,BI}}, - ::Colon,col) where {T,MM,DS,RS,BI} + DS,RS,BI}}, + ::Colon,col) where {T,DS,RS,BI} if col ≤ QR.ncols return QR end @@ -127,8 +127,8 @@ end function resizedata!(QR::QROperator{CachedOperator{T,RaggedMatrix{T}, - MM,DS,RS,BI}}, -::Colon,col) where {T<:BlasFloat,MM,DS,RS,BI} + DS,RS,BI}}, +::Colon,col) where {T<:BlasFloat,DS,RS,BI} if col ≤ QR.ncols return QR end diff --git a/src/Fun.jl b/src/Fun.jl index e6be4613..11b0a12c 100644 --- a/src/Fun.jl +++ b/src/Fun.jl @@ -20,7 +20,6 @@ end const VFun{S,T} = Fun{S,T,Vector{T}} Fun(sp::Space,coeff::AbstractVector) = Fun{typeof(sp),eltype(coeff),typeof(coeff)}(sp,coeff) -Fun() = Fun(identity) Fun(d::Domain) = Fun(identity,d) Fun(d::Space) = Fun(identity,d) diff --git a/src/LinearAlgebra/AlmostBandedMatrix.jl b/src/LinearAlgebra/AlmostBandedMatrix.jl index 914b1349..e66c04b0 100644 --- a/src/LinearAlgebra/AlmostBandedMatrix.jl +++ b/src/LinearAlgebra/AlmostBandedMatrix.jl @@ -3,7 +3,7 @@ struct AlmostBandedMatrix{T} <: AbstractMatrix{T} - bands::BandedMatrix{T} + bands::BandedMatrix{T,Matrix{T},Base.OneTo{Int}} fill::LowRankMatrix{T} function AlmostBandedMatrix{T}(bands::BandedMatrix{T}, fill::LowRankMatrix{T}) where T if size(bands) ≠ size(fill) diff --git a/src/LinearAlgebra/RaggedMatrix.jl b/src/LinearAlgebra/RaggedMatrix.jl index 65fea5d3..9de093c9 100644 --- a/src/LinearAlgebra/RaggedMatrix.jl +++ b/src/LinearAlgebra/RaggedMatrix.jl @@ -28,8 +28,14 @@ end RaggedMatrix(dat::Vector,cols::Vector{Int},m::Int) = RaggedMatrix{eltype(dat)}(dat,cols,m) -RaggedMatrix{T}(::UndefInitializer, m::Int, colns::AbstractVector{Int}) where {T} = - RaggedMatrix(Vector{T}(undef, sum(colns)),Int[1;1 .+ cumsum(colns)],m) +function RaggedMatrix{T}(::UndefInitializer, m::Int, colns::AbstractVector{Int}) where {T} + cs = Vector{Int}(undef, length(colns)+1) + cs[1] = 1 + for j=2:length(cs) + cs[j] = cs[j-1] + colns[j-1] + end + RaggedMatrix(Vector{T}(undef, cs[end]),cs,m) +end Base.size(A::RaggedMatrix) = (A.m,length(A.cols)-1) diff --git a/src/LinearAlgebra/helper.jl b/src/LinearAlgebra/helper.jl index 13771ec0..38e53279 100644 --- a/src/LinearAlgebra/helper.jl +++ b/src/LinearAlgebra/helper.jl @@ -733,7 +733,7 @@ end const TrivialInterlacer{d} = BlockInterlacer{NTuple{d,<:Ones}} -BlockInterlacer(v::AbstractVector) = BlockInterlacer(tuple(v...)) +BlockInterlacer(v::AbstractVector) = uninfer(BlockInterlacer(tuple(v...))) Base.eltype(it::BlockInterlacer) = Tuple{Int,Int} @@ -792,4 +792,9 @@ function iterate(it::BlockInterlacer, (N,k,blkst,lngs)) return (N,lngs[N]),(N,k+1,blkst,lngs) end -cache(Q::BlockInterlacer) = CachedIterator(Q) \ No newline at end of file +cache(Q::BlockInterlacer) = CachedIterator(Q) + + + + +@noinline uninfer(@nospecialize(x)) = Ref{Any}(x)[] \ No newline at end of file diff --git a/src/Operators/Operator.jl b/src/Operators/Operator.jl index d5de7981..8782c10e 100644 --- a/src/Operators/Operator.jl +++ b/src/Operators/Operator.jl @@ -688,9 +688,13 @@ BlockBandedMatrix(::Type{Zeros}, V::Operator) = (AbstractVector{Int}(blocklengths(rangespace(V))), AbstractVector{Int}(blocklengths(domainspace(V)))), blockbandwidths(V)) -RaggedMatrix(::Type{Zeros}, V::Operator) = - RaggedMatrix(Zeros{eltype(V)}(size(V)), - Int[max(0,colstop(V,j)) for j=1:size(V,2)]) +function RaggedMatrix(::Type{Zeros}, V::Operator) + cs = Vector{Int}(undef, size(V,2)) + for j = 1:length(cs) + cs[j] = max(0,colstop(V,j)) + end + RaggedMatrix(Zeros{eltype(V)}(size(V)), cs) +end convert_axpy!(::Type{MT}, S::Operator) where {MT <: AbstractMatrix} = @@ -772,10 +776,10 @@ end function arraytype(V::SubOperator) P = parent(V) - isbanded(P) && return BandedMatrix + isbanded(P) && return BandedMatrix |> uninfer # isbandedblockbanded(P) && return BandedBlockBandedMatrix - isinf(size(P,1)) && israggedbelow(P) && return RaggedMatrix - return Matrix + isinf(size(P,1)) && israggedbelow(P) && return RaggedMatrix|> uninfer + return Matrix|> uninfer end AbstractMatrix(V::Operator) = arraytype(V)(V) diff --git a/src/Operators/banded/CalculusOperator.jl b/src/Operators/banded/CalculusOperator.jl index fcdb4891..a75e0a5d 100644 --- a/src/Operators/banded/CalculusOperator.jl +++ b/src/Operators/banded/CalculusOperator.jl @@ -18,8 +18,8 @@ macro calculus_operator(Op) space::S # the domain space order::OT end - struct $WrappOp{BT<:Operator,S<:Space,OT,T} <: $Op{S,OT,T} - op::BT + struct $WrappOp{S<:Space,OT,T} <: $Op{S,OT,T} + op::Operator{T} order::OT end @@ -60,7 +60,7 @@ macro calculus_operator(Op) end $WrappOp(op::Operator,order) = - $WrappOp{typeof(op),typeof(domainspace(op)),typeof(order),eltype(op)}(op,order) + $WrappOp{typeof(domainspace(op)),typeof(order),eltype(op)}(op,order) $WrappOp(op::Operator) = $WrappOp(op,1) function Base.convert(::Type{Operator{T}},D::$WrappOp) where T diff --git a/src/Operators/banded/Conversion.jl b/src/Operators/banded/Conversion.jl index e418e187..db7b5e8c 100644 --- a/src/Operators/banded/Conversion.jl +++ b/src/Operators/banded/Conversion.jl @@ -26,7 +26,7 @@ rangespace(C::ConcreteConversion)=C.rangespace -function defaultConversion(a::Space,b::Space)::Any +function defaultConversion(a::Space,b::Space) if a==b Conversion(a) elseif conversion_type(a,b)==NoSpace() @@ -40,7 +40,7 @@ function defaultConversion(a::Space,b::Space)::Any end else error("Implement Conversion from " * string(typeof(a)) * " to " * string(typeof(b))) - end + end |> uninfer end Conversion(a::Space,b::Space) = defaultConversion(a,b) @@ -53,16 +53,14 @@ Conversion() = ConversionWrapper(Operator(I,UnsetSpace())) # the domain and range space # but continue to know its a derivative -struct ConversionWrapper{S<:Operator,T} <: Conversion{T} - op::S +struct ConversionWrapper{T} <: Conversion{T} + op::Operator{T} end @wrapper ConversionWrapper -ConversionWrapper(::Type{T},op) where {T} = ConversionWrapper{typeof(op),T}(op) -ConversionWrapper(B::Operator) = - ConversionWrapper{typeof(B),eltype(B)}(B) +ConversionWrapper(::Type{T},op) where {T} = ConversionWrapper{T}(op) Conversion(A::Space,B::Space,C::Space) = ConversionWrapper(Conversion(B,C)*Conversion(A,B)) Conversion(A::Space,B::Space,C::Space,D::Space...) = @@ -76,7 +74,7 @@ function convert(::Type{Operator{T}},D::ConversionWrapper) where T D else BO=convert(Operator{T},D.op) - ConversionWrapper{typeof(BO),T}(BO) + ConversionWrapper{T}(BO) end end diff --git a/src/Operators/banded/Multiplication.jl b/src/Operators/banded/Multiplication.jl index 29a117d4..cc0602f2 100644 --- a/src/Operators/banded/Multiplication.jl +++ b/src/Operators/banded/Multiplication.jl @@ -88,17 +88,17 @@ choosedomainspace(M::ConcreteMultiplication{D,UnsetSpace},sp::Space) where {D} = diagm(a::Fun) = Multiplication(a) -struct MultiplicationWrapper{D<:Space,S<:Space,O<:Operator,T} <: Multiplication{D,S,T} +struct MultiplicationWrapper{D<:Space,S<:Space,T} <: Multiplication{D,S,T} f::VFun{D,T} - op::O + op::Operator{T} end -MultiplicationWrapper(T::Type,f::Fun{D,V},op::Operator) where {D<:Space,V} = MultiplicationWrapper{D,typeof(domainspace(op)),typeof(op),T}(f,op) +MultiplicationWrapper(T::Type,f::Fun{D,V},op::Operator) where {D<:Space,V} = MultiplicationWrapper{D,typeof(domainspace(op)),T}(f,op) MultiplicationWrapper(f::Fun{D,V},op::Operator) where {D<:Space,V} = MultiplicationWrapper(eltype(op),f,op) @wrapper MultiplicationWrapper -function convert(::Type{Operator{TT}},C::MultiplicationWrapper{S,V,O,T}) where {TT,S,V,O,T} +function convert(::Type{Operator{TT}},C::MultiplicationWrapper{S,V,T}) where {TT,S,V,T} if TT==T C else diff --git a/src/Operators/banded/Reverse.jl b/src/Operators/banded/Reverse.jl index d4f725c8..1653d6ce 100644 --- a/src/Operators/banded/Reverse.jl +++ b/src/Operators/banded/Reverse.jl @@ -5,11 +5,10 @@ for TYP in (:ReverseOrientation,:Reverse) @eval begin abstract type $TYP{T} <: Operator{T} end - struct $WRAP{OS,T} <: Operator{T} - op::OS + struct $WRAP{T} <: Operator{T} + op::Operator{T} end - $WRAP(op::Operator) = $WRAP{typeof(op),eltype(op)}(op) convert(::Type{Operator{T}},op::$TYP) where {T} = $TYP{T}() convert(::Type{Operator{T}},op::$WRAP) where {T} = $WRAP(Operator{T}(op.op))::Operator{T} diff --git a/src/Operators/functionals/CalculusFunctional.jl b/src/Operators/functionals/CalculusFunctional.jl index e7318477..a0a6f700 100644 --- a/src/Operators/functionals/CalculusFunctional.jl +++ b/src/Operators/functionals/CalculusFunctional.jl @@ -14,8 +14,8 @@ macro calculus_functional(Op) struct $ConcOp{S,T} <: $Op{S,T} domainspace::S end - struct $WrappOp{BT<:Operator,S<:Space,T} <: $Op{S,T} - op::BT + struct $WrappOp{S<:Space,T} <: $Op{S,T} + op::Operator{T} end @wrapper $WrappOp diff --git a/src/Operators/functionals/Evaluation.jl b/src/Operators/functionals/Evaluation.jl index a5eb6b6c..95b4d6a3 100644 --- a/src/Operators/functionals/Evaluation.jl +++ b/src/Operators/functionals/Evaluation.jl @@ -20,7 +20,7 @@ ConcreteEvaluation(sp::Space,x,o::Number) = Evaluation(::Type{T},sp::Space,x,order) where {T} = ConcreteEvaluation{typeof(sp),typeof(x),typeof(order),T}(sp,x,order) # TODO: This seems like a bad idea: if you are specifying x, just go with the generic version -function Evaluation(::Type{T},sp::UnivariateSpace,x::Number,order) where {T} +function Evaluation(::Type{T}, @nospecialize(sp::UnivariateSpace), x::Number, order) where {T} d=domain(sp) if isa(d,IntervalOrSegment) && isapprox(leftendpoint(d),x) Evaluation(T,sp,leftendpoint,order) @@ -28,7 +28,7 @@ function Evaluation(::Type{T},sp::UnivariateSpace,x::Number,order) where {T} Evaluation(T,sp,rightendpoint,order) else ConcreteEvaluation{typeof(sp),typeof(x),typeof(order),T}(sp,x,order) - end + end |> uninfer end Evaluation(sp::Space,x,order) = Evaluation(rangetype(sp),sp,x,order) @@ -83,17 +83,17 @@ end ## EvaluationWrapper -struct EvaluationWrapper{S<:Space,M,FS<:Operator,OT,T<:Number} <: Evaluation{T} +struct EvaluationWrapper{S<:Space,M,OT,T<:Number} <: Evaluation{T} space::S x::M order::OT - op::FS + op::Operator{T} end @wrapper EvaluationWrapper EvaluationWrapper(sp::Space,x,order,func::Operator) = - EvaluationWrapper{typeof(sp),typeof(x),typeof(func),typeof(order),eltype(func)}(sp,x,order,func) + EvaluationWrapper{typeof(sp),typeof(x),typeof(order),eltype(func)}(sp,x,order,func) domainspace(E::Evaluation) = E.space @@ -166,13 +166,13 @@ convert(::Type{Operator{T}},B::ConcreteDirichlet{S,V}) where {S,V,T} = struct DirichletWrapper{S,T} <: Dirichlet{S,T} - op::S + op::Operator{T} order::Int end @wrapper DirichletWrapper -DirichletWrapper(B::Operator,λ=0) = DirichletWrapper{typeof(B),eltype(B)}(B,λ) +DirichletWrapper(B::Operator,λ=0) = DirichletWrapper{typeof(domainspace(B)),eltype(B)}(B,λ) convert(::Type{Operator{T}},B::DirichletWrapper) where {T} = DirichletWrapper(Operator{T}(B.op),B.order)::Operator{T} diff --git a/src/Operators/general/CachedOperator.jl b/src/Operators/general/CachedOperator.jl index 38905538..6c85479d 100644 --- a/src/Operators/general/CachedOperator.jl +++ b/src/Operators/general/CachedOperator.jl @@ -8,8 +8,8 @@ export cache ## CachedOperator -mutable struct CachedOperator{T<:Number,DM<:AbstractMatrix,M<:Operator,DS,RS,BI} <: Operator{T} - op::M +mutable struct CachedOperator{T<:Number,DM<:AbstractMatrix,DS,RS,BI} <: Operator{T} + op::Operator{T} data::DM datasize::Tuple{Int,Int} domainspace::DS @@ -21,8 +21,7 @@ mutable struct CachedOperator{T<:Number,DM<:AbstractMatrix,M<:Operator,DS,RS,BI} end CachedOperator(op::Operator,data::AbstractMatrix,sz::Tuple{Int,Int},ds,rs,bi,pd=false) = - CachedOperator{eltype(data),typeof(data),typeof(op), - typeof(ds),typeof(rs),typeof(bi)}(op,data,sz,ds,rs,bi,pd) + CachedOperator{eltype(data),typeof(data),typeof(ds),typeof(rs),typeof(bi)}(op,data,sz,ds,rs,bi,pd) |> uninfer CachedOperator(op::Operator,data::AbstractMatrix,sz::Tuple{Int,Int},pd=false) = @@ -31,7 +30,7 @@ CachedOperator(op::Operator,data::AbstractMatrix,padding=false) = CachedOperator -function default_CachedOperator(op::Operator;padding::Bool=false) +@inline function default_CachedOperator(op::Operator;padding::Bool=false) if isbanded(op) CachedOperator(BandedMatrix,op;padding=padding) elseif isbandedblockbanded(op) && !padding @@ -42,17 +41,18 @@ function default_CachedOperator(op::Operator;padding::Bool=false) CachedOperator(RaggedMatrix,op;padding=padding) else CachedOperator(Matrix,op;padding=padding) - end + end |> uninfer end -CachedOperator(op::Operator;padding::Bool=false) = default_CachedOperator(op;padding=padding) +@inline CachedOperator(op::Operator; padding::Bool=false) = + default_CachedOperator(op;padding=padding) |> uninfer """ cache(op::Operator) Caches the entries of an operator, to speed up multiplying a Fun by the operator. """ -cache(O::Operator;kwds...) = CachedOperator(O;kwds...) +@inline cache(O::Operator;kwds...) = CachedOperator(O;kwds...) cache(::Type{MT},O::Operator;kwds...) where {MT<:AbstractMatrix} = CachedOperator(MT,O;kwds...) convert(::Type{Operator{T}},S::CachedOperator{T}) where {T} = S @@ -63,8 +63,8 @@ convert(::Type{Operator{T}},S::CachedOperator) where {T} = domainspace(C::CachedOperator) = C.domainspace rangespace(C::CachedOperator) = C.rangespace -bandwidths(C::CachedOperator{T,BM,M}) where {T<:Number,BM<:BandedMatrix,M<:Operator} = C.bandwidths -blockbandwidths(C::CachedOperator{T,BM,M}) where {T<:Number,BM<:BandedMatrix,M<:Operator} = C.bandwidths +bandwidths(C::CachedOperator{T,BM}) where {T<:Number,BM<:BandedMatrix} = C.bandwidths +blockbandwidths(C::CachedOperator{T,BM}) where {T<:Number,BM<:BandedMatrix} = C.bandwidths # TODO: cache this information as well diff --git a/src/Operators/general/InterlaceOperator.jl b/src/Operators/general/InterlaceOperator.jl index 46d75f87..fbc121cc 100644 --- a/src/Operators/general/InterlaceOperator.jl +++ b/src/Operators/general/InterlaceOperator.jl @@ -345,7 +345,7 @@ for TYP in (:BandedMatrix, :BlockBandedMatrix, :BandedBlockBandedMatrix, :Ragged if !isempty(ret_kr) sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2] - LinearAlgebra.axpy!(1.0,view(L.ops[ν],sub_kr,jr),view(ret,ret_kr,:)) + LinearAlgebra.axpy!(one(T),view(L.ops[ν],sub_kr,jr),view(ret,ret_kr,:)) end end ret @@ -460,7 +460,7 @@ function _vcat(A::OperatorTypes...) push!(Av,a) end end - InterlaceOperator(vnocat(Av...)) + uninfer(InterlaceOperator(vnocat(Av...))) end diff --git a/src/Operators/general/algebra.jl b/src/Operators/general/algebra.jl index 3d87a259..940d0955 100644 --- a/src/Operators/general/algebra.jl +++ b/src/Operators/general/algebra.jl @@ -45,7 +45,7 @@ function PlusOperator(ops::Vector) b1=max(br[1],b1) b2=max(br[end],b2) end - PlusOperator(ops,(b1,b2)) + uninfer(PlusOperator(ops,(b1,b2))) end function convert(::Type{Operator{T}},P::PlusOperator) where T @@ -378,14 +378,12 @@ function getindex(P::TimesOperator,k::AbstractVector) vec(Matrix(P[1:1,k])) end -for TYP in (:Matrix, :BandedMatrix, :RaggedMatrix) +for TYP in (:Matrix, :RaggedMatrix) @eval function $TYP(V::SubOperator{T,TO,Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,TO<:TimesOperator} P = parent(V) if isbanded(P) - if $TYP ≠ BandedMatrix - return $TYP(BandedMatrix(V)) - end + return $TYP(BandedMatrix(V)) elseif isbandedblockbanded(P) N = block(rangespace(P), last(parentindices(V)[1])) M = block(domainspace(P), last(parentindices(V)[2])) @@ -438,9 +436,9 @@ for TYP in (:Matrix, :BandedMatrix, :RaggedMatrix) # The following returns a banded Matrix with all rows # for large k its upper triangular - BA = convert($TYP{T}, P.ops[end][krl[end,1]:krl[end,2],jr]) + BA = uninfer(convert($TYP{T}, P.ops[end][krl[end,1]:krl[end,2],jr])) for m = (length(P.ops)-1):-1:1 - BA = convert($TYP{T}, P.ops[m][krl[m,1]:krl[m,2],krl[m+1,1]:krl[m+1,2]])*BA + BA = uninfer(convert($TYP{T}, P.ops[m][krl[m,1]:krl[m,2],krl[m+1,1]:krl[m+1,2]])*BA) end $TYP{T}(BA) @@ -502,6 +500,71 @@ for TYP in (:BlockBandedMatrix, :BandedBlockBandedMatrix) end end +function _krlin(ops, kr, jr) + krlin = Matrix{Int}(undef,length(ops),2) + + krlin[1,1],krlin[1,2]=kr[1],kr[end] + for m=1:length(ops)-1 + krlin[m+1,1]=rowstart(ops[m],krlin[m,1]) + krlin[m+1,2]=rowstop(ops[m],krlin[m,2]) + end + krlin[end,1]=max(krlin[end,1],colstart(ops[end],jr[1])) + krlin[end,2]=min(krlin[end,2],colstop(ops[end],jr[end])) + for m=length(ops)-1:-1:2 + krlin[m,1]=max(krlin[m,1],colstart(ops[m],krlin[m+1,1])) + krlin[m,2]=min(krlin[m,2],colstop(ops[m],krlin[m+1,2])) + end + krlin +end + +function BandedMatrix(V::SubOperator{T,TO,Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,TO<:TimesOperator} + P = parent(V) + + if isbandedblockbanded(P) + N = block(rangespace(P), last(parentindices(V)[1]))::Block{1,Int64} + M = block(domainspace(P), last(parentindices(V)[2]))::Block{1,Int64} + B = P[Block(1):N, Block(1):M]::DefaultBandedBlockBandedMatrix{T} + return BandedMatrix(view(B, parentindices(V)...), _colstops(V)) + end + + kr,jr = parentindices(V) + + (isempty(kr) || isempty(jr)) && return BandedMatrix(Zeros, V) + + if maximum(kr) > size(P,1) || maximum(jr) > size(P,2) || + minimum(kr) < 1 || minimum(jr) < 1 + throw(BoundsError()) + end + + @assert length(P.ops) ≥ 2 + if size(V,1)==0 + return BandedMatrix(Zeros, V) + end + + + # find optimal truncations for each operator + # by finding the non-zero entries + krl = _krlin(P.ops, kr, jr) + + # Check if any range is invalid, in which case return zero + for m=1:length(P.ops) + if krl[m,1]>krl[m,2] + return BandedMatrix(Zeros, V) + end + end + + + + # The following returns a banded Matrix with all rows + # for large k its upper triangular + BA = convert(BandedMatrix{T,Matrix{T},Base.OneTo{Int}}, P.ops[end][krl[end,1]:krl[end,2],jr]) + for m = (length(P.ops)-1):-1:1 + BA = convert(BandedMatrix{T,Matrix{T},Base.OneTo{Int}}, P.ops[m][krl[m,1]:krl[m,2],krl[m+1,1]:krl[m+1,2]])*BA + end + + convert(BandedMatrix{T,Matrix{T},Base.OneTo{Int}},BA) +end + ## Algebra: assume we promote @@ -540,11 +603,11 @@ end # Conversions we always assume are intentional: no need to promote -*(A::ConversionWrapper{TO1},B::ConversionWrapper{TO}) where {TO1<:TimesOperator,TO<:TimesOperator} = +*(A::ConversionWrapper,B::ConversionWrapper) = ConversionWrapper(TimesOperator(A.op,B.op)) -*(A::ConversionWrapper{TO},B::Conversion) where {TO<:TimesOperator} = +*(A::ConversionWrapper, B::Conversion) = ConversionWrapper(TimesOperator(A.op,B)) -*(A::Conversion,B::ConversionWrapper{TO}) where {TO<:TimesOperator} = +*(A::Conversion,B::ConversionWrapper) = ConversionWrapper(TimesOperator(A,B.op)) *(A::Conversion,B::Conversion) = ConversionWrapper(TimesOperator(A,B)) diff --git a/src/Operators/ldiv.jl b/src/Operators/ldiv.jl index abbbb090..ceb74faa 100644 --- a/src/Operators/ldiv.jl +++ b/src/Operators/ldiv.jl @@ -14,6 +14,8 @@ for TYP in (:Fun,:StridedVector,:AbstractVector,:Any) end end + + """ \\(A::Operator,b;tolerance=tol,maxlength=n) diff --git a/src/Operators/spacepromotion.jl b/src/Operators/spacepromotion.jl index 7490e332..dd32f5d1 100644 --- a/src/Operators/spacepromotion.jl +++ b/src/Operators/spacepromotion.jl @@ -47,27 +47,26 @@ rangespace(S::SpaceOperator) = S.rangespace ##TODO: Do we need both max and min? -function findmindomainspace(ops::AbstractVector)::Any +function findmindomainspace(ops::AbstractVector) sp = UnsetSpace() for op in ops sp = union(sp,domainspace(op)) end - sp + sp |> uninfer end -function findmaxrangespace(ops::AbstractVector)::Any +function findmaxrangespace(ops::AbstractVector) sp = UnsetSpace() for op in ops sp = maxspace(sp,rangespace(op)) end - sp + sp |> uninfer end - # The coolest definitions ever!! # supports Derivative():Chebyshev()→Ultraspherical(1) (:)(A::Operator,b::Space) = promotedomainspace(A,b) @@ -79,12 +78,10 @@ end promoterangespace(P::Operator,sp::Space) = promoterangespace(P,sp,rangespace(P)) promotedomainspace(P::Operator,sp::Space) = promotedomainspace(P,sp,domainspace(P)) - -promoterangespace(P::Operator,sp::Space,cursp::Space)::Any = - (sp==cursp) ? P : Conversion(cursp,sp)*P -promotedomainspace(P::Operator,sp::Space,cursp::Space)::Any = - (sp==cursp) ? P : P*Conversion(sp,cursp) - +promoterangespace(P::Operator,sp::Space,cursp::Space) = + ((sp==cursp) ? P : Conversion(cursp,sp)*P) |> uninfer +promotedomainspace(P::Operator,sp::Space,cursp::Space) = + ((sp==cursp) ? P : P*Conversion(sp,cursp)) |> uninfer @@ -134,27 +131,28 @@ choosedomainspace(A::Operator,_) = choosedomainspace(A) choosedomainspace(A) = choosedomainspace(A,UnsetSpace()) -function choosedomainspace(ops::AbstractVector,spin)::Any + +function choosedomainspace(@nospecialize(ops::AbstractVector), @nospecialize(spin)) sp = UnsetSpace() for op in ops sp = conversion_type(sp,choosedomainspace(op,spin)) end - sp + sp |> uninfer end -choosespaces(A::Operator,b) = promotedomainspace(A,choosedomainspace(A,b)) +choosespaces(A::Operator, b) = promotedomainspace(A,choosedomainspace(A,b)) -spacescompatible(A::Operator,B::Operator) = +spacescompatible(A::Operator, B::Operator) = spacescompatible(domainspace(A),domainspace(B)) && spacescompatible(rangespace(A),rangespace(B)) #It's important that domain space is promoted first as it might impact range space promotespaces(ops::AbstractVector) = promoterangespace(promotedomainspace(ops)) -function promotespaces(ops::AbstractVector,b::Fun) +function promotespaces(@nospecialize(ops::AbstractVector), @nospecialize(b::Fun)) A=promotespaces(ops) if isa(rangespace(A),AmbiguousSpace) # try setting the domain space @@ -164,7 +162,7 @@ function promotespaces(ops::AbstractVector,b::Fun) end -function promotespaces(A::Operator,B::Operator) +function promotespaces(@nospecialize(A::Operator), @nospecialize(B::Operator)) if spacescompatible(A,B) A,B else diff --git a/src/PDE/KroneckerOperator.jl b/src/PDE/KroneckerOperator.jl index 21fac4c2..8dd856d9 100644 --- a/src/PDE/KroneckerOperator.jl +++ b/src/PDE/KroneckerOperator.jl @@ -149,7 +149,7 @@ const Wrappers = Union{ConversionWrapper,MultiplicationWrapper,DerivativeWrapper -isbandedblockbanded(P::Union{PlusOperator,TimesOperator}) = all(isbandedblockbanded,P.ops) +isbandedblockbanded(P::Union{PlusOperator,TimesOperator}) = all(isbandedblockbanded,P.ops)::Bool diff --git a/src/Space.jl b/src/Space.jl index 2bde43da..cac23d95 100644 --- a/src/Space.jl +++ b/src/Space.jl @@ -147,7 +147,7 @@ mappoint(a::Domain,b::Space,x)=mappoint(a,domain(b),x) for FUNC in (:conversion_rule,:maxspace_rule,:union_rule) @eval begin - function $FUNC(a,b)::Any + function $FUNC(@nospecialize(a), @nospecialize(b)) if spacescompatible(a,b) a else @@ -169,7 +169,7 @@ end # gives a space c that has a banded conversion operator TO a and b -function conversion_type(a,b)::Any +function conversion_type(@nospecialize(a), @nospecialize(b)) if spacescompatible(a,b) a elseif !domainscompatible(a,b) @@ -188,7 +188,7 @@ end # gives a space c that has a banded conversion operator FROM a and b maxspace(a,b) = NoSpace() # TODO: this fixes weird bug with Nothing -function maxspace(a::Space, b::Space)::Any +function maxspace(@nospecialize(a::Space), @nospecialize(b::Space)) if spacescompatible(a,b) return a elseif !domainscompatible(a,b) @@ -248,7 +248,7 @@ union(a::AmbiguousSpace, b::Space) = b union(a::Space, b::AmbiguousSpace) = a -function union_by_union_rule(a::Space,b::Space)::Any +function union_by_union_rule(@nospecialize(a::Space), @nospecialize(b::Space)) if spacescompatible(a,b) if isambiguous(domain(a)) return b @@ -263,7 +263,7 @@ function union_by_union_rule(a::Space,b::Space)::Any union_rule(b,a) end -function union(a::Space, b::Space)::Any +function union(@nospecialize(a::Space), @nospecialize(b::Space)) cr = union_by_union_rule(a,b) cr isa NoSpace || return cr @@ -278,7 +278,7 @@ function union(a::Space, b::Space)::Any return cr end - a ⊕ b + a ⊕ b |> uninfer end union(a::Space) = a diff --git a/src/Spaces/ArraySpace.jl b/src/Spaces/ArraySpace.jl index cf59d860..486a2e53 100644 --- a/src/Spaces/ArraySpace.jl +++ b/src/Spaces/ArraySpace.jl @@ -34,7 +34,7 @@ Vector(sp::VectorSpace) = sp.spaces Matrix(sp::MatrixSpace) = sp.spaces -BlockInterlacer(sp::ArraySpace) = BlockInterlacer(blocklengths.(tuple(sp.spaces...))) +BlockInterlacer(sp::ArraySpace) = uninfer(BlockInterlacer(map(blocklengths,sp.spaces))) interlacer(sp::ArraySpace) = BlockInterlacer(sp) for OP in (:length,:firstindex,:lastindex,:size) diff --git a/src/Spaces/ProductSpaceOperators.jl b/src/Spaces/ProductSpaceOperators.jl index ccbe87c2..283782f6 100644 --- a/src/Spaces/ProductSpaceOperators.jl +++ b/src/Spaces/ProductSpaceOperators.jl @@ -5,14 +5,15 @@ export continuity for TYP in (:PiecewiseSpace,:ArraySpace) @eval begin - function promotedomainspace(A::InterlaceOperator{T,2},sp::$TYP)::Any where T + @nospecialize + function promotedomainspace(A::InterlaceOperator{T,2},sp::$TYP) where T if domainspace(A) == sp return A end @assert size(A.ops,2) == length(sp) InterlaceOperator([promotedomainspace(A.ops[k,j],sp[j]) for k=1:size(A.ops,1),j=1:size(A.ops,2)],$TYP) end - function interlace_choosedomainspace(ops,rs::$TYP)::Any + function interlace_choosedomainspace(ops,rs::$TYP) @assert length(ops) == length(rs) # this ensures correct dispatch for unino sps = Array{Space}( @@ -23,6 +24,7 @@ for TYP in (:PiecewiseSpace,:ArraySpace) union(sps...) end end + @specialize end end @@ -93,7 +95,7 @@ end # Sum Space and PiecewiseSpace need to allow permutation of space orders for TYP in (:SumSpace,:PiecewiseSpace) - @eval function Conversion(S1::$TYP,S2::$TYP)::Any + @eval function Conversion(@nospecialize(S1::$TYP), @nospecialize(S2::$TYP)) v1 = collect(S1.spaces) v2 = collect(S2.spaces) @@ -180,7 +182,7 @@ end for (OPrule,OP) in ((:conversion_rule,:conversion_type),(:maxspace_rule,:maxspace), (:union_rule,:union)) for TYP in (:SumSpace,:PiecewiseSpace) - @eval function $OPrule(S1sp::$TYP,S2sp::$TYP)::Any + @eval function $OPrule(@nospecialize(S1sp::$TYP), @nospecialize(S2sp::$TYP)) S1 = components(S1sp) S2 = components(S2sp) cs1,cs2=map(canonicalspace,S1),map(canonicalspace,S2) @@ -239,7 +241,7 @@ for (Op,OpWrap) in ((:Derivative,:DerivativeWrapper),(:Integral,:IntegralWrapper end end -function Derivative(S::SumSpace,k::Integer)::Any +function Derivative(@nospecialize(S::SumSpace), k::Integer) # we want to map before we decompose, as the map may introduce # mixed bases. if typeof(canonicaldomain(S))==typeof(domain(S)) diff --git a/src/Spaces/Spaces.jl b/src/Spaces/Spaces.jl index 7169e40a..4c9b2935 100644 --- a/src/Spaces/Spaces.jl +++ b/src/Spaces/Spaces.jl @@ -8,7 +8,7 @@ include("SubSpace.jl") include("QuotientSpace.jl") -⊕(A::Space,B::Space)::Any = domainscompatible(A,B) ? SumSpace(A,B) : PiecewiseSpace(A,B) +⊕(@nospecialize(A::Space), @nospecialize(B::Space)) = domainscompatible(A,B) ? SumSpace(A,B) : PiecewiseSpace(A,B) ⊕(f::Fun,g::Fun) = Fun(space(f) ⊕ space(g), interlace(coefficients(f),coefficients(g))) ⊕(f::Fun,g::Fun,h::Fun...) = ⊕((f ⊕ g), h...) diff --git a/src/hacks.jl b/src/hacks.jl index 337caa9c..32927a3a 100644 --- a/src/hacks.jl +++ b/src/hacks.jl @@ -11,6 +11,7 @@ end # linear algebra +Fun() = Fun(identity,ChebyshevInterval()) ## Constructors that involve MultivariateFun Fun(f::Fun) = f # Fun of Fun should be like a conversion diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl index a1a3af23..734b5fd6 100644 --- a/src/specialfunctions.jl +++ b/src/specialfunctions.jl @@ -173,6 +173,7 @@ for (op,ODE,RHS,growth) in ((:(exp),"D-f'","0",:(real)), # g=exp(Fun(f.coefficients,space(f).space)) # Fun(g.coefficients,MappedSpace(domain(f),space(g))) # end + @nospecialize function $op(fin::Fun{S,T}) where {S,T} f=setcanonicaldomain(fin) # removes possible issues with roots @@ -184,8 +185,9 @@ for (op,ODE,RHS,growth) in ((:(exp),"D-f'","0",:(real)), B=Evaluation(domainspace(D),xmax) u=\([B,eval($L)],Any[opfxmax,eval($R)];tolerance=eps(T)*opmax) - setdomain(u,domain(fin)) + setdomain(u,domain(fin)) end + @specialize end end