Add lowering operators

Author: dlfivefifty

src/QuasiArrays/matmul.jl

@@ -32,7 +32,7 @@ end
 QuasiMatMulMat{styleA, styleB, T, V} = QuasiArrayMulArray{styleA, styleB, 2, 2, T, V}
 QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styleB, 2, 2, T, V}
 
-*(A::AbstractQuasiArray, B::AbstractQuasiArray, C::AbstractQuasiArray) = materialize(Mul(A,B,C))
+*(A::AbstractQuasiArray, B::AbstractQuasiArray, C::AbstractQuasiArray, D::AbstractQuasiArray...) = materialize(Mul(A,B,C,D...))
 *(A::AbstractQuasiArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
 *(A::AbstractQuasiArray, B::AbstractArray) = materialize(Mul(A,B))
 *(A::AbstractArray, B::AbstractQuasiArray) = materialize(Mul(A,B))

src/bases/jacobi.jl

@@ -26,18 +26,85 @@ axes(::AbstractJacobi) = (ChebyshevInterval(), OneTo(∞))
 materialize(M::Mul2{<:Any,<:Any,<:QuasiAdjoint{<:Any,<:Legendre},<:Legendre}) =
     Diagonal(2 ./ (2(0:∞) .+ 1))
 
+# pinv(Jacobi(b+1,a+1))D*W*Jacobi(a,b)
+function materialize(M::Mul{<:Tuple,<:Tuple{<:PInv{<:Any,<:Jacobi},
+                                        <:Derivative{<:Any,<:ChebyshevInterval},
+                                        <:Jacobi}})
+    Ji, D, S = M.factors
+    J = parent(Ji)
+    (J.b == S.b+1 && J.a == S.a+1) || throw(ArgumentError())
+    _BandedMatrix(((1:∞ .+ 1 .+ S.a .+ S.b)/2)', ∞, -1,1)
+end
+
 
 function materialize(M::Mul2{<:Any,<:Any,<:Derivative{<:Any,<:ChebyshevInterval},<:Jacobi})
-    _, S = M.factors
-    D = _BandedMatrix(((1:∞ .+ 1 .+ S.a .+ S.b)/2)', ∞, -1,1)
-    Mul(Jacobi(S.a+1,S.b+1), D)
+    D, S = M.factors
+    A = PInv(Jacobi(S.b+1,S.a+1))*D*S
+    Mul(Jacobi(S.b+1,S.a+1), A)
+end
+
+# pinv(Legendre())D*W*Jacobi(true,true)
+function materialize(M::Mul{<:Tuple,<:Tuple{<:PInv{<:Any,<:Legendre},
+                                        <:Derivative{<:Any,<:ChebyshevInterval},
+                                        QuasiDiagonal{Bool,JacobiWeight{Bool}},Jacobi{Bool}}})
+    Li, _, W, S = M.factors
+    w = parent(W)
+    (w.a && S.a && w.b && S.b) || throw(ArgumentError())
+    _BandedMatrix((-2*(1:∞))', ∞, 1,-1)
 end
 
+# reduce to Legendre
 function materialize(M::Mul{<:Tuple,<:Tuple{<:Derivative{<:Any,<:ChebyshevInterval},
                                         QuasiDiagonal{Bool,JacobiWeight{Bool}},Jacobi{Bool}}})
-    _, W, S = M.factors
+    D, W, S = M.factors
     w = parent(W)
     (w.a && S.a && w.b && S.b) || throw(ArgumentError())
-    D = _BandedMatrix((-2*(1:∞))', ∞, 1,-1)
-    Mul(Legendre(), D)
+    A = pinv(Legendre{eltype(M)}())*D*W*S
+    Mul(Legendre(), A)
+end
+
+function materialize(M::Mul{<:Tuple,<:Tuple{<:PInv{<:Any,<:Jacobi{Bool}},
+                            QuasiDiagonal{Bool,JacobiWeight{Bool}},Jacobi{Bool}}})
+    Ji, W, S = M.factors
+    J,w = parent(Ji),parent(W)
+    @assert  S.b && S.a
+    if w.b && !w.a
+        @assert !J.b && J.a
+        _BandedMatrix(Vcat(((2:2:∞)./(3:2:∞))',((2:2:∞)./(3:2:∞))'), ∞, 1,0)
+    elseif !w.b && w.a
+        @assert J.b && !J.a
+        _BandedMatrix(Vcat(((2:2:∞)./(3:2:∞))',(-(2:2:∞)./(3:2:∞))'), ∞, 1,0)
+    else
+        error("Not implemented")
+    end
+end
+
+function materialize(M::Mul{<:Tuple,<:Tuple{<:PInv{<:Any,<:Legendre},
+                            QuasiDiagonal{Bool,JacobiWeight{Bool}},Jacobi{Bool}}})
+    Li, W, S = M.factors
+    L,w = parent(Li),parent(W)
+    if w.b && w.a
+        @assert S.b && S.a
+        _BandedMatrix(Vcat(((2:2:∞)./(3:2:∞))', Zeros(1,∞), (-(2:2:∞)./(3:2:∞))'), ∞, 2,0)
+    elseif w.b && !w.a
+        @assert S.b && !S.a
+        _BandedMatrix(Ones{eltype(M)}(2,∞), ∞, 1,0)
+    elseif !w.b && w.a
+        @assert !S.b && S.a
+        _BandedMatrix(Vcat(Ones{eltype(M)}(1,∞),-Ones{eltype(M)}(1,∞)), ∞, 1,0)
+    else
+        error("Not implemented")
+    end
+end
+
+function materialize(M::Mul{<:Tuple,<:Tuple{QuasiAdjoint{Bool,Jacobi{Bool}},
+                                        QuasiDiagonal{Int,JacobiWeight{Int}},Jacobi{Bool}}})
+    St, W, S = M.factors
+
+    w = parent(W)
+    (w.b == 2 && S.b && w.a == 2 && S.a && parent(St) == S) || throw(ArgumentError())
+    W_sqrt = Diagonal(JacobiWeight(true,true))
+    L = Legendre()
+    A = PInv(L)*W_sqrt*S
+    A'*(L'L)*A
 end

test/runtests.jl

@@ -158,6 +158,82 @@ S = Jacobi(true,true)
 W = Diagonal(JacobiWeight(true,true))
 D = Derivative(axes(W,1))
 
+A = @inferred(PInv(Jacobi(2,2))*D*S)
+@test A isa BandedMatrix
+@test size(A) == (∞,∞)
+@test A[1:10,1:10] == diagm(1 => 1:0.5:5)
+
+M = @inferred(D*S)
+@test M isa Mul
+@test M.factors[1] == Jacobi(2,2)
+@test M.factors[2][1:10,1:10] == A[1:10,1:10]
+
+
+L = Diagonal(JacobiWeight(true,false))
+A = @inferred(PInv(Jacobi(false,true))*L*S)
+@test A isa BandedMatrix
+@test size(A) == (∞,∞)
+
+L = Diagonal(JacobiWeight(false,true))
+A = @inferred(PInv(Jacobi(true,false))*L*S)
+@test A isa BandedMatrix
+@test size(A) == (∞,∞)
+
+
+L = Diagonal(JacobiWeight(true,true))
+
+
+L = Legendre()
+
+A, B = (L'L),(PInv(L)*W*S)
+
+M = A*B
+    @which M.mul[1,1]
+
+@which materialize(Mul(A,B))
+
+M = LazyArrays.MulArray(Mul(A,B))
+    axes(M)
+
+@time A[1:100000,1:100000]
+@profiler A.data[:,1:10_000]
+
+V = view(A.data,:,1:100)
+
+N = 10_000; M = randn(3,N)
+
+A.data.arrays[2]
+
+@time begin
+    M[1,:] .= view(A.data.arrays[1]', 1:N)
+    M[2,:] .= view(A.data.arrays[2], 1, 1:N)
+    M[3,:] .= view(A.data.arrays[3]', 1:N)
+end
+M
+
+@which copyto!(M, V)
+
+@time
+
+typeof(V)
+
+using Profile
+Profile.clear()
+@time randn(3,10_000)
+
+W*S
+
+
+
+
+M = Mul(S',W',W,S)
+materialize(M)
+materialize(Mul(S',W',W))
+
+@which materialize(M)
+
+W*W
+
 S'W'W*S
 
 N = 10
@@ -192,7 +268,7 @@ A.*B
 _Vcat(broadcast((a,b) -> broadcast(op,a,b), A_arrays, B.arrays))
 
 f = Legendre() * Vcat(randn(20), Zeros(∞))
-@time L'f
+@time Vector(L'f)
 
 A,v=(L'f).factors[end-1:end]