Move more here

Author: dlfivefifty

src/ApproxFunBase.jl

@@ -25,7 +25,7 @@ import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !,
                 getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                 @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                 Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
-                split
+                split, iszero
 
 import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle, broadcastable,
                         DefaultArrayStyle, broadcasted

src/LinearAlgebra/helper.jl

@@ -620,11 +620,6 @@ end
 
 length(A::CachedIterator) = length(A.iterator)
 
-
-# The following don't need caching
-cache(A::AbstractVector{T}) where {T<:Number} = A
-
-
 ## nocat
 vnocat(A...) = Base.vect(A...)
 hnocat(A...) = Base.typed_hcat(mapreduce(typeof,promote_type,A),A...)

src/Operators/general/PartialInverseOperator.jl

@@ -37,3 +37,33 @@ function getindex(P::PartialInverseOperator,k::Integer,j::Integer)
         return zero(eltype(P))
     end
 end
+
+
+## These are both hacks that apparently work
+
+function BandedMatrix(S::SubOperator{T,PP,Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,PP<:PartialInverseOperator}
+    kr,jr = parentindices(S)
+    P = parent(S)
+    #ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
+    ret = BandedMatrix{eltype(S)}(undef, (last(kr),last(jr)), bandwidths(P))
+    b = bandwidth(P, 2)
+    #@assert first(kr) == first(jr) == 1
+
+    @inbounds for j in 1:last(jr)
+        kk = colrange(ret, j)
+        if j in kk
+            ret[j,j] = inv(P.cache[j,j])
+        end
+        for k in first(kk):min(last(kk),j-1)
+            t = zero(T)
+            for i = max(k,j-b-1):j-1
+                t += ret[k,i]*P.cache[i,j]
+            end
+            ret[k,j] = -t/P.cache[j,j]
+        end
+    end
+
+    ret[kr,jr]
+end
+
+

src/Spaces/QuotientSpace.jl

@@ -260,3 +260,40 @@ for (gesdd, elty, relty) in ((:dgesdd_,:Float64,:Float64),
         end
     end
 end
+
+
+function BandedMatrix(S::SubOperator{T,ConcreteConversion{QuotientSpace{SP,O,D,R},SP,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {SP,O,D,R,T}
+    kr,jr = parentindices(S)
+    C = parent(S)
+    #ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
+    ret = BandedMatrix{eltype(S)}(undef, (last(kr),last(jr)), bandwidths(C))
+    #@assert first(kr) == first(jr) == 1
+
+    sp = domainspace(C)
+    F = sp.F
+    A = F.factors
+    n = size(A, 1)
+    B = sp.bcs[1:n,1:last(jr)+n]
+    x = sp.x
+    @inbounds for j in 1:last(jr)
+        kk = colrange(ret, j)
+        if j in kk
+            ret[j,j] = one(R)
+        end
+        for jj = 1:n, ii = 1:n
+            A[ii,jj] = B[ii,j+jj]
+        end
+        for ii = 1:n
+            x[ii] = -B[ii,j]
+        end
+        if norm(x) > 8*norm(A)*eps(R)
+            mutable_lu!(F)
+            ldiv!(F, x)
+        end
+        for k = first(kk)+1:last(kk)
+            ret[k,j] = x[k-j]
+        end
+    end
+
+    ret[kr,jr]
+end

src/Spaces/Spaces.jl

@@ -140,10 +140,10 @@ abstract type WeightSpace{S,DD,RR} <: Space{DD,RR} end
 
 
 domain(S::WeightSpace) = domain(S.space)
-
-
 points(sp::WeightSpace,n) = points(sp.space,n)
 
+weight(S::WeightSpace, x...) = error("Override `weight(::$S, $x)`")
+
 
 struct WeightSpacePlan{S,P,T,V}
     space::S