blockbandwidths in InterlaceOperator (#410)

* blockbandwidths in InterlaceOperator

* version bump to v0.8.7

* accept blockbandwidths in inner constructor

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunBase"
 uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
-version = "0.8.6"
+version = "0.8.7"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/Operators/Operator.jl

@@ -137,9 +137,9 @@ Base.ndims(::Operator) = 2
 
 
 ## bandrange and indexrange
-isbandedbelow(A::Operator) = isfinite(bandwidth(A,1))
-isbandedabove(A::Operator) = isfinite(bandwidth(A,2))
-isbanded(A::Operator) = isbandedbelow(A) && isbandedabove(A)
+isbandedbelow(A::Operator) = isfinite(bandwidth(A,1))::Bool
+isbandedabove(A::Operator) = isfinite(bandwidth(A,2))::Bool
+isbanded(A::Operator) = all(isfinite, bandwidths(A))::Bool
 
 
 isbandedblockbandedbelow(_) = false
@@ -178,11 +178,11 @@ end
 # assume dense blocks
 subblockbandwidths(K::Operator) = maximum(blocklengths(rangespace(K)))-1, maximum(blocklengths(domainspace(K)))-1
 
-isblockbandedbelow(A) = isfinite(blockbandwidth(A,1))
-isblockbandedabove(A) = isfinite(blockbandwidth(A,2))
-isblockbanded(A::Operator) = isblockbandedbelow(A) && isblockbandedabove(A)
+isblockbandedbelow(A::Operator) = isfinite(blockbandwidth(A,1))::Bool
+isblockbandedabove(A::Operator) = isfinite(blockbandwidth(A,2))::Bool
+isblockbanded(A::Operator) = all(isfinite, blockbandwidths(A))::Bool
 
-israggedbelow(A::Operator) = isbandedbelow(A) || isbandedblockbanded(A) || isblockbandedbelow(A)
+israggedbelow(A::Operator) = isbandedbelow(A)::Bool || isbandedblockbanded(A)::Bool || isblockbandedbelow(A)::Bool
 
 
 blockbandwidth(K::Operator, k::Integer) = blockbandwidths(K)[k]

src/Operators/general/InterlaceOperator.jl

@@ -69,13 +69,22 @@ end
 
 ## Interlace operator
 
-struct InterlaceOperator{T,p,DS,RS,DI,RI,BI} <: Operator{T}
+struct InterlaceOperator{T,p,DS,RS,DI,RI,BI,BBW} <: Operator{T}
     ops::Array{Operator{T},p}
     domainspace::DS
     rangespace::RS
     domaininterlacer::DI
     rangeinterlacer::RI
     bandwidths::BI
+    blockbandwidths::BBW
+    israggedbelow::Bool
+
+    function InterlaceOperator(ops::Array{Operator{T},p}, ds::DS, rs::RS, dsi::DI, rsi::RI, bw::BI,
+        blockbandwidths::BBW = bandwidthsmax(ops, blockbandwidths),
+        israggedbelow::Bool = all(israggedbelow, ops)) where {T,p,DS,RS,DI,RI,BI,BBW}
+
+        new{T,p,DS,RS,DI,RI,BI,BBW}(ops, ds, rs, dsi, rsi, bw, blockbandwidths, israggedbelow)
+    end
 end
 
 const VectorInterlaceOperator = InterlaceOperator{T,1,DS,RS} where {T,DS,RS<:Space{D,R}} where {D,R<:AbstractVector}
@@ -133,7 +142,9 @@ end
 
 function InterlaceOperator(ops::AbstractMatrix{<:Operator},ds::Space,rs::Space;
         # calculate bandwidths TODO: generalize
-        bandwidths = interlace_bandwidths(ops, ds, rs))
+        bandwidths = interlace_bandwidths(ops, ds, rs),
+        blockbandwidths = bandwidthsmax(ops, blockbandwidths),
+        israggedbelow = all(israggedbelow, ops))
 
     dsi = interlacer(ds)
     rsi = interlacer(rs)
@@ -143,7 +154,9 @@ function InterlaceOperator(ops::AbstractMatrix{<:Operator},ds::Space,rs::Space;
     InterlaceOperator(opsm,ds,rs,
                         cache(dsi),
                         cache(rsi),
-                        bandwidths)
+                        bandwidths,
+                        blockbandwidths,
+                        israggedbelow)
 end
 
 @inline function _interlace_bandwidths(ops::VectorOrTupleOfOp, ds, rs, allbanded = all(isbanded, ops))
@@ -165,14 +178,18 @@ end
 
 function InterlaceOperator(ops::VectorOrTupleOfOp, ds::Space, rs::Space;
         # calculate bandwidths
-        bandwidths = interlace_bandwidths(ops, ds, rs))
+        bandwidths = interlace_bandwidths(ops, ds, rs),
+        blockbandwidths = bandwidthsmax(ops, blockbandwidths),
+        israggedbelow = all(israggedbelow, ops))
 
     VT = Vector{Operator{promote_eltypeof(ops)}}
     opsv = strictconvert(VT, convert_vector(ops))
     InterlaceOperator(opsv,ds,rs,
                         cache(BlockInterlacer(tuple(blocklengths(ds)))),
                         cache(interlacer(rs)),
-                        bandwidths)
+                        bandwidths,
+                        blockbandwidths,
+                        israggedbelow)
 end
 
 interlace_domainspace(ops::AbstractMatrix, ::Type{NoSpace}) = domainspace(ops)
@@ -217,7 +234,8 @@ function convert(::Type{Operator{T}},S::InterlaceOperator) where T
     else
         ops = convert(AbstractArray{Operator{T}}, S.ops)
         InterlaceOperator(ops,domainspace(S),rangespace(S),
-                            S.domaininterlacer,S.rangeinterlacer,S.bandwidths)
+                            S.domaininterlacer,S.rangeinterlacer,S.bandwidths,
+                            S.blockbandwidths, S.israggedbelow)
     end
 end
 
@@ -226,11 +244,7 @@ end
 #TODO: More efficient to save bandwidth
 bandwidths(M::InterlaceOperator) = M.bandwidths
 
-blockbandwidths(M::InterlaceOperator) =
-    (mapreduce(op->blockbandwidth(op,1),max,M.ops),
-     mapreduce(op->blockbandwidth(op,2),max,M.ops))
-
-isblockbanded(M::InterlaceOperator) = all(isblockbanded,M.ops)
+blockbandwidths(M::InterlaceOperator) = M.blockbandwidths
 
 function blockcolstop(M::InterlaceOperator,J::Integer)
     if isblockbandedbelow(M)
@@ -262,7 +276,7 @@ function colstop(M::InterlaceOperator, j::Integer)
     end
 end
 
-israggedbelow(M::InterlaceOperator) = all(israggedbelow,M.ops)
+israggedbelow(M::InterlaceOperator) = M.israggedbelow
 
 getindex(op::InterlaceOperator,k::Integer,j::Integer) =
     error("Higher tensor InterlaceOperators not supported")

src/Operators/general/algebra.jl

@@ -21,7 +21,7 @@ struct PlusOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW} <: Operator{T}
     end
 end
 
-bandwidthsmax(ops, f=bandwidths) = mapreduce(f, (t1, t2) -> max.(t1, t2), ops, init=(-720, -720)) #= approximate (-∞,-∞) =#
+bandwidthsmax(ops, f=bandwidths) = mapfoldl(f, (t1, t2) -> max.(t1, t2), ops, init=(-720, -720)) #= approximate (-∞,-∞) =#
 
 function PlusOperator(ops::Vector{O}, args...) where {O<:Operator}
     PlusOperator{eltype(O)}(ops, args...)

src/Spaces/ProductSpaceOperators.jl

@@ -162,10 +162,19 @@ function Conversion(a::SumSpace, b::Space)
     end
 
     m=zeros(Operator{promote_type(prectype(a), prectype(b))},1,length(a.spaces))
-    for n=1:length(a.spaces)
-        m[1,n]=Conversion(a.spaces[n],b)
-    end
-    return ConversionWrapper(InterlaceOperator(m, a, b, cache(interlacer(a)), cache(BlockInterlacer((Fill(1,∞),))), (1-dimension(b),dimension(a)-1)))
+    ops = map(x -> Conversion(x,b), a.spaces)
+    copyto!(m, ops)
+    bbw = bandwidthsmax(ops, blockbandwidths)
+    irb = any(israggedbelow, ops)
+    return ConversionWrapper(
+        InterlaceOperator(m, a, b,
+            cache(interlacer(a)),
+            cache(BlockInterlacer((Fill(1,∞),))),
+            (1-dimension(b),dimension(a)-1),
+            bbw,
+            irb,
+        )
+    )
 end
 
 
@@ -229,7 +238,10 @@ _spacename(::PiecewiseSpace) = PiecewiseSpace
     opbw = map(bandwidths, t)
     D = Diagonal(convert_vector_or_svector(t))
     iopbw = interlace_bandwidths(D, ds, rs, allbanded, opbw)
-    InterlaceOperator(D, ds, rs, bandwidths = iopbw)
+    iopbbw = bandwidthsmax(t, blockbandwidths)
+    irb = all(israggedbelow, t)
+    InterlaceOperator(D, ds, rs,
+        bandwidths = iopbw, blockbandwidths = iopbbw, israggedbelow = irb)
 end
 
 for (Op,OpWrap) in ((:Derivative,:DerivativeWrapper),(:Integral,:IntegralWrapper))