Lazy MemoryLayout for QR/Cholesky Jacobi matrices (#143)

* QR and Cholesky jacobimatrix with lazy memorylayout

* add X multiplication tests

* increase version number

Author: TSGut

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.10"
+version = "0.10.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/choleskyQR.jl

@@ -56,7 +56,7 @@ function cholesky_jacobimatrix(W::AbstractMatrix, Q)
 end
 
 # The generated Jacobi operators are symmetric tridiagonal, so we store their data in two cached bands which are generated in tandem but can be accessed separately.
-mutable struct CholeskyJacobiData{T} <: AbstractMatrix{T}
+mutable struct CholeskyJacobiData{T} <: LazyMatrix{T}
     dv::AbstractVector{T} # store diagonal band entries in adaptively sized vector
     ev::AbstractVector{T} # store off-diagonal band entries in adaptively sized vector
     U::UpperTriangular{T} # store upper triangular conversion matrix (needed to extend available entries)
@@ -135,7 +135,7 @@ function qr_jacobimatrix(sqrtW::AbstractMatrix{T}, Q, method = :Q) where T
 end
 
 # The generated Jacobi operators are symmetric tridiagonal, so we store their data in two cached bands which are generated in tandem but can be accessed separately.
-mutable struct QRJacobiData{method,T} <: AbstractMatrix{T}
+mutable struct QRJacobiData{method,T} <: LazyMatrix{T}
     dv::AbstractVector{T} # store diagonal band entries in adaptively sized vector
     ev::AbstractVector{T} # store off-diagonal band entries in adaptively sized vector
     U                     # store conversion, Q method: stores QR object. R method: only stores R.

test/test_choleskyQR.jl

@@ -116,6 +116,21 @@ import LazyArrays: AbstractCachedMatrix, resizedata!
                 # Comparison with Lanczos
                 @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
             end
+
+            @testset "Jacobi matrix multiplication" begin
+                P = Normalized(legendre(0..1))
+                x = axes(P,1)
+                J = jacobimatrix(P)
+                wf(x) = (1-x)^2
+                sqrtwf(x) = (1-x)
+                # compute Jacobi matrix via decomp
+                Jchol = cholesky_jacobimatrix(wf, P)
+                JqrQ = qr_jacobimatrix(sqrtwf, P)
+                JqrR = qr_jacobimatrix(sqrtwf, P, :R)
+                @test (Jchol*Jchol)[1:10,1:10] ≈ ApplyArray(*,Jchol,Jchol)[1:10,1:10]
+                @test (Jchol*Jchol)[1:10,1:10] ≈ (JqrQ*JqrQ)[1:10,1:10]
+                @test (Jchol*Jchol)[1:10,1:10] ≈ (JqrR*JqrR)[1:10,1:10]
+            end
         end
     end