Supports approximating functions and solving differential equations on various higher dimensional domains such as disks and triangles

License: Other

Julia 100.00%

multivariateorthogonalpolynomials.jl's Introduction

MultivariateOrthogonalPolynomials.jl

This is an experimental package to add support for multivariate orthogonal polynomials on disks, spheres, triangles, and other simple geometries to ContinuumArrays.jl. At the moment it primarily supports triangles. For example, we can solve variable coefficient Helmholtz on the triangle with zero Dirichlet conditions as follows:

julia> using MultivariateOrthogonalPolynomials, StaticArrays, LinearAlgebra

julia> P = JacobiTriangle()
JacobiTriangle(0, 0, 0)

julia> x,y = first.(axes(P,1)), last.(axes(P,1));

julia> u = P * (P \ (exp.(x) .* cos.(y))) # Expand in Triangle OPs
JacobiTriangle(0, 0, 0) * [1.3365085377830084, 0.5687967596428205, -0.22812040274224554, 0.07733064070637755, 0.016169744493985644, -0.08714886622738759, 0.00338435674992512, 0.01220019521126353, -0.016867598915573725, 0.003930461395801074  …  ]

julia> u[SVector(0.1,0.2)] # Evaluate expansion
1.083141079608063

See the examples folder for more examples, including non-zero Dirichlet conditions, Neumann conditions, and piecing together multiple triangles. In particular, the examples from Olver, Townsend & Vasil 2019.

This code relies on Slevinsky's FastTransforms C library for calculating transforms between values and coefficients. At the moment the path to the compiled FastTransforms library is hard coded in c_transforms.jl.

References

S. Olver, A. Townsend & G.M. Vasil (2019), A sparse spectral method on triangles, arXiv:1902.04
S. Olver & Y. Xuan (2019), Orthogonal polynomials in and on a quadratic surface of revolution, arXiv:1906.12305
G.M. Vasil, K.J. Burns, D. Lecoanet, S. Olver, B.P. Brown & J.S. Oishi (2016), Tensor calculus in polar coordinates using Jacobi polynomials, J. Comp. Phys., 325: 53–73

multivariateorthogonalpolynomials.jl's People

Contributors

Stargazers

Watchers

Forkers

julienmln tsgut ikaroruan ioannispapapadopoulos danielvandh

multivariateorthogonalpolynomials.jl's Issues

Weak Laplacians on disks?

It would be nice to support weak laplacians. For example, it leads to a nice way of achieving Neumann condition via natural conditions: we would use Weighted(Zernike(1)) combined with the harmonic polynomials real(z^m) and imag(z^m) to pick up the rest of the polynomials. (see the new Neumann p-FEM test for an example in 1D) But it seems like there are a few options:

Support gradients of Zernike. But how do we incorporate the rotational invariance? And in a way so dot(∇, ∇*Weighted(Zernike(1))) still gives a diagonal operator?
Support tensor calculus a la Vasil et al. If I recall correctly this works on ∇_r and ∇_θ . This will successfully capture rotational invariance, but technically leaves the world of polynomials. That is, we can't directly view it as a vector orthogonal polynomial basis. Though perhaps (1) can be built from these.
Do weak Laplacian's directly without constructing gradients. In the case of m = 0 this seems straightforward via:

Enable CI

Should help keep the code alive. Will help identify breaking changes in dependencies

Better syntax for x,y = first.(xy),last.(xy)

Right now we can construct the basic monomials via e.g.

P = Zernike()
xy = axes(P,1)
x,y = first.(xy),last.(xy)

which already isn't very nice but even worse if we do 3D do we really want to write:

P = ZernikeBall()
xyz = axes(P,1)
x,y,z = first.(xyz),getindex.(xyz,2),last.(xyz)

The catch: what's a good name???? Best I could come up with is

P = Zernike()
x,y,z = coordinateaxes(P,1)

Lowering bug?

@dlfivefifty The lowering Zernike(0,1) \ Weighted(Zernike(0,1)) is a factor of 2 off. This also seems to be confirmed in the tests..?

MultivariateOrthogonalPolynomials.jl/test/test_disk.jl

Line 303 in cb0fb24

 @test 2w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike(1)[xy,Block.(1:7)]'*L[Block.(1:7),Block.(1:5)] 

Is there a reason for this?

Jacobi matrices for RectPolynomial

I just added these in commit b2f95fe, but they cannot be instantiated like T² \ (x .* T²), like those on the Dunkl-Xu disk or the JacobiTriangle for reasons I don't understand. The code is so simple that the error must be somewhere else. I suspect it's because the Jacobi matrices for the basis are lazy tridiagonals instead of banded-block-banded matrices for the nontrivial Koornwinder classes, which means the error is not really in this code.

julia> using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials

julia> J = JacobiTriangle(0, 0, 0)
JacobiTriangle(0, 0, 0)

julia> jacobimatrix(Val{1}(), J); # works

julia> x,y = first.(axes(J,1)), last.(axes(J,1))
(first.(Inclusion(UnitSimplex(Val(2)))), last.(Inclusion(UnitSimplex(Val(2)))))

julia> jacobimatrix(Val{1}(), J)^C

julia> J \ (x .* J); # also works

julia> T² = RectPolynomial(ChebyshevT(), ChebyshevT())
ChebyshevT() ⊗ ChebyshevT()

julia> jacobimatrix(Val{1}(), T²); # works. This is what I just added

julia> x,y = first.(axes(T²,1)), last.(axes(T²,1))
(first.(Inclusion((-1.0 .. 1.0 (Chebyshev)) × (-1.0 .. 1.0 (Chebyshev)))), last.(Inclusion((-1.0 .. 1.0 (Chebyshev)) × (-1.0 .. 1.0 (Chebyshev)))))

julia> T² \ (x .* T²) # errs
ERROR: DimensionMismatch: Second axis of A, 1:1:ℵ₀, and first axis of B, Inclusion((-1.0 .. 1.0 (Chebyshev)) × (-1.0 .. 1.0 (Chebyshev))) must match
Stacktrace:
 [1] _check_mul_axes(A::RectPolynomial{Float64, Tuple{…}}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:95
 [2] check_mul_axes(::RectPolynomial{Float64, Tuple{…}}, ::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:111
 [3] instantiate
   @ ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:123 [inlined]
 [4] materialize
   @ ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:127 [inlined]
 [5] mul
   @ ~/.julia/packages/ArrayLayouts/i3v9G/src/mul.jl:128 [inlined]
 [6] *(A::RectPolynomial{Float64, Tuple{…}}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/ZCn0F/src/matmul.jl:23
 [7] broadcasted(::QuasiArrays.LazyQuasiArrayStyle{…}, ::typeof(*), x::QuasiArrays.BroadcastQuasiVector{…}, P::RectPolynomial{…})
   @ HarmonicOrthogonalPolynomials ~/.julia/packages/HarmonicOrthogonalPolynomials/ViYwn/src/multivariateops.jl:45
 [8] broadcasted(::typeof(*), ::QuasiArrays.BroadcastQuasiVector{…}, ::RectPolynomial{…})
   @ Base.Broadcast ./broadcast.jl:1347
 [9] top-level scope
   @ REPL[119]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia>

Ridge polynomials (on the disk)?

These look like P_{n,k}(x,y) = U_n(x * cos(k*pi/(n+1)) + y * sin(k*pi/(n+1))), orthogonal on the unweighted unit disk.

Sparse differentiation is trivial (using ultraspherical generalizations), but the Jacobi operators are not obvious.

How to instantiate more complicated continuous operators?

Such a cool repository has no open issues?

Here's a MWE of what I envision working:

julia> using MultivariateOrthogonalPolynomials

julia> a, b, c = 0, 0, 0
(0, 0, 0)

julia> P = JacobiTriangle(a, b, c)
JacobiTriangle(0, 0, 0)

julia> Q = JacobiTriangle(a+1, b+1, c+1)
JacobiTriangle(1, 1, 1)

julia> ∂x = PartialDerivative{1}(P)
PartialDerivative{1, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}}

julia> ∂y = PartialDerivative{2}(P)
PartialDerivative{2, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}}

julia> x, y = first.(axes(P, 1)), last.(axes(P, 1))
(first.(Inclusion(DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}())), last.(Inclusion(DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}())))

julia> x*∂y - y*∂x
ERROR: DimensionMismatch("Second axis of A, Base.OneTo(1), and first axis of B, Inclusion(JacobiTriangle(0, 0, 0)) must match")
Stacktrace:
 [1] _check_mul_axes(A::QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}}, B::PartialDerivative{2, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:87
 [2] check_mul_axes(::QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}}, ::PartialDerivative{2, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:89
 [3] instantiate
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:101 [inlined]
 [4] materialize(M::ArrayLayouts.Mul{LazyArrays.BroadcastLayout{typeof(first)}, QuasiArrays.QuasiLazyLayout, QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}}, PartialDerivative{2, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:105
 [5] mul
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:106 [inlined]
 [6] *(A::QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}}, B::PartialDerivative{2, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/bxqu9/src/matmul.jl:23
 [7] top-level scope
   @ REPL[48]:1

Dots don't seem to work but the error message of an ambiguity warning appears a bit more hopeful:

julia> x.*∂y - y.*∂x
ERROR: MethodError: *(::Infinities.InfiniteCardinal{1}, ::Infinities.InfiniteCardinal{0}) is ambiguous. Candidates:
  *(a::Integer, b::Infinities.InfiniteCardinal) in Infinities at /Users/mikael/.julia/packages/Infinities/lXuMq/src/cardinality.jl:134
  *(a::Number, b::Infinities.InfiniteCardinal) in Infinities at /Users/mikael/.julia/packages/Infinities/lXuMq/src/cardinality.jl:139
  *(a::Infinities.InfiniteCardinal, b::Integer) in Infinities at /Users/mikael/.julia/packages/Infinities/lXuMq/src/cardinality.jl:141
  *(a::Integer, b::Integer) in Base at int.jl:919
  *(a::Infinities.InfiniteCardinal, b::Number) in Infinities at /Users/mikael/.julia/packages/Infinities/lXuMq/src/cardinality.jl:142
Possible fix, define
  *(::Infinities.InfiniteCardinal, ::Infinities.InfiniteCardinal)
Stacktrace:
  [1] prod(x::Tuple{Infinities.InfiniteCardinal{1}, Infinities.InfiniteCardinal{0}})
    @ Base ./tuple.jl:480
  [2] length
    @ ~/.julia/packages/QuasiArrays/bxqu9/src/abstractquasiarray.jl:153 [inlined]
  [3] cardinality(d::JacobiTriangle{Float64, Int64})
    @ QuasiArrays ~/.julia/packages/QuasiArrays/bxqu9/src/QuasiArrays.jl:75
  [4] length(S::QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}})
    @ QuasiArrays ~/.julia/packages/QuasiArrays/bxqu9/src/indices.jl:182
  [5] _bcsm(a::QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}, b::QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}})
    @ Base.Broadcast ./broadcast.jl:503
  [6] _bcs1(a::QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}, b::QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}})
    @ Base.Broadcast ./broadcast.jl:501
  [7] _bcs(shape::Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}, newshape::Tuple{QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}})
    @ Base.Broadcast ./broadcast.jl:495
  [8] broadcast_shape(::Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}, ::Tuple{QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}})
    @ Base.Broadcast ./broadcast.jl:489
  [9] combine_axes
    @ ./broadcast.jl:484 [inlined]
 [10] instantiate
    @ ./broadcast.jl:266 [inlined]
 [11] materialize(bc::Base.Broadcast.Broadcasted{QuasiArrays.LazyQuasiArrayStyle{2}, Nothing, typeof(*), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArrays.SVector{2, Float64}, DomainSets.FixedUnitSimplex{StaticArrays.SVector{2, Float64}, :closed}}}}, PartialDerivative{2, Float64, QuasiArrays.Inclusion{Float64, JacobiTriangle{Float64, Int64}}}}})
    @ Base.Broadcast ./broadcast.jl:883
 [12] top-level scope
    @ REPL[49]:1

Fine, I go along and materialize my four operators independently: partial derivatives and multiplication

julia> Dx = Q \ (∂x * P)
(ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) * (ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyBandedMatrices.BlockBroadcastArray{Float64, 2, typeof(hcat), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) with indices 1:1:ℵ₀×1:1:ℵ₀:
  ⋅   │  3.0  1.0  │  -1.33333  0.666667
 ─────┼────────────┼────────────────────  …  
  ⋅   │   ⋅    ⋅   │   3.33333  0.833333     
  ⋅   │   ⋅    ⋅   │    ⋅       2.5     
 ─────┼────────────┼────────────────────     
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅           
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅           
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅      
 ─────┼────────────┼────────────────────  …  
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅           
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅           
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅           
  ⋅   │   ⋅    ⋅   │    ⋅        ⋅      
 ─────┼────────────┼────────────────────     
 ⋮                                  ⋱  

julia> Dy = Q \ (∂y * P)
(ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, LazyBandedMatrices.BlockVcat{Float64, 2, Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#7#9", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64, Int64}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#8#10", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) * (ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Float64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) with indices 1:1:ℵ₀×1:1:ℵ₀:
  ⋅   │   ⋅   2.0  │   ⋅   1.33333   ⋅   │   ⋅    ⋅  
 ─────┼────────────┼─────────────────────┼───────────  …  
  ⋅   │   ⋅    ⋅   │   ⋅   1.66667   ⋅   │   ⋅   1.25     
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅       3.0  │   ⋅    ⋅  
 ─────┼────────────┼─────────────────────┼───────────     
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅   1.5      
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅    ⋅       
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅    ⋅  
 ─────┼────────────┼─────────────────────┼───────────  …  
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅    ⋅       
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅    ⋅       
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅    ⋅       
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅        ⋅   │   ⋅    ⋅  
 ─────┼────────────┼─────────────────────┼───────────     
 ⋮                            ⋮               ⋱  

julia> Mx = Q \ (x .* Q)
(ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) * (ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, LazyBandedMatrices.BlockVcat{Float64, 2, Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#7#9", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64, Int64}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#8#10", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) with indices 1:1:ℵ₀×1:1:ℵ₀:
 0.333333  │  0.190476   ⋅     │   ⋅      
 ──────────┼───────────────────┼──────────  …  
 0.166667  │  0.416667   ⋅     │  0.208333     
  ⋅        │   ⋅        0.25   │   ⋅      
 ──────────┼───────────────────┼──────────     
  ⋅        │  0.214286   ⋅     │  0.45         
  ⋅        │   ⋅        0.125  │   ⋅           
  ⋅        │   ⋅         ⋅     │   ⋅      
 ──────────┼───────────────────┼──────────  …  
  ⋅        │   ⋅         ⋅     │  0.233333     
  ⋅        │   ⋅         ⋅     │   ⋅           
  ⋅        │   ⋅         ⋅     │   ⋅           
  ⋅        │   ⋅         ⋅     │   ⋅      
 ──────────┼───────────────────┼──────────     
 ⋮                                    ⋱  

julia> My = Q \ (y .* Q)
(ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) * (ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}) with indices 1:1:ℵ₀×1:1:ℵ₀:
  0.333333   │  -0.0952381   0.0952381   │
 ────────────┼───────────────────────────┼  …  
 -0.0833333  │   0.291667   -0.0416667   │     
  0.25       │  -0.125       0.375       │
 ────────────┼───────────────────────────┼     
   ⋅         │  -0.107143    0.00714286  │     
   ⋅         │   0.1875     -0.0625      │     
   ⋅         │    ⋅          0.266667    │
 ────────────┼───────────────────────────┼  …  
   ⋅         │    ⋅           ⋅          │     
   ⋅         │    ⋅           ⋅          │     
   ⋅         │    ⋅           ⋅          │     
   ⋅         │    ⋅           ⋅          │
 ────────────┼───────────────────────────┼     
  ⋮                                   ⋱

But Julia hangs on taking the difference below:

julia> Mx*Dy - My*Dx
^CERROR: InterruptException:
Stacktrace:
  [1] getindex
    @ ~/.julia/packages/LazyArrays/CpDkh/src/lazyapplying.jl:231 [inlined]
  [2] _getindex
    @ ./abstractarray.jl:1214 [inlined]
  [3] getindex
    @ ./abstractarray.jl:1170 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:614 [inlined]
  [5] _getindex
    @ ./broadcast.jl:644 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:620 [inlined]
  [7] getindex
    @ ./broadcast.jl:575 [inlined]
  [8] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] copyto!
    @ ./broadcast.jl:983 [inlined]
 [11] copyto!
    @ ./broadcast.jl:936 [inlined]
 [12] copy
    @ ./broadcast.jl:908 [inlined]
 [13] materialize
    @ ./broadcast.jl:883 [inlined]
 [14] broadcast_preserving_zero_d
    @ ./broadcast.jl:872 [inlined]
 [15] -(A::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}, BlockBandedMatrices.BandedBlockBandedMatrix{Float64, LazyBandedMatrices.BlockVcat{Float64, 2, Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#7#9", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64, Int64}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#8#10", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}, BlockBandedMatrices.BandedBlockBandedMatrix{Float64, LazyBandedMatrices.BlockVcat{Float64, 2, Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#7#9", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64, Int64}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, MultivariateOrthogonalPolynomials.var"#8#10", Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}, BlockBandedMatrices.BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Float64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}, B::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}, BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}, BlockBandedMatrices.BandedBlockBandedMatrix{Float64, BlockArrays.PseudoBlockMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}}}, Tuple{BlockArrays.BlockedUnitRange{StepRange{Int64, Int64}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}, BlockBandedMatrices.BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyBandedMatrices.BlockBroadcastArray{Float64, 2, typeof(hcat), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}})
    @ Base ./arraymath.jl:39
 [16] top-level scope
    @ REPL[54]:1

Of course, the fallback of indexing the four operators to create arrays works:

julia> Mx[1:10,1:10]*Dy[1:10,1:10] - My[1:10,1:10]*Dx[1:10,1:10]
10×10 Matrix{Float64}:
 0.0  -1.0    0.333333   0.761905  …   0.369048  -0.0357143  -0.107143
 0.0   0.25   0.416667  -1.08333       0.479167  -0.3125      0.0625
 0.0  -0.75  -0.25       0.75          0.5625     0.4375     -0.6875
 0.0   0.0    0.0        0.357143      0.966429   0.233571   -0.0192857
 0.0   0.0    0.0       -0.625        -0.91875    0.76875     0.20625
 0.0   0.0    0.0        0.0       …   0.48      -1.28       -0.32
 0.0   0.0    0.0        0.0           0.4025    -0.0525      0.0025
 0.0   0.0    0.0        0.0           0.116667   0.616667   -0.0261905
 0.0   0.0    0.0        0.0          -0.56      -0.14        0.54
 0.0   0.0    0.0        0.0           0.0       -0.75       -0.535714

But now I'm confused how I've thrown away all structure and ended up with a matrix (rather than a BandedBlockBandedMatrix).

Fix plotting

Plotting on the disk using Plots.jl is basically broken:

I suspect the best thing to do is forget about Plots.jl and focus on Makie.jl which I believe has advanced a lot. The catch is we probably don't want to make Makie.jl a dependency and there isn't a lightweight recipes package like in Plots.jl.

We also need to get it to use inverse transforms to compute the values.

Laplacian vs. Laplace?

This is a minor issue, but most other operators are nouns: Derivative, Integral, etc... This convention would get more strange with union-ed names: LaplacianBeltramian?

Dependency with HOPs reversed?

I'm struggling to understand the dependency order between HarmonicOrthogonalPolynomials and MultivariateOrthogonalPolynomials. To me, a harmonic OP is an OP before it's harmonic, and since harmonicity is trivial for univariate polynomials that means multivariate orthogonality among polynomials is more fundamental an idea.

Maybe this is just a temporary defect in the way the repositories evolved?

JacobiTriangle() expansions failing in strange way

There is some strange behavior when expanding functions on triangles, apparently when it comes to globally defined variables? Here is a minimal example:

julia> using MultivariateOrthogonalPolynomials

julia> T = JacobiTriangle()
JacobiTriangle(0, 0, 0)

julia> xy = axes(T,1)
Inclusion(UnitSimplex(Val(2)))

julia> x, y = first.(xy), last.(xy)
(first.(Inclusion(UnitSimplex(Val(2)))), last.(Inclusion(UnitSimplex(Val(2)))))

julia> α = 2.
2.0

julia> K(x,y) = α*x*y
K (generic function with 1 method)

julia> K(α,x,y) = α*x*y
K (generic function with 2 methods)

julia> T \ (@. α*x*y)  # this works
ℵ₀-blocked ℵ₀-element LazyBandedMatrices.DiagTrav{Float64, 2, LazyArrays.ApplyArray{Float64, 2, typeof(Base.setindex), Tuple{FillArrays.Zeros{Float64, 2, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Matrix{Float64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}:
  0.1666666666666667    
 ───────────────────────
  0.06666666666666667   
  0.19999999999999998   
 ───────────────────────
 -0.09999999999999999   
  0.19999999999999993   
 -6.699811370772778e-18 
 ───────────────────────
  0.0                   
 -6.148984834775253e-18 
 -1.6452652161292338e-18
  ⋮

julia> T \ (@. K(α,x,y))  # this works
ℵ₀-blocked ℵ₀-element LazyBandedMatrices.DiagTrav{Float64, 2, LazyArrays.ApplyArray{Float64, 2, typeof(Base.setindex), Tuple{FillArrays.Zeros{Float64, 2, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Matrix{Float64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}}:
  0.1666666666666667    
 ───────────────────────
  0.06666666666666667   
  0.19999999999999998   
 ───────────────────────
 -0.09999999999999999   
  0.19999999999999993   
 -6.699811370772778e-18 
 ───────────────────────
  0.0                   
 -6.148984834775253e-18 
 -1.6452652161292338e-18
  ⋮

julia> T \ (@. K(x,y))  # this fails
ERROR: MethodError: no method matching *(::FastTransforms.FTPlan{Float64, 2, FastTransforms.TRIANGLEANALYSIS}, ::Matrix{Any})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:591
  *(::Union{SparseArrays.AbstractSparseMatrixCSC{TA, Ti}, SubArray{TA, 2, <:SparseArrays.AbstractSparseMatrixCSC{TA, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, I}} where I<:AbstractUnitRange} where Ti, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}) where TA at ~/Documents/Backends/julia-1.8.5/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:52
  *(::LinearAlgebra.Hermitian{<:Any, <:ArrayLayouts.LayoutMatrix}, ::AbstractMatrix) at ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:263
  ...
Stacktrace:
  [1] *(T::MultivariateOrthogonalPolynomials.TriPlan{Float64}, F::Matrix{Any})
    @ MultivariateOrthogonalPolynomials ~/.julia/packages/MultivariateOrthogonalPolynomials/tNdGo/src/triangle.jl:661
  [2] \(a::ContinuumArrays.TransformFactorization{Any, Matrix{StaticArraysCore.SVector{2, Float64}}, MultivariateOrthogonalPolynomials.TriPlan{Float64}}, b::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/plans.jl:26
  [3] transform_ldiv(V::QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64})
    @ HarmonicOrthogonalPolynomials ~/.julia/packages/HarmonicOrthogonalPolynomials/51vXS/src/multivariateops.jl:91
  [4] transform_ldiv(A::QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:258
  [5] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}}, #unused#::Base.OneTo{Int64})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:289
  [6] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:290
  [7] copy(L::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:291
  [8] materialize(M::ArrayLayouts.Ldiv{ContinuumArrays.SubBasisLayout, LazyArrays.BroadcastLayout{typeof(K)}, QuasiArrays.SubQuasiArray{Float64, 2, JacobiTriangle{Float64, Int64}, Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, BlockArrays.BlockSlice{BlockArrays.BlockRange{1, Tuple{Base.OneTo{Int64}}}, BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, Base.OneTo{Int64}}}}}, false}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22
  [9] materialize
    @ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22 [inlined]
 [10] #ldiv#18
    @ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86 [inlined]
 [11] ldiv
    @ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86 [inlined]
 [12] \
    @ ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:34 [inlined]
 [13] adaptivetransform_ldiv(A::JacobiTriangle{Float64, Int64}, f::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
    @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/ix48P/src/adaptivetransform.jl:35
 [14] transform_ldiv(A::JacobiTriangle{Float64, Int64}, f::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}, #unused#::Tuple{Infinities.InfiniteCardinal{1}, Infinities.InfiniteCardinal{0}})
    @ ClassicalOrthogonalPolynomials ~/.julia/packages/ClassicalOrthogonalPolynomials/ix48P/src/adaptivetransform.jl:2
 [15] transform_ldiv(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:258
 [16] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}}, #unused#::Base.OneTo{Int64})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:289
 [17] transform_ldiv_if_columns(L::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:290
 [18] copy(L::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
    @ ContinuumArrays ~/.julia/packages/ContinuumArrays/pWRm2/src/bases/bases.jl:291
 [19] #materialize#11
    @ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22 [inlined]
 [20] materialize(M::ArrayLayouts.Ldiv{HarmonicOrthogonalPolynomials.MultivariateOPLayout{2}, LazyArrays.BroadcastLayout{typeof(K)}, JacobiTriangle{Float64, Int64}, QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:22
 [21] ldiv(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}}; kwds::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86
 [22] ldiv
    @ ~/.julia/packages/ArrayLayouts/4IG3b/src/ldiv.jl:86 [inlined]
 [23] \(A::JacobiTriangle{Float64, Int64}, B::QuasiArrays.BroadcastQuasiVector{Any, typeof(K), Tuple{QuasiArrays.BroadcastQuasiVector{Float64, typeof(first), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}, QuasiArrays.BroadcastQuasiVector{Float64, typeof(last), Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}}}})
    @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/matmul.jl:34
 [24] top-level scope
    @ ~/Documents/Projects/SparseVolterraExamples.jl/examples/Example 0 - Some linear Volterra integral equation examples.jl:82

Any idea what is happening here?

Add ZeroTo (somewhere)?

One of the greatest things is the notation to fill in multivariate operators: given the mortared types of n and k, broadcasted operations Just Work (ok, probably that took a lot of effort).

Wouldn't it be nice to index from zero to infinity to match the math formulas given that polynomial degree starts at 0?

New Operator: AngularMomentum

This one is just about as fundamental as the Laplacian. For complex spherical harmonics and Zernike polynomials, it's diagonal. Other multivariate OPs may have nontrivial representations. In polar coordinates, it's ∂/∂θ; in spherical coordinates, it's ∂/∂φ. Or in both systems, it's x ∂/∂y - y ∂/∂x. Since it's spectrum "separates orders," it is complementary to the Laplacian which "separates degrees."

Missing triangle functionality (Laplacian, restriction)

Am I correct in believing that the Laplacian and edge restriction operators on the triangle have not been implemented in the ClassicalOPs version? I distinctly remember that functionality working in ApproxFun.

Are these actually the adjoint Jacobi matrices? I.e. multiplication operators?

MultivariateOrthogonalPolynomials.jl/src/triangle.jl

Lines 189 to 192 in 04704a8

 ## Jacobi matrix 

 jacobimatrix(::Val{1}, P::JacobiTriangle) = Lx(P.a+1, P.b, P.c) * Rx(P.a, P.b, P.c) 

 jacobimatrix(::Val{2}, P::JacobiTriangle) = Ly(P.a, P.b+1, P.c) * Ry(P.a, P.b, P.c)

Just wondering. (I don't really care, do you?)

Make evaluation of triangle functions type stable

Examples using old ApproxFun package structure

In some places the package still attempts to import functions which have been moved from ApproxFun to ApproxFunBase and so on. A few examples I quickly found are in

https://github.com/JuliaApproximation/MultivariateOrthogonalPolynomials.jl/blob/038850f57d2fd1aa9e50f460e2b391fcb9b9b652/examples/L-shaped%20Laplace.jl

and

https://github.com/JuliaApproximation/MultivariateOrthogonalPolynomials.jl/blob/ea84cf3f15a18f2eccbb510938bc1ff45beee423/examples/triangleexamples.jl

That means even a consistent fresh install won't be able to run them as it is and it will also throw warnings upon compiling the package. We should probably go through and update them to the new ApproxFun structure.

How to make rotationally invariant solves fast?

Rotationally invariant operators are typically represented by ModalInterlace which combines matrices acting only on Fourier modes into a diagonal-block-banded matrix (that is, block banded matrix with diagonal blocks, i.e. subblockbandwidths are (0,0)).

At the moment we aren't taking advantage of this structure. There are two approaches:

Do it at the level of factorizing a BandedBlockBandedMatrix by checking whether the subblockbandwidths are (0,0). The easiest way would be to copy to a banded matrix, do a QR, and then copy the data back.
Introduce a DiagonalBlockBandedMatrix though this doesn't seem to have any obvious benefits over (1).
Do it directly on a ModalInterlace. This has the benefit of being able to call qr on each of the banded operators separately (and potentially taking advantage of @threads). Though probably (1) could also be parallelised.

@TSGut @ioannisPApapadopoulos any thoughts? Do you think this is important right now? It won't help with the variable coefficients problem but could be important for fractional DEs and time-stepping.

Triangle examples in examples/triangleexamples.jl fail

The examples here fail (below is Lines 1-12 from that file)

julia> using ApproxFun, MultivariateOrthogonalPolynomials, BlockArrays, SpecialFunctions, FillArrays, Plots
       import ApproxFun: blockbandwidths, Vec, PiecewiseSegment
       import MultivariateOrthogonalPolynomials: DirichletTriangle
       S = TriangleWeight(1,1,1,JacobiTriangle(1,1,1))
       Δ = Laplacian(S)
WARNING: could not import ApproxFun.blockbandwidths into Main
WARNING: could not import ApproxFun.Vec into Main
WARNING: could not import MultivariateOrthogonalPolynomials.DirichletTriangle into Main
ERROR: MethodError: no method matching TriangleWeight(::Int64, ::Int64, ::Int64, ::JacobiTriangle{Float64, Int64})

Closest candidates are:
  TriangleWeight(::T, ::T, ::T) where T
   @ MultivariateOrthogonalPolynomials C:\Users\User\.julia\packages\MultivariateOrthogonalPolynomials\dZw1Q\src\triangle.jl:111

Stacktrace:
 [1] top-level scope
   @ Untitled-1:4

I believe in this case the fix is to do S = WeightedTriangle(1, 1, 1)? This error appears later in the script as well (might be some others too - Laplacian(S) also errors -, but I haven't run it that far). I could try and fix it up later once I've looked into it all a bit more.

Use fast triangle transform

@MikaelSlevinsky If you have example code that I could copy, that would be great

Either: PartialDerivative{k}(<:Basis) or Inclusion(::Basis) should err (or both)

julia> P = JacobiTriangle()
JacobiTriangle(0, 0, 0)

julia> Q = JacobiTriangle(2, 2, 2)
JacobiTriangle(2, 2, 2)

julia> ∂x = PartialDerivative{1}(P)
PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}

julia> Q \ (∂x^2 * P)
ERROR: DimensionMismatch("Second axis of A, Inclusion(JacobiTriangle(0, 0, 0)), and first axis of B, Inclusion(DomainSets.FixedUnitSimplex{SVector{2, Float64}, :closed}()) must match")
Stacktrace:
 [1] _check_mul_axes(A::PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, B::ApplyQuasiMatrix{Float64, typeof(*), Tuple{JacobiTriangle{Float64, Int64}, BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyBandedMatrices.BlockBroadcastArray{Float64, 2, typeof(hcat), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:87
 [2] check_mul_axes(::PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, ::ApplyQuasiMatrix{Float64, typeof(*), Tuple{JacobiTriangle{Float64, Int64}, BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyBandedMatrices.BlockBroadcastArray{Float64, 2, typeof(hcat), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:89
 [3] instantiate
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:101 [inlined]
 [4] materialize
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:105 [inlined]
 [5] mul
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:106 [inlined]
 [6] *(A::PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, B::ApplyQuasiMatrix{Float64, typeof(*), Tuple{JacobiTriangle{Float64, Int64}, BandedBlockBandedMatrix{Float64, Adjoint{Float64, LazyBandedMatrices.BlockBroadcastArray{Float64, 2, typeof(hcat), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockVector{Int64, LazyArrays.BroadcastVector{StepRange{Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, BlockedUnitRange{RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/bxqu9/src/matmul.jl:23
 [7] *(App::ApplyQuasiMatrix{Float64, typeof(^), Tuple{PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, Int64}}, b::JacobiTriangle{Float64, Int64})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/bxqu9/src/lazyquasiarrays.jl:231
 [8] top-level scope
   @ REPL[72]:1

julia> Q \ (∂x.^2 * P)
ERROR: DimensionMismatch("Second axis of A, Inclusion(JacobiTriangle(0, 0, 0)), and first axis of B, Inclusion(DomainSets.FixedUnitSimplex{SVector{2, Float64}, :closed}()) must match")
Stacktrace:
 [1] _check_mul_axes(A::BroadcastQuasiMatrix{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, Base.RefValue{Val{2}}}}, B::JacobiTriangle{Float64, Int64})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:87
 [2] check_mul_axes(::BroadcastQuasiMatrix{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, Base.RefValue{Val{2}}}}, ::JacobiTriangle{Float64, Int64})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:89
 [3] instantiate
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:101 [inlined]
 [4] materialize(M::ArrayLayouts.Mul{LazyArrays.BroadcastLayout{typeof(Base.literal_pow)}, ContinuumArrays.BasisLayout, BroadcastQuasiMatrix{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, Base.RefValue{Val{2}}}}, JacobiTriangle{Float64, Int64}})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:105
 [5] mul
   @ ~/.julia/packages/ArrayLayouts/mTQvg/src/mul.jl:106 [inlined]
 [6] *(A::BroadcastQuasiMatrix{Float64, typeof(Base.literal_pow), Tuple{Base.RefValue{typeof(^)}, PartialDerivative{1, Float64, Inclusion{Float64, JacobiTriangle{Float64, Int64}}}, Base.RefValue{Val{2}}}}, B::JacobiTriangle{Float64, Int64})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/bxqu9/src/matmul.jl:23
 [7] top-level scope
   @ REPL[73]:1

julia>

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

Weighted ∂ₓ on triangles errs

julia> using MultivariateOrthogonalPolynomials
[ Info: Precompiling MultivariateOrthogonalPolynomials [4f6956fd-4f93-5457-9149-7bfc4b2ce06d]

julia> MultivariateOrthogonalPolynomials.Wx(1,1,1)
ERROR: ArgumentError: Data matrix must have number of row blocks equal to number of block bands
Stacktrace:
 [1] check_data_sizes(data::LazyBandedMatrices.BlockVcat{Float64, 2, Tuple{LinearAlgebra.Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRangeLen{Int64, Int64, Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}, LinearAlgebra.Adjoint{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{LazyArrays.BroadcastVector{Int64, typeof(*), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{LazyArrays.BroadcastVector{Int64, typeof(-), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}, Type{FillArrays.Fill}, Tuple{InfiniteArrays.OneToInf{Int64}, InfiniteArrays.OneToInf{Int64}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}}}, Int64}}}}, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{BlockArrays.BlockVector{Int64, LazyArrays.BroadcastVector{StepRangeLen{Int64, Int64, Int64, Int64}, typeof(*), Tuple{Int64, LazyArrays.BroadcastVector{Base.OneTo{Int64}, Type{Base.OneTo}, Tuple{InfiniteArrays.OneToInf{Int64}}}}}, Tuple{BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Int64}}}}}}}, raxis::BlockArrays.BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}, ::Tuple{Int64, Int64}, ::Tuple{Int64, Int64})
   @ BlockBandedMatrices ~/.julia/packages/BlockBandedMatrices/sYxap/src/BandedBlockBandedMatrix.jl:3
 [2] _BandedBlockBandedMatrix
   @ ~/.julia/packages/BlockBandedMatrices/sYxap/src/BandedBlockBandedMatrix.jl:37 [inlined]
 [3] Wx(a::Int64, b::Int64, c::Int64)
   @ MultivariateOrthogonalPolynomials ~/.julia/dev/MultivariateOrthogonalPolynomials/src/triangle.jl:120
 [4] top-level scope
   @ REPL[2]:1

Which notation is best?

There are multiple notations used to create multivariate operators:

MultivariateOrthogonalPolynomials.jl/src/triangle.jl

Lines 89 to 128 in 04704a8

 function Wx(a,b,c) 

 n = mortar(Fill.(oneto(∞),oneto(∞))) 

 k = mortar(Base.OneTo.(oneto(∞))) 

 dat = BlockVcat( 

 ((k .+ (c-1)) .* ( k .- n .- 1 ) ./ (2k .+ (b+c-1)))', 

 (k .* (k .- n .- a) ./ (2k .+ (b+c-1)))' 

 ) 

 _BandedBlockBandedMatrix(dat, axes(k,1), (1,-1), (1,0)) 

 end 

 function Rx(a,b,c) 

 n = mortar(Fill.(oneto(∞),oneto(∞))) 

 k = mortar(Base.OneTo.(oneto(∞))) 

 dat = BlockHcat( 

 BroadcastVector((n,k,bc1,abc) -> (n + k + bc1) / (2n + abc), n, k, b+c-1, a+b+c), 

 BroadcastVector((n,k,abc) -> (n + k + abc) / (2n + abc), n, k, a+b+c) 

 ) 

 _BandedBlockBandedMatrix(dat', axes(k,1), (0,1), (0,0)) 

 end 

 function Ry(a,b,c) 

 n = mortar(Fill.(oneto(∞),oneto(∞))) 

 k = mortar(Base.OneTo.(oneto(∞))) 

 dat = PseudoBlockArray(Vcat( 

 ((k .+ (c-1) ) .* (n .+ k .+ (b+c-1)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))', 

 ((k .- n .- a ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))', 

 ((k .+ (c-1) ) .* (k .- n .- 1) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))', 

 ((n .+ k .+ (a+b+c) ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))' 

 ), (blockedrange(Fill(2,2)), axes(n,1))) 

 _BandedBlockBandedMatrix(dat, axes(k,1), (0,1), (0,1)) 

 end 

 function Lx(a,b,c) 

 n = mortar(Fill.(oneto(∞),oneto(∞))) 

 k = mortar(Base.OneTo.(oneto(∞))) 

 dat = PseudoBlockArray(Vcat( 

 ((n .- k .+ a) ./ (2n .+ (a+b+c)))', 

 ((n .- k .+ 1) ./ (2n .+ (a+b+c)))' 

 ), (blockedrange(Fill(1,2)), axes(n,1))) 

 _BandedBlockBandedMatrix(dat, axes(k,1), (1,0), (0,0)) 

 end

When to use the dot notation and when to use BroadcastVector? Is it better to use @. or is that wasteful because it creates some unnecessary broadcasting?

(My impression is that these are due to successive invasions of new notation and not all operators are fully updated, but I could be way off.)

Generalized eigenvalue problems on a disc

Timon @TSGut has done the heavy lifting and has implemented Jacobi matrices for Zernike polynomials on a disc.

Alberto @APaganini could you remind us what exactly you required? If I remember correctly, as a first step you wanted to solve

-Laplace(u)(x,y) = q(x,y) u(x,y)

on a disc with zero Dirichlet bcs. What should we pick for q(x,y)? (@dlfivefifty)

Encountering BoundsErrors with ClenshawMultiplication operators

Running simple things like this in fresh install:

T = JacobiTriangle()
(x,y) = Fun(identity, T)
Multiplication(y, T)

returns

BoundsError: attempt to access 1×2 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}} at index [2, 1]

The error seems to originate somewhere in the Clenshaw algorithm step.

Would like to know if this reproduces for other people as well. Was this bug introduced during the LazyArrays version support change? I'm running all packages at the highest consistent version, i.e. LazyArrays v0.10.0 etc..

My older working instances run the above code without problems.

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.

	## Jacobi matrix

	jacobimatrix(::Val{1}, P::JacobiTriangle) = Lx(P.a+1, P.b, P.c) * Rx(P.a, P.b, P.c)
	jacobimatrix(::Val{2}, P::JacobiTriangle) = Ly(P.a, P.b+1, P.c) * Ry(P.a, P.b, P.c)

	function Wx(a,b,c)
	n = mortar(Fill.(oneto(∞),oneto(∞)))
	k = mortar(Base.OneTo.(oneto(∞)))
	dat = BlockVcat(
	((k .+ (c-1)) .* ( k .- n .- 1 ) ./ (2k .+ (b+c-1)))',
	(k .* (k .- n .- a) ./ (2k .+ (b+c-1)))'
	)
	_BandedBlockBandedMatrix(dat, axes(k,1), (1,-1), (1,0))
	end

	function Rx(a,b,c)
	n = mortar(Fill.(oneto(∞),oneto(∞)))
	k = mortar(Base.OneTo.(oneto(∞)))
	dat = BlockHcat(
	BroadcastVector((n,k,bc1,abc) -> (n + k + bc1) / (2n + abc), n, k, b+c-1, a+b+c),
	BroadcastVector((n,k,abc) -> (n + k + abc) / (2n + abc), n, k, a+b+c)
	)
	_BandedBlockBandedMatrix(dat', axes(k,1), (0,1), (0,0))
	end
	function Ry(a,b,c)
	n = mortar(Fill.(oneto(∞),oneto(∞)))
	k = mortar(Base.OneTo.(oneto(∞)))
	dat = PseudoBlockArray(Vcat(
	((k .+ (c-1) ) .* (n .+ k .+ (b+c-1)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
	((k .- n .- a ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
	((k .+ (c-1) ) .* (k .- n .- 1) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
	((n .+ k .+ (a+b+c) ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))'
	), (blockedrange(Fill(2,2)), axes(n,1)))
	_BandedBlockBandedMatrix(dat, axes(k,1), (0,1), (0,1))
	end

	function Lx(a,b,c)
	n = mortar(Fill.(oneto(∞),oneto(∞)))
	k = mortar(Base.OneTo.(oneto(∞)))
	dat = PseudoBlockArray(Vcat(
	((n .- k .+ a) ./ (2n .+ (a+b+c)))',
	((n .- k .+ 1) ./ (2n .+ (a+b+c)))'
	), (blockedrange(Fill(1,2)), axes(n,1)))
	_BandedBlockBandedMatrix(dat, axes(k,1), (1,0), (0,0))
	end

juliaapproximation / multivariateorthogonalpolynomials.jl Goto Github PK