Cholesky (#63)

* Add cholesky to factorizations

* Update factorizations.jl

* Update factorizations.jl

* add tests, cholesky default

* Update test_layoutarray.jl

Author: dlfivefifty

src/ArrayLayouts.jl

@@ -35,7 +35,7 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcas
 
 import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, factorize, qr, lu, cholesky,
                         norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!, AdjOrTrans, HermOrSym, copy_oftype,
-                        AdjointAbsVec, TransposeAbsVec
+                        AdjointAbsVec, TransposeAbsVec, cholcopy
 
 import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
 

src/factorizations.jl

@@ -281,10 +281,13 @@ end
 __qr(layout, lengths, A; kwds...) = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
 _qr(layout, axes, A; kwds...) = __qr(layout, map(length, axes), A; kwds...)
 _qr(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(qr, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
+_qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
 _lu(layout, axes, A; kwds...) = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)}}, A; kwds...)
 _lu(layout, axes, A, pivot::P; kwds...) where P = Base.invoke(lu, Tuple{AbstractMatrix{eltype(A)},P}, A, pivot; kwds...)
-_qr!(layout, axes, A, args...; kwds...) = error("Overload _qr!(::$(typeof(layout)), axes, A)")
 _lu!(layout, axes, A, args...; kwds...) = error("Overload _lu!(::$(typeof(layout)), axes, A)")
+_cholesky(layout, axes, A, ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
+_cholesky(layout, axes, A, ::Val{true}; tol = 0.0, check::Bool = true) = cholesky!(cholcopy(A), Val(true); tol = tol, check = check)
+_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
 
 _inv_eye(_, ::Type{T}, axs::NTuple{2,Base.OneTo{Int}}) where T = Matrix{T}(I, map(length,axs)...)
@@ -307,6 +310,8 @@ end
 
 macro _layoutfactorizations(Typ)
     esc(quote
+        LinearAlgebra.cholesky(A::$Typ, args...; kwds...) = ArrayLayouts._cholesky(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
+        LinearAlgebra.cholesky!(A::$Typ, v::Val{false}=Val(false); check::Bool = true) = ArrayLayouts._cholesky!(ArrayLayouts.MemoryLayout(A), axes(A), A, v; check=check)
         LinearAlgebra.qr(A::$Typ, args...; kwds...) = ArrayLayouts._qr(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
         LinearAlgebra.qr!(A::$Typ, args...; kwds...) = ArrayLayouts._qr!(ArrayLayouts.MemoryLayout(A), axes(A), A, args...; kwds...)
         LinearAlgebra.lu(A::$Typ, pivot::Union{Val{false}, Val{true}}; kwds...) = ArrayLayouts._lu(ArrayLayouts.MemoryLayout(A), axes(A), A, pivot; kwds...)
@@ -321,5 +326,17 @@ macro layoutfactorizations(Typ)
     esc(quote
         ArrayLayouts.@_layoutfactorizations $Typ
         ArrayLayouts.@_layoutfactorizations SubArray{<:Any,2,<:$Typ}
+        ArrayLayouts.@_layoutfactorizations LinearAlgebra.RealHermSymComplexHerm{<:Real,<:$Typ}
+        ArrayLayouts.@_layoutfactorizations LinearAlgebra.RealHermSymComplexHerm{<:Real,<:SubArray{<:Real,2,<:$Typ}}
     end)
-end
+end
+
+function ldiv!(C::Cholesky{<:Any,<:AbstractMatrix}, B::LayoutArray)
+    if C.uplo == 'L'
+        return ldiv!(adjoint(LowerTriangular(C.factors)), ldiv!(LowerTriangular(C.factors), B))
+    else
+        return ldiv!(UpperTriangular(C.factors), ldiv!(adjoint(UpperTriangular(C.factors)), B))
+    end
+end
+
+

test/test_layoutarray.jl

@@ -81,16 +81,21 @@ MemoryLayout(::Type{MyVector}) = DenseColumnMajor()
 
         @test copyto!(view(MyMatrix(Array{Float64}(undef,5,5)),:,:), view(A',:,:)) == A'
 
-        @test qr(A).factors ≈ qr(A.A).factors
-        @test qr(A,Val(true)).factors ≈ qr(A.A,Val(true)).factors
-        @test lu(A).factors ≈ lu(A.A).factors
-        @test lu(A,Val(true)).factors ≈ lu(A.A,Val(true)).factors
-        @test_throws ErrorException qr!(A)
-        @test_throws ErrorException lu!(A)
-
-        @test qr(A) isa LinearAlgebra.QRCompactWY
-        @test inv(A) ≈ inv(A.A)
-
+        @testset "factorizations" begin
+            @test qr(A).factors ≈ qr(A.A).factors
+            @test qr(A,Val(true)).factors ≈ qr(A.A,Val(true)).factors
+            @test lu(A).factors ≈ lu(A.A).factors
+            @test lu(A,Val(true)).factors ≈ lu(A.A,Val(true)).factors
+            @test_throws ErrorException qr!(A)
+            @test_throws ErrorException lu!(A)
+
+            @test qr(A) isa LinearAlgebra.QRCompactWY
+            @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,Val(true)).U ≈ cholesky(Matrix(S),Val(true)).U
+        end
         Bin = randn(5,5)
         B = MyMatrix(copy(Bin))
         muladd!(1.0, A, A, 2.0, B)