Workaround inference of Cholesky factorisations of BandedMatrices (#67)

* Copy cholesky! definition from LinearAlgebra.jl

The _cholesky! definition using `invoke` is left as a fallback.

* Add cholesky inference test using BandedMatrix

* Don't use BandedMatrices in tests

* Bump patch version

Author: jipolanco

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.7.7"
+version = "0.7.8"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/factorizations.jl

@@ -290,6 +290,13 @@ _cholesky(layout, axes, A, ::Val{true}; tol = 0.0, check::Bool = true) = cholesk
 _cholesky!(layout, axes, A, v::Val{tf}; kwds...) where tf = Base.invoke(cholesky!, Tuple{LinearAlgebra.RealHermSymComplexHerm,Val{tf}}, A, v; kwds...)
 _factorize(layout, axes, A) = qr(A) # Default to QR
 
+# Cholesky factorization without pivoting (copied from stdlib/LinearAlgebra).
+function _cholesky!(layout, axes, A::LinearAlgebra.RealHermSymComplexHerm, ::Val{false}; check::Bool = true)
+    C, info = LinearAlgebra._chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+    check && LinearAlgebra.checkpositivedefinite(info)
+    return Cholesky(C.data, A.uplo, info)
+end
+
 _inv_eye(_, ::Type{T}, axs::NTuple{2,Base.OneTo{Int}}) where T = Matrix{T}(I, map(length,axs)...)
 function _inv_eye(A, ::Type{T}, (rows,cols)) where T
     dest = zero!(similar(A, T, (cols,rows)))

test/test_layoutarray.jl

@@ -99,7 +99,7 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
             @test inv(A) ≈ inv(A.A)
 
             S = Symmetric(MyMatrix(reshape(inv.(1:25),5,5) + 10I))
-            @test cholesky(S).U ≈ cholesky!(deepcopy(S)).U
+            @test cholesky(S).U ≈ @inferred(cholesky!(deepcopy(S))).U
             @test cholesky(S,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
         end
         Bin = randn(5,5)
@@ -349,4 +349,4 @@ triangulardata(A::MyUpperTriangular) = triangulardata(A.A)
     @test lmul!(U, copy(B)) ≈ U * B
 
     @test_skip lmul!(U,view(copy(B),collect(1:5),1:5)) ≈ U * B
-end
+end