Skip to content

Commit efe368f

Browse files
jishnubjishnub
authored
despecialize union, split default_mult (#150)
* despecialize union * missing return in union * inference in Fun constructor * revert coefficient coercion
1 parent 00b3e94 commit efe368f

File tree

7 files changed

+114
-84
lines changed

7 files changed

+114
-84
lines changed

‎src/Fun.jl‎

Lines changed: 9 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -43,10 +43,11 @@ hasnumargs(f::Fun,k) = k == 1 || domaindimension(f) == k # all funs take a sing
4343

4444
function coefficients(f::Fun,msp::Space)
4545
#zero can always be converted
46-
if ncoefficients(f) == 0 || (ncoefficients(f) == 1 && f.coefficients[1] == 0)
47-
f.coefficients
46+
fc = f.coefficients
47+
if ncoefficients(f) == 0 || (ncoefficients(f) == 1 && fc[1] == 0)
48+
fc
4849
else
49-
coefficients(f.coefficients, space(f), msp)
50+
coefficients(fc, space(f), msp)
5051
end
5152
end
5253
coefficients(f::Fun,::Type{T}) where {T<:Space} = coefficients(f,T(domain(f)))
@@ -65,12 +66,7 @@ end
6566

6667
function coefficient(f::Fun,kr::AbstractRange)
6768
b = maximum(kr)
68-
69-
if b ≤ ncoefficients(f)
70-
f.coefficients[kr]
71-
else
72-
[coefficient(f,k) for k in kr]
73-
end
69+
f.coefficients[first(kr):min(b, end)]
7470
end
7571

7672
coefficient(f::Fun,K::Block) = coefficient(f,blockrange(space(f),K.n[1]))
@@ -241,9 +237,11 @@ extrapolate(f::Fun,x,y,z...) = extrapolate(f.coefficients,f.space,Vec(x,y,z...))
241237
##Data routines
242238

243239

244-
values(f::Fun,dat...) = itransform(f.space,f.coefficients,dat...)
240+
values(f::Fun,dat...) = _values(f.space, f.coefficients, dat...)
241+
_values(sp, v, dat...) = itransform(sp, v, dat...)
242+
_values(sp, v::Vector{T}, dat...) where {T} = itransform(sp, v, dat...)::Vector{float(T)}
245243
points(f::Fun) = points(f.space,ncoefficients(f))
246-
ncoefficients(f::Fun) = length(f.coefficients)
244+
ncoefficients(f::Fun)::Int = length(f.coefficients)
247245
blocksize(f::Fun) = (block(space(f),ncoefficients(f)).n[1],)
248246

249247
function stride(f::Fun)

‎src/Operators/SubOperator.jl‎

Lines changed: 19 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -48,12 +48,12 @@ end
4848
# work around strange bug with bool size
4949
SubOperator(A,inds,dims::Tuple{Bool,Bool},lu) = SubOperator(A,inds,Int.(dims),lu)
5050

51-
function SubOperator(A,inds::Tuple{Block,Block},lu)
51+
function SubOperator(A,inds::NTuple{2,Block},lu)
5252
checkbounds(A,inds...)
5353
SubOperator(A,inds,(blocklengths(rangespace(A))[inds[1].n[1]],blocklengths(domainspace(A))[inds[2].n[1]]),lu)
5454
end
5555

56-
SubOperator(A, inds::Tuple{Block,Block}) = SubOperator(A,inds,subblockbandwidths(A))
56+
SubOperator(A, inds::NTuple{2,Block}) = SubOperator(A,inds,subblockbandwidths(A))
5757
function SubOperator(A, inds::Tuple{BlockRange{1,R},BlockRange{1,R}}) where R
5858
checkbounds(A,inds...)
5959
dims = (sum(blocklengths(rangespace(A))[inds[1].indices[1]]),
@@ -169,7 +169,7 @@ view(V::SubOperator, kr, jr) = view(V.parent,reindex(V,parentindices(V),(kr,jr))
169169
view(V::SubOperator,kr::InfRanges,jr::InfRanges) = view(V.parent,reindex(V,parentindices(V),(kr,jr))...)
170170

171171
bandwidths(S::SubOperator) = S.bandwidths
172-
function colstop(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP}
172+
function colstop(S::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}},j::Integer)
173173
cs = colstop(parent(S),parentindices(S)[2][j])
174174
kr = parentindices(S)[1]
175175
n = size(S,1)
@@ -181,34 +181,34 @@ function colstop(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::In
181181
min(n,findfirst(isequal(cs),kr))
182182
end
183183
end
184-
colstart(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP} =
184+
colstart(S::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}},j::Integer) =
185185
max(findfirst(parentindices(S)[1],colstart(parent(S),parentindices(S)[2][j])),1)
186-
rowstart(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP} =
186+
rowstart(S::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}},j::Integer) =
187187
max(1,findfirst(parentindices(S)[2],rowstart(parent(S),parentindices(S)[1][j])))
188-
rowstop(S::SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},j::Integer) where {T,OP} =
188+
rowstop(S::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}},j::Integer) =
189189
findfirst(parentindices(S)[2],rowstop(parent(S),parentindices(S)[1][j]))
190190

191191

192192
# blocks don't change
193-
blockcolstop(S::SubOperator{T,OP,Tuple{II,JJ}},J::Integer) where {T,OP,II<:AbstractRange{Int},JJ<:AbstractRange{Int}} =
193+
blockcolstop(S::SubOperator{<:Any,<:Any,Tuple{AbstractRange{Int},AbstractRange{Int}}},J::Integer) =
194194
blockcolstop(parent(S),J)
195195

196196
israggedbelow(S::SubOperator) = israggedbelow(parent(S))
197197

198198
# since blocks don't change with indexex, neither do blockbandwidths
199-
blockbandwidths(S::SubOperator{T,OP,Tuple{II,JJ}}) where {T,OP,II<:AbstractRange{Int},JJ<:AbstractRange{Int}} =
199+
blockbandwidths(S::SubOperator{<:Any,<:Any,NTuple{2,AbstractRange{Int}}}) =
200200
blockbandwidths(parent(S))
201-
function blockbandwidths(S::SubOperator{T,B,Tuple{BlockRange1,BlockRange1}}) where {T,B}
201+
function blockbandwidths(S::SubOperator{<:Any,<:Any,NTuple{2,BlockRange1}})
202202
KR,JR = parentindices(S)
203203
l,u = blockbandwidths(parent(S))
204204
sh = first(KR).n[1]-first(JR).n[1]
205205
l-sh,u+sh
206206
end
207207

208-
isblockbanded(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = false
209-
isbanded(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = isbandedblockbanded(parent(S))
210-
bandwidths(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = subblockbandwidths(parent(S))
211-
blockbandwidths(S::SubOperator{T,B,Tuple{Block,Block}}) where {T,B} = 0,0
208+
isblockbanded(S::SubOperator{<:Any,<:Any,NTuple{2,Block}}) = false
209+
isbanded(S::SubOperator{<:Any,<:Any,NTuple{2,Block}}) = isbandedblockbanded(parent(S))
210+
bandwidths(S::SubOperator{<:Any,<:Any,NTuple{2,Block}}) = subblockbandwidths(parent(S))
211+
blockbandwidths(S::SubOperator{<:Any,<:Any,NTuple{2,Block}}) = 0,0
212212

213213
function BandedBlockBandedMatrix(::Type{Zeros}, S::SubOperator)
214214
kr,jr=parentindices(S)
@@ -234,7 +234,7 @@ function BandedBlockBandedMatrix(::Type{Zeros}, S::SubOperator)
234234
rows,cols, (l,u), (λ-jsh,μ+ksh))
235235
end
236236

237-
function BandedBlockBandedMatrix(::Type{Zeros}, S::SubOperator{T,B,Tuple{BlockRange1,BlockRange1}}) where {T,B}
237+
function BandedBlockBandedMatrix(::Type{Zeros}, S::SubOperator{<:Any,<:Any,Tuple{BlockRange1,BlockRange1}})
238238
KR,JR = parentindices(S)
239239
KO = parent(S)
240240
l,u = blockbandwidths(KO)::Tuple{Int,Int}
@@ -297,7 +297,7 @@ end
297297
_colstops(V) = Int[max(0,colstop(V,j)) for j=1:size(V,2)]
298298

299299
for TYP in (:RaggedMatrix, :Matrix)
300-
def_TYP = Meta.parse("default_" * string(TYP))
300+
def_TYP = Symbol(:default_, TYP)
301301
@eval begin
302302
function $TYP(V::SubOperator)
303303
if isinf(size(V,1)) || isinf(size(V,2))
@@ -311,7 +311,7 @@ for TYP in (:RaggedMatrix, :Matrix)
311311
end
312312
end
313313

314-
function $TYP(V::SubOperator{T,BB,NTuple{2,UnitRange{Int}}}) where {T,BB}
314+
function $TYP(V::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}})
315315
if isinf(size(V,1)) || isinf(size(V,2))
316316
error("Cannot convert $V to a $TYP")
317317
end
@@ -331,7 +331,7 @@ for TYP in (:RaggedMatrix, :Matrix)
331331
end
332332

333333
# fast converts to banded matrices would be based on indices, not blocks
334-
function BandedMatrix(S::SubOperator{T,B,Tuple{BlockRange1,BlockRange1}}) where {T,B}
334+
function BandedMatrix(S::SubOperator{<:Any,<:Any,Tuple{BlockRange1,BlockRange1}})
335335
A = parent(S)
336336
ds = domainspace(A)
337337
rs = rangespace(A)
@@ -344,14 +344,14 @@ end
344344

345345

346346

347-
function mul_coefficients(A::SubOperator{T,B,Tuple{UnitRange{Int},UnitRange{Int}}},b) where {T,B}
347+
function mul_coefficients(A::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}}, b)
348348
if size(A,2) == length(b)
349349
AbstractMatrix(A)*b
350350
else
351351
view(A,:,1:length(b))*b
352352
end
353353
end
354-
function mul_coefficients!(A::SubOperator{T,B,Tuple{UnitRange{Int},UnitRange{Int}}},b) where {T,B}
354+
function mul_coefficients!(A::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}}, b)
355355
if size(A,2) == length(b)
356356
mul!(b, AbstractMatrix(A), b)
357357
else

‎src/Operators/banded/Multiplication.jl‎

Lines changed: 26 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -114,43 +114,48 @@ end
114114

115115

116116
hasfasttransform(_) = false
117-
hasfasttransform(f::Fun) = hasfasttransform(space(f))
118-
hasfasttransformtimes(f,g) = pointscompatible(f,g) && hasfasttransform(f) && hasfasttransform(g)
117+
hasfasttransform(f::Fun)::Bool = hasfasttransform(space(f))
118+
hasfasttransformtimes(f,g)::Bool = pointscompatible(f,g) && hasfasttransform(f) && hasfasttransform(g)
119119

120120

121+
function default_mult_compatible(f::Fun, g::Fun)
122+
m,n = ncoefficients(f), ncoefficients(g)
123+
# Heuristic division of parameter space between value-space and coefficient-space multiplication.
124+
if log10(m)*log10(n)>4 && hasfasttransformtimes(f,g)
125+
transformtimes(f,g)
126+
elseif m≤n
127+
coefficienttimes(f,g)
128+
else
129+
coefficienttimes(g,f)
130+
end
131+
end
121132
# This should be overriden whenever the multiplication space is different
122133
function default_mult(f::Fun,g::Fun)
123134
# When the spaces differ we promote and multiply
124135
if domainscompatible(space(f),space(g))
125-
m,n = ncoefficients(f),ncoefficients(g)
126-
# Heuristic division of parameter space between value-space and coefficient-space multiplication.
127-
if hasfasttransformtimes(f,g) && log10(m)*log10(n)>4
128-
transformtimes(f,g)
129-
elseif m≤n
130-
coefficienttimes(f,g)
131-
else
132-
coefficienttimes(g,f)
133-
end
136+
default_mult_compatible(f, g)
134137
else
135-
sp=union(space(f),space(g))
136-
Fun(f,sp)*Fun(g,sp)
138+
sp=union(space(f),space(g))::Space
139+
fnew = Fun(f,sp)
140+
gnew = Fun(g,sp)
141+
default_mult_compatible(fnew, gnew)
137142
end
138143
end
139144

140145
*(f::Fun,g::Fun) = default_mult(f,g)
141146

142147
coefficienttimes(f::Fun,g::Fun) = Multiplication(f,space(g))*g
143148

144-
function transformtimes(f::Fun,g::Fun,n)
145-
@assert pointscompatible(space(f),space(g))
146-
isempty(f.coefficients) && return f
147-
isempty(g.coefficients) && return g
148-
f2,g2,sp = pad(f,n),pad(g,n),space(f)
149-
hc = transform(sp,values(f2).*values(g2))
149+
function transformtimes(f::Fun,g::Fun, n = ncoefficients(f) + ncoefficients(g) - 1, sp = space(f))
150+
@assert pointscompatible(sp,space(g))::Bool
151+
iszero(ncoefficients(f)) && return f
152+
iszero(ncoefficients(g)) && return g
153+
f2,g2 = pad(f,n), pad(g,n)
154+
v = values(f2)
155+
v .*= values(g2)
156+
hc = transform(sp, v)
150157
chop!(Fun(sp,hc),10eps(eltype(hc)))
151158
end
152-
transformtimes(f::Fun,g::Fun) = transformtimes(f,g,ncoefficients(f) + ncoefficients(g) - 1)
153-
154159

155160

156161
*(a::Fun,L::UniformScaling) = Multiplication(a*L.λ,UnsetSpace())

‎src/Operators/qr.jl‎

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -146,27 +146,27 @@ mul_coefficients(At::Transpose{T,<:QROperatorQ{T}},B::AbstractMatrix{T}) where {
146146
mul_coefficients!(At::Transpose{T,<:QROperatorQ{T}},B::AbstractVector{T}) where {T<:Real} = mul!(B, parent(At)', B)
147147
mul_coefficients!(At::Transpose{T,<:QROperatorQ{T}},B::AbstractMatrix{T}) where {T<:Real} = mul!(B, parent(At)', B)
148148

149-
function mul_coefficients(Ac::Adjoint{T,<:QROperatorQ{QR,T}},B::AbstractVector{T};
150-
tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {QR,T}
149+
function mul_coefficients(Ac::Adjoint{T,<:QROperatorQ{<:Any,T}},B::AbstractVector{T};
150+
tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {T}
151151
mulpars(Ac,B,tolerance,maxlength)
152152
end
153-
function mul_coefficients!(Ac::Adjoint{T,<:QROperatorQ{QR,T}},B::AbstractVector{T};
154-
tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {QR,T}
153+
function mul_coefficients!(Ac::Adjoint{T,<:QROperatorQ{<:Any,T}},B::AbstractVector{T};
154+
tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {T}
155155
mulpars(Ac,B,tolerance,maxlength,Val(true))
156156
end
157157

158-
mul_coefficients(Ac::Adjoint{T,<:QROperatorQ{QR,T}},B::AbstractVector{V};opts...) where {QR,T,V} =
159-
mul_coefficients(Ac,AbstractVector{T}(B); opts...)
158+
mul_coefficients(Ac::Adjoint{<:Any,<:QROperatorQ}, B::AbstractVector; opts...) =
159+
mul_coefficients(Ac, AbstractVector{eltype(Ac)}(B); opts...)
160160

161161

162162
function *(Ac::Adjoint{<:Any,<:QROperatorQ}, b::AbstractVector; kwds...)
163163
A = parent(Ac)
164-
Fun(domainspace(A),mul_coefficients(Ac,coefficients(b,rangespace(A));kwds...))
164+
Fun(domainspace(A), mul_coefficients(Ac,coefficients(b,rangespace(A)); kwds...))
165165
end
166166

167167
function *(Ac::Adjoint{<:Any,<:QROperatorQ}, b; kwds...)
168168
A = parent(Ac)
169-
Fun(domainspace(A),mul_coefficients(Ac,coefficients(b,rangespace(A));kwds...))
169+
Fun(domainspace(A), mul_coefficients(Ac,coefficients(b,rangespace(A)); kwds...))
170170
end
171171

172172

@@ -200,26 +200,26 @@ end
200200
# QR
201201

202202
for TYP in (:Real,:Complex,:Number)
203-
@eval ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{T}; kwds...) where {CO,MT,T<:$TYP} =
204-
ldiv_coefficients(QR.R, mul_coefficients(QR.Q',b;kwds...))
205-
@eval ldiv_coefficients!(QR::QROperator{CO,MT,T},b::AbstractVector{T}; kwds...) where {CO,MT,T<:$TYP} =
206-
ldiv_coefficients!(QR.R, mul_coefficients!(QR.Q',b;kwds...))
203+
@eval ldiv_coefficients(QR::QROperator{<:Any,<:Any,T}, b::AbstractVector{T}; kwds...) where {T<:$TYP} =
204+
ldiv_coefficients(QR.R, mul_coefficients(QR.Q', b; kwds...))
205+
@eval ldiv_coefficients!(QR::QROperator{<:Any,<:Any,T}, b::AbstractVector{T}; kwds...) where {T<:$TYP} =
206+
ldiv_coefficients!(QR.R, mul_coefficients!(QR.Q', b; kwds...))
207207
end
208208

209209

210-
function ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{V};kwds...) where {CO,MT,T,V<:Number}
211-
TV = promote_type(T,V)
212-
ldiv_coefficients(convert(Operator{TV}, QR),convert(Vector{TV}, b);kwds...)
210+
function ldiv_coefficients(QR::QROperator, b::AbstractVector{<:Number}; kwds...)
211+
TV = promote_type(eltype(QR), eltype(b))
212+
ldiv_coefficients(convert(Operator{TV}, QR), convert(Vector{TV}, b); kwds...)
213213
end
214214

215-
function ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{V};kwds...) where {CO,MT,T<:Real,V<:Complex}
215+
function ldiv_coefficients(QR::QROperator{<:Any,<:Any,<:Real}, b::AbstractVector{<:Complex}; kwds...)
216216
a=ldiv_coefficients(QR,real(b);kwds...)
217217
b=im*ldiv_coefficients(QR,imag(b);kwds...)
218218
n=max(length(a),length(b))
219219
pad!(a,n)+pad!(b,n)
220220
end
221-
ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{V};kwds...) where {CO,MT,T<:Complex,V<:Real} =
222-
ldiv_coefficients(QR,Vector{T}(b);kwds...)
221+
ldiv_coefficients(QR::QROperator{<:Any,<:Any,<:Complex}, b::AbstractVector{<:Real}; kwds...) =
222+
ldiv_coefficients(QR, Vector{eltype(QR)}(b);kwds...)
223223

224224

225225
\(A::QROperator,b::Fun;kwds...) =

‎src/Space.jl‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -244,7 +244,7 @@ union_by_union_rule(a::AmbiguousSpace, b::Space) = b
244244
union_by_union_rule(a::Space, b::AmbiguousSpace) = a
245245

246246

247-
function union_by_union_rule(a::Space,b::Space)
247+
function union_by_union_rule(@nospecialize(a::Space), @nospecialize(b::Space))
248248
if spacescompatible(a,b)
249249
if isambiguous(domain(a))
250250
return b
@@ -259,20 +259,19 @@ function union_by_union_rule(a::Space,b::Space)
259259
union_rule(b,a)
260260
end
261261

262-
function union(a::Space, b::Space)
262+
function union(@nospecialize(a::Space), @nospecialize(b::Space))
263263
cr = union_by_union_rule(a,b)
264264
cr isa NoSpace || return cr
265265

266266
cspa=canonicalspace(a)
267267
cspb=canonicalspace(b)
268268
if cspa!=a || cspb!=b
269-
cr = union_by_union_rule(cspa,cspb)
269+
crc = union_by_union_rule(cspa,cspb)
270+
crc isa NoSpace || return crc
270271
end
271272
# TODO: Uncomment when Julia bug is fixed
272-
cr=maxspace(a,b) #Max space since we can convert both to it
273-
if !isa(cr,NoSpace)
274-
return cr
275-
end
273+
cr2=maxspace(a,b) #Max space since we can convert both to it
274+
cr2 isa NoSpace || return cr2
276275

277276
a ⊕ b
278277
end

‎src/testing.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ function testtransforms(S::Space;minpoints=1,invertibletransform=true)
3434
end
3535

3636
function testcalculus(S::Space;haslineintegral=true,hasintegral=true)
37-
for k=1:min(5,dimension(S))
37+
@testset for k=1:min(5,dimension(S))
3838
v = [zeros(k-1);1.0]
3939
f = Fun(S,v)
4040
@test abs(DefiniteIntegral()*f-sum(f)) < 100eps()
@@ -51,7 +51,7 @@ function testcalculus(S::Space;haslineintegral=true,hasintegral=true)
5151
end
5252

5353
function testmultiplication(spa,spb)
54-
for k=1:10
54+
@testset for k=1:10
5555
a = Fun(spa,[zeros(k-1);1.])
5656
M = Multiplication(a,spb)
5757
pts = ApproxFunBase.checkpoints(rangespace(M))

0 commit comments

Comments
 (0)