Implement T'U (and U'T) (#97)

* Implement T'U (and U'T)

* Add T'T

* Update test_chebyshev.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.6.2"
+version = "0.6.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -39,7 +39,7 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.3"
 InfiniteLinearAlgebra = "0.6.7"
 IntervalSets = "0.5, 0.6, 0.7"
-LazyArrays = "0.22"
+LazyArrays = "0.22.11"
 LazyBandedMatrices = "0.7.14"
 QuasiArrays = "0.9"
 SpecialFunctions = "1.0, 2"

src/classical/chebyshev.jl

@@ -148,6 +148,20 @@ massmatrix(::ChebyshevU{V}) where V = Diagonal(Fill(convert(V,π)/2,∞))
 @simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevT}}, B::ChebyshevT) = massmatrix(ChebyshevT{promote_type(eltype(A),eltype(B))}())
 @simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevU}}, B::ChebyshevU) = massmatrix(ChebyshevU{promote_type(eltype(A),eltype(B))}())
 
+@simplify function *(A::QuasiAdjoint{<:Any,<:ChebyshevT}, B::ChebyshevT)
+    T = promote_type(eltype(A), eltype(B))
+    f = (k,j) -> isodd(j-k) ? zero(T) : -(((1 + (-1)^(j + k))*(-1 + j^2 + k^2))/(j^4 + (-1 + k^2)^2 - 2j^2*(1 + k^2)))
+    BroadcastMatrix{T}(f, 0:∞, (0:∞)')
+end
+
+@simplify function *(A::QuasiAdjoint{<:Any,<:ChebyshevT}, B::ChebyshevU)
+    T = promote_type(eltype(A), eltype(B))
+    f = (k,j) -> isodd(j-k) ? zero(T) : ((one(T) + (-1)^(j + k))*(1 + j))/((1 + j - k)*(1 + j + k))
+    BroadcastMatrix{T}(f, 0:∞, (0:∞)')
+end
+
+
+
 ##########
 # Derivatives
 ##########

test/test_chebyshev.jl

@@ -456,6 +456,18 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
         T = ChebyshevT()
         @test ChebyshevT{ComplexF64}() ≡ convert(AbstractQuasiArray{ComplexF64}, T) ≡ convert(AbstractQuasiMatrix{ComplexF64}, T)
     end
+
+    @testset "innerproduct" begin
+        T = ChebyshevT()
+        U = ChebyshevU()
+        x = axes(T,1)
+
+        @test (T'U)[1:10,1:10] ≈ (U'T)[1:10,1:10]'
+        @test (T'T)[1:10,1:10] ≈ ((T'U)*(U\T))[1:10,1:10]
+        @test_broken (U'U)[1:10,1:10] ≈ ((U'T)*(T\U))[1:10,1:10]
+        f,g = (T/T\exp.(x)),(U/U\exp.(x))
+        @test f'g ≈ g'f ≈ f'f ≈ dot(f,g) ≈ dot(f,f) ≈ exp(2)/2 - exp(-2)/2
+    end
 end
 
 struct QuadraticMap{T} <: Map{T} end