Refinements

Author: dlfivefifty

src/ApproxFunBase.jl

@@ -4,7 +4,7 @@ module ApproxFunBase
     using Base, BlockArrays, BandedMatrices, BlockBandedMatrices, DomainSets, IntervalSets,
             SpecialFunctions, AbstractFFTs, FFTW, SpecialFunctions, DSP, DualNumbers, FastTransforms,
             LinearAlgebra, SparseArrays, LowRankApprox, FillArrays, InfiniteArrays #, Arpack
-import StaticArrays
+import StaticArrays, Calculus
 
 import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
             Interval, ChebyshevInterval, boundary, ∂, rightendpoint, leftendpoint,

src/Fun.jl

@@ -20,6 +20,10 @@ 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)
+
 
 function Fun(sp::Space,v::AbstractVector{Any})
     if isempty(v)

src/Operators/Operator.jl

@@ -486,11 +486,11 @@ macro wrappergetindex(Wrap)
         BLAS.axpy!(α,P::ApproxFunBase.SubOperator{T,OP},A::AbstractMatrix) where {T,OP<:$Wrap} =
             ApproxFunBase.unwrap_axpy!(α,P,A)
 
-        ApproxFunBase.mul_coefficients(A::$Wrap,b) = mul_coefficients(A.op,b)
+        ApproxFunBase.mul_coefficients(A::$Wrap,b) = ApproxFunBase.mul_coefficients(A.op,b)
         ApproxFunBase.mul_coefficients(A::ApproxFunBase.SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},b) where {T,OP<:$Wrap} =
-            mul_coefficients(view(parent(A).op,S.indexes[1],S.indexes[2]),b)
+            ApproxFunBase.mul_coefficients(view(parent(A).op,S.indexes[1],S.indexes[2]),b)
         ApproxFunBase.mul_coefficients(A::ApproxFunBase.SubOperator{T,OP},b) where {T,OP<:$Wrap} =
-            mul_coefficients(view(parent(A).op,S.indexes[1],S.indexes[2]),b)
+            ApproxFunBase.mul_coefficients(view(parent(A).op,S.indexes[1],S.indexes[2]),b)
     end
 
     for TYP in (:(BandedMatrices.BandedMatrix),:(ApproxFunBase.RaggedMatrix),
@@ -779,3 +779,64 @@ end
 
 AbstractMatrix(V::Operator) = arraytype(V)(V)
 AbstractVector(S::Operator) = Vector(S)
+
+
+
+
+# default copy is to loop through
+# override this for most operators.
+function default_BandedMatrix(S::Operator)
+    Y=BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
+
+    for j=1:size(S,2),k=colrange(Y,j)
+        @inbounds inbands_setindex!(Y,S[k,j],k,j)
+    end
+
+    Y
+end
+
+
+# default copy is to loop through
+# override this for most operators.
+function default_RaggedMatrix(S::Operator)
+    data=Array{eltype(S)}(undef, 0)
+    cols=Array{Int}(undef, size(S,2)+1)
+    cols[1]=1
+    for j=1:size(S,2)
+        cs=colstop(S,j)
+        K=cols[j]-1
+        cols[j+1]=cs+cols[j]
+        resize!(data,cols[j+1]-1)
+
+        for k=1:cs
+            data[K+k]=S[k,j]
+        end
+    end
+
+    RaggedMatrix(data,cols,size(S,1))
+end
+
+function default_Matrix(S::Operator)
+    n, m = size(S)
+    if isinf(n) || isinf(m)
+        error("Cannot convert $S to a Matrix")
+    end
+
+    eltype(S)[S[k,j] for k=1:n, j=1:m]
+end
+
+
+
+
+# The diagonal of the operator may not be the diagonal of the sub
+# banded matrix, so the following calculates the row of the
+# Banded matrix corresponding to the diagonal of the original operator
+
+
+diagindshift(S,kr,jr) = first(kr)-first(jr)
+diagindshift(S::SubOperator) = diagindshift(S,parentindices(S)[1],parentindices(S)[2])
+
+
+#TODO: Remove
+diagindrow(S,kr,jr) = bandwidth(S,2)+first(jr)-first(kr)+1
+diagindrow(S::SubOperator) = diagindrow(S,parentindices(S)[1],parentindices(S)[2])

src/Operators/banded/ToeplitzOperator.jl

@@ -170,3 +170,112 @@ function Fun(T::ToeplitzOperator)
         Fun(Laurent(Circle()),interlace(T.nonnegative,T.negative))
     end
 end
+
+
+
+####
+# Toeplitz/Hankel
+####
+
+
+# αn,α0,αp give the constant for the negative, diagonal and positive
+# entries.  The usual case is 1,2,1
+function toeplitz_axpy!(αn,α0,αp,neg,pos,kr,jr,ret)
+    dat=ret.data
+    m=size(dat,2)
+
+    dg=diagindrow(ret,kr,jr)
+
+    # diagonal
+    if dg ≥ 1
+        α0p=α0*pos[1]
+        @simd for j=1:m
+            @inbounds dat[dg,j]+=α0p
+        end
+    end
+
+    # positive entries
+    for k=2:min(length(pos),dg)
+        αpp=αp*pos[k]
+        @simd for j=1:m
+            @inbounds dat[dg-k+1,j]+=αpp
+        end
+    end
+
+    # negative entries
+    for k=1:min(length(neg),size(dat,1)-dg)
+        αnn=αn*neg[k]
+        @simd for j=1:m
+            @inbounds dat[dg+k,j]+=αnn
+        end
+    end
+
+    ret
+end
+
+# this routine is for when we want a symmetric toeplitz
+function sym_toeplitz_axpy!(αn,α0,αp,cfs,kr,jr,ret)
+    dg=diagindrow(ret,kr,jr)
+
+    dat=ret.data
+    m=size(dat,2)
+    # diagonal
+    if dg ≥ 1
+        α0p=α0*cfs[1]
+        @simd for j=1:m
+            @inbounds dat[dg,j]+=α0p
+        end
+    end
+
+    # positive entries
+    for k=2:min(length(cfs),dg)
+        αpp=αp*cfs[k]
+        @simd for j=1:m
+            @inbounds dat[dg-k+1,j]+=αpp
+        end
+    end
+
+    # negative entries
+    for k=2:min(length(cfs),size(dat,1)-dg+1)
+        αnn=αn*cfs[k]
+        @simd for j=1:m
+            @inbounds dat[dg+k-1,j]+=αnn
+        end
+    end
+
+    ret
+end
+
+toeplitz_axpy!(α,neg,pos,kr,jr,ret) =
+    toeplitz_axpy!(α,α,α,neg,pos,kr,jr,ret)
+
+sym_toeplitz_axpy!(α0,αp,cfs,kr,jr,ret) =
+    sym_toeplitz_axpy!(αp,α0,αp,cfs,kr,jr,ret)
+
+
+function hankel_axpy!(α,cfs,kr,jr,ret)
+    dat=ret.data
+
+    st=stride(dat,2)-2
+    mink=first(kr)+first(jr)-1
+
+    N=length(dat)
+
+    # dg gives the row corresponding to the diagonal of the original operator
+    dg=diagindrow(ret,kr,jr)
+
+    # we need the entry where the first entry is written to
+    # this is going to be a shift of the diagonal of the true operator
+    dg1=dg+first(kr)-first(jr)
+
+    for k=mink:min(length(cfs),size(dat,1)+mink-dg1)
+        dk=k-mink+dg1
+        nk=k-mink
+        αc=α*cfs[k]
+        @simd for j=dk:st:min(dk+nk*st,N)
+            @inbounds dat[j]+=αc
+        end
+    end
+
+    ret
+end

src/Spaces/ConstantSpace.jl

@@ -240,3 +240,69 @@ choosedomainspace(M::ConcreteMultiplication{D,UnsetSpace},sp::UnsetSpace) where
 choosedomainspace(M::ConcreteMultiplication{D,UnsetSpace},sp::Space) where {D<:ConstantSpace} = space(M.f)
 
 Base.isfinite(f::Fun{CS}) where {CS<:ConstantSpace} = isfinite(Number(f))
+
+
+
+
+
+## ConstantSpace and PointSpace default overrides
+
+for SP in (:ConstantSpace,)
+    for OP in (:abs,:sign,:exp,:sqrt,:angle)
+        @eval begin
+            $OP(z::Fun{<:$SP,<:Complex}) = Fun(space(z),$OP.(coefficients(z)))
+            $OP(z::Fun{<:$SP,<:Real}) = Fun(space(z),$OP.(coefficients(z)))
+            $OP(z::Fun{<:$SP}) = Fun(space(z),$OP.(coefficients(z)))
+        end
+    end
+
+    # we need to pad coefficients since 0^0 == 1
+    for OP in (:^,)
+        @eval begin
+            function $OP(z::Fun{<:$SP},k::Integer)
+                k ≠ 0 && return Fun(space(z),$OP.(coefficients(z),k))
+                Fun(space(z),$OP.(pad(coefficients(z),dimension(space(z))),k))
+            end
+            function $OP(z::Fun{<:$SP},k::Number)
+                k ≠ 0 && return Fun(space(z),$OP.(coefficients(z),k))
+                Fun(space(z),$OP.(pad(coefficients(z),dimension(space(z))),k))
+            end
+        end
+    end
+end
+
+
+for OP in (:(Base.max),:(Base.min))
+    @eval begin
+        $OP(a::Fun{CS1,T},b::Fun{CS2,V}) where {CS1<:ConstantSpace,CS2<:ConstantSpace,T<:Real,V<:Real} =
+            Fun($OP(Number(a),Number(b)),space(a) ∪ space(b))
+        $OP(a::Fun{CS,T},b::Real) where {CS<:ConstantSpace,T<:Real} =
+            Fun($OP(Number(a),b),space(a))
+        $OP(a::Real,b::Fun{CS,T}) where {CS<:ConstantSpace,T<:Real} =
+            Fun($OP(a,Number(b)),space(b))
+    end
+end
+
+for OP in (:<,:(Base.isless),:(<=),:>,:(>=))
+    @eval begin
+        $OP(a::Fun{CS},b::Fun{CS}) where {CS<:ConstantSpace} = $OP(convert(Number,a),Number(b))
+        $OP(a::Fun{CS},b::Number) where {CS<:ConstantSpace} = $OP(convert(Number,a),b)
+        $OP(a::Number,b::Fun{CS}) where {CS<:ConstantSpace} = $OP(a,convert(Number,b))
+    end
+end
+
+# from DualNumbers
+for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
+    funsym == :abs && continue
+    funsym == :sign && continue
+    funsym == :exp && continue
+    funsym == :sqrt && continue
+    @eval begin
+        $(funsym)(z::Fun{CS,T}) where {CS<:ConstantSpace,T<:Real} =
+            Fun($(funsym)(Number(z)),space(z))
+        $(funsym)(z::Fun{CS,T}) where {CS<:ConstantSpace,T<:Complex} =
+            Fun($(funsym)(Number(z)),space(z))
+        $(funsym)(z::Fun{CS}) where {CS<:ConstantSpace} =
+            Fun($(funsym)(Number(z)),space(z))
+    end
+end