Out-of-place pad for BandedMatrix in resizedata (#430)

Author: jishnub

src/Caching/almostbanded.jl

@@ -181,7 +181,7 @@ function resizedata!(co::CachedOperator{T,<:AlmostBandedMatrix{T},<:InterlaceOpe
     end
 
     (l,u)=bandwidths(co.data.bands)
-    pad!(co.data,n,n+u)
+    co.data = pad(co.data, n, n+u)
 
     r = rank(co.data.fill)
     ind = findfirst(op->isinf(size(op,1)),co.op.ops)
@@ -227,8 +227,7 @@ function resizedata!(co::CachedOperator{T,<:AlmostBandedMatrix{T},<:InterlaceOpe
     p = length(d∞)
 
     (l,u)=bandwidths(co.data.bands)
-    pad!(co.data,n,n+u)
-    co.data
+    co.data = pad(co.data,n,n+u)
     # r is number of extra rows, ncols is number of extra columns
     r = rank(co.data.fill)
     ncols = mapreduce(d->isfinite(d) ? d : 0,+,ddims)
@@ -286,36 +285,38 @@ function resizedata!(QR::QROperator{<:CachedOperator{T,<:AlmostBandedMatrix{T}}}
     MO = QR.R_cache
     W = QR.H
 
-    R = MO.data.bands
-    M = R.l+1   # number of diag+subdiagonal bands
+    Rl, Ru = bandwidths(MO.data.bands)
+    M = Rl + 1   # number of diag+subdiagonal bands
 
     if col+M-1 ≥ MO.datasize[1]
         resizedata!(MO,(col+M-1)+100,:)  # double the last rows
     end
 
+    R = MO.data.bands # has to be accessed after the resizedata!
+
     if col > size(W,2)
         W = QR.H = unsafe_resize!(W,:,2col)
     end
 
     F = MO.data.fill.U
 
     for k = QR.ncols+1:col
-        W[:,k] = view(R.data, (R.u+1).+(0:R.l), k) # diagonal and below
+        W[:,k] = view(R.data, (Ru+1).+(0:Rl), k) # diagonal and below
         wp = view(W,:,k)
         W[1,k]+= flipsign(norm(wp),W[1,k])
         normalize!(wp)
 
         # scale banded entries
-        for j = k:k+R.u
-            dind = R.u+1+k-j
+        for j = k:k+Ru
+            dind = Ru+1+k-j
             v = view(R.data, range(dind, length=M), j)
             dt = dot(wp,v)
             axpy!(-2dt,wp,v)
         end
 
         # scale banded/filled entries
-        for j = (k+R.u).+(1:M-1)
-            p = j-k-R.u
+        for j = (k+Ru).+(1:M-1)
+            p = j-k-Ru
             v = view(R.data,1:M-p,j)  # shift down each time
             wp2=view(wp,p+1:M)
             dt = dot(wp2,v)

src/Caching/banded.jl

@@ -17,7 +17,7 @@ function resizedata!(B::CachedOperator{T,<:BandedMatrix{T}},n::Integer,m_in::Int
     n = min(n, N)
 
     if n > B.datasize[1]
-        pad!(B.data,min(N,2n),m)
+        B.data = pad(B.data,min(N,2n),m)
 
         kr=B.datasize[1]+1:n
         jr=max(B.datasize[1]+1-B.data.l,1):min(n+B.data.u,M)
@@ -52,34 +52,36 @@ function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T}}}, ::Co
     MO=QR.R_cache
     W=QR.H
 
-    R=MO.data
-    M=R.l+1   # number of diag+subdiagonal bands
+    Rl,Ru = bandwidths(MO.data)
+    M = Rl + 1   # number of diag+subdiagonal bands
 
     if col+M-1 ≥ MO.datasize[1]
         resizedata!(MO,(col+M-1)+100,:)  # double the last rows
     end
 
+    R = MO.data # has to be accessed after resizedata!, as the matrix might change
+
     if col > size(W,2)
         W=QR.H=unsafe_resize!(W,:,2col)
     end
 
     for k=QR.ncols+1:col
-        W[:,k] = view(R.data, (R.u+1).+(0:R.l), k) # diagonal and below
+        W[:,k] = view(R.data, (Ru+1).+(0:Rl), k) # diagonal and below
         wp=view(W,:,k)
         W[1,k]+= flipsign(norm(wp),W[1,k])
         normalize!(wp)
 
         # scale banded entries
-        for j=k:k+R.u
-            dind=R.u+1+k-j
+        for j=k:k+Ru
+            dind=Ru+1+k-j
             v=view(R.data, range(dind, length=M), j)
             dt=dot(wp,v)
             axpy!(-2dt,wp,v)
         end
 
         # scale banded/filled entries
-        for j = (k+R.u).+(1:M-1)
-            p=j-k-R.u
+        for j = (k+Ru).+(1:M-1)
+            p=j-k-Ru
             v=view(R.data,1:M-p,j)  # shift down each time
             wp2=view(wp,p+1:M)
             dt=dot(wp2,v)

src/Fun.jl

@@ -86,6 +86,7 @@ function Fun(sp::Space,v::AbstractVector{Any})
     end
 end
 
+Fun(f::Fun) = f # Fun of Fun should be like a conversion
 
 hasnumargs(f::Fun,k) = k == 1 || domaindimension(f) == k  # all funs take a single argument as a SVector
 

src/LinearAlgebra/AlmostBandedMatrix.jl

@@ -49,8 +49,8 @@ function setindex!(B::AlmostBandedMatrix,v,k::Integer,j::Integer)
 end
 
 
-function pad!(B::AlmostBandedMatrix,n::Integer,m::Integer)
-    pad!(B.bands,n,m)
-    pad!(B.fill,n,m)
-    B
+function pad(B::AlmostBandedMatrix,n::Integer,m::Integer)
+    bands = pad(B.bands,n,m)
+    fill = pad(B.fill,n,m)
+    AlmostBandedMatrix(bands, fill)
 end

src/LinearAlgebra/LowRankMatrix.jl

@@ -34,12 +34,12 @@ end
 
 # constructors
 
-function pad!(L::LowRankMatrix,n::Integer,::Colon)
-    L.U=pad(L.U,n,:)
-    L
+function pad(L::LowRankMatrix,n::Integer,::Colon)
+    U = pad(L.U, n, :)
+    LowRankMatrix(U, L.V)
 end
-function pad!(L::LowRankMatrix,::Colon,m::Integer)
-    L.V=pad(L.V,m,:)
-    L
+function pad(L::LowRankMatrix,::Colon,m::Integer)
+    V = pad(L.V,m,:)
+    LowRankMatrix(L.U, V)
 end
-pad!(L::LowRankMatrix,n::Integer,m::Integer) = pad!(pad!(L,n,:),:,m)
+pad(L::LowRankMatrix,n::Integer,m::Integer) = pad(pad(L,n,:),:,m)

src/LinearAlgebra/helper.jl

@@ -260,6 +260,13 @@ function pad(A::AbstractMatrix,n::Integer,m::Integer)
 	end
 end
 
+function pad(A::BandedMatrix, n::Integer, m::Integer)
+    B = BandedMatrix{eltype(A)}(undef, (n,m), bandwidths(A))
+    copyto!(B.data, A.data)
+    B.data[length(A.data)+1:end] .= 0
+    return B
+end
+
 pad(A::AbstractMatrix,::Colon,m::Integer) = pad(A,size(A,1),m)
 pad(A::AbstractMatrix,n::Integer,::Colon) = pad(A,n,size(A,2))
 
@@ -534,7 +541,7 @@ Base.isless(x::PosInfinity, y::Block{1}) = isless(x, Int(y))
 
 
 
-pad!(A::BandedMatrix,n,::Colon) = pad!(A,n,n+A.u)  # Default is to get all columns
+pad(A::BandedMatrix,n,::Colon) = pad(A,n,n+A.u)  # Default is to get all columns
 columnrange(A,row::Integer) = max(1,row-bandwidth(A,1)):row+bandwidth(A,2)
 
 

src/constructors.jl

@@ -242,3 +242,47 @@ Fun(f,n::Integer) = Fun(f,ChebyshevInterval(),n)
 Fun(T::Type,d::AbstractVector) = Fun(T(),d)
 
 Fun(f::Fun{SequenceSpace},s::Space) = Fun(s,f.coefficients)
+
+"""
+    Fun(f)
+
+Return `Fun(f, space)` by choosing an appropriate `space` for the function.
+For univariate functions, `space` is chosen to be `Chebyshev()`, whereas for
+multivariate functions, it is a tensor product of `Chebyshev()` spaces.
+
+# Examples
+```jldoctest
+julia> f = Fun(x -> x^2)
+Fun(Chebyshev(), [0.5, 0.0, 0.5])
+
+julia> f(0.1) == (0.1)^2
+true
+
+julia> f = Fun((x,y) -> x + y);
+
+julia> f(0.1, 0.2) ≈ 0.3
+true
+```
+"""
+function Fun(f::Function)
+    if hasonearg(f)
+        # check for tuple
+        try
+            f(0)
+        catch ex
+            if ex isa BoundsError
+                # assume its a tuple
+                return Fun(f,ChebyshevInterval()^2)
+            else
+                rethrow()
+            end
+        end
+
+        Fun(f,ChebyshevInterval())
+    elseif hasnumargs(f,2)
+            Fun(f,ChebyshevInterval()^2)
+    else
+        error("Function not defined on interval or square")
+    end
+end
+

src/hacks.jl

@@ -1,66 +1,3 @@
-## Functions that depend on the structure of BandedMatrix
-
-
-function pad!(A::BandedMatrix,n,m)
-    A.data = pad(A.data,size(A.data,1),m)
-    A.raxis = Base.OneTo(n)
-    A
-end
-
-
-
-# linear algebra
-
-
-## Constructors that involve MultivariateFun
-Fun(f::Fun) = f # Fun of Fun should be like a conversion
-
-"""
-    Fun(f)
-
-Return `Fun(f, space)` by choosing an appropriate `space` for the function.
-For univariate functions, `space` is chosen to be `Chebyshev()`, whereas for
-multivariate functions, it is a tensor product of `Chebyshev()` spaces.
-
-# Examples
-```jldoctest
-julia> f = Fun(x -> x^2)
-Fun(Chebyshev(), [0.5, 0.0, 0.5])
-
-julia> f(0.1) == (0.1)^2
-true
-
-julia> f = Fun((x,y) -> x + y);
-
-julia> f(0.1, 0.2) ≈ 0.3
-true
-```
-"""
-function Fun(f::Function)
-    if hasonearg(f)
-        # check for tuple
-        try
-            f(0)
-        catch ex
-            if isa(ex,BoundsError)
-                # assume its a tuple
-                return Fun(f,ChebyshevInterval()^2)
-            else
-                throw(ex)
-            end
-        end
-
-        Fun(f,ChebyshevInterval())
-    elseif hasnumargs(f,2)
-            Fun(f,ChebyshevInterval()^2)
-    else
-        error("Function not defined on interval or square")
-    end
-end
-
-
-
-
 ## These hacks support PDEs with block matrices