Create blockcholesky.jl

src/BlockArrays.jl

@@ -35,7 +35,7 @@ using Base: ReshapedArray, dataids
 import Base: (:), IteratorSize, iterate, axes1, strides, isempty
 import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted, broadcastable
 import LinearAlgebra: lmul!, rmul!, AbstractTriangular, HermOrSym, AdjOrTrans,
-                        StructuredMatrixStyle
+                        StructuredMatrixStyle, cholesky, cholesky!, cholcopy
 import ArrayLayouts: _fill_lmul!, MatMulVecAdd, MatMulMatAdd, MatLmulVec, MatLdivVec,
                         materialize!, MemoryLayout, sublayout, transposelayout, conjlayout, 
                         triangularlayout, triangulardata, _inv, _copyto!, axes_print_matrix_row,
@@ -59,5 +59,6 @@ include("blockproduct.jl")
 include("show.jl")
 include("blockreduce.jl")
 include("blockdeque.jl")
+include("blockcholesky.jl")
 
 end # module

src/blockcholesky.jl

@@ -0,0 +1,83 @@
+
+
+##########################
+# Cholesky Factorization #
+##########################
+
+"""
+Funtions to do
+2.check square
+3.check positive definite
+"""
+
+cholesky(A::Symmetric{<:Real,<:BlockArray},
+    ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
+
+
+function b_chol!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real
+    n = blocksize(A)[1]
+
+    @inbounds begin
+        # Initializing the first role of blocks
+        cholesky!(Symmetric(getblock(A,1,1)))
+        for j = 2:n
+            ldiv!(UpperTriangular(getblock(A,1,1))', getblock(A,1,j))
+        end
+
+        # For the left blocks
+        for i = 2:n
+            Pii = getblock(A,i,i) 
+            for k = 1:i-1
+                muladd!(-1.0, getblock(A,k,i)', getblock(A,k,i), 1.0, Pii)
+            end
+            cholesky!(Symmetric(Pii))
+
+            for j = i+1:n
+                Pij = getblock(A,i,j)
+                for k = 1:i-1
+                    muladd!(-1.0, getblock(A,k,i)', getblock(A,k,j), 1.0, Pij)
+                end
+                ldiv!(UpperTriangular(getblock(A,i,i))', Pij)
+            end
+        end
+    end
+
+    return UpperTriangular(A)
+end
+function b_chol!(A::BlockArray{T}, ::Type{LowerTriangular}) where T<:Real
+    n = blocksize(A)[1]
+
+    @inbounds begin
+        # Initializing the first role of blocks
+        cholesky!(Symmetric(getblock(A,1,1), :L))
+        for j = 2:n
+            rdiv!(getblock(A,j,1), LowerTriangular(getblock(A,1,1))')
+        end
+
+        # For the left blocks
+        for i = 2:n
+            Pii = getblock(A,i,i) 
+            for k = 1:i-1
+                muladd!(-1.0, getblock(A,i,k), getblock(A,i,k)', 1.0, Pii)
+            end
+            cholesky!(Symmetric(Pii, :L))
+
+            for j = i+1:n
+                Pij = getblock(A,j,i)
+                for k = 1:i-1
+                    muladd!(-1.0, getblock(A,j,k), getblock(A,i,k)', 1.0, Pij)
+                end
+                rdiv!(Pij, LowerTriangular(getblock(A,i,i))')
+            end
+        end
+    end
+
+    return LowerTriangular(A)
+end
+
+function cholesky!(A::Symmetric{<:Real,<:BlockArray}, ::Val{false}=Val(false); check::Bool = true)
+    C = b_chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+    #check && checkpositivedefinite(info)
+    return Cholesky(C.data, A.uplo, 0)
+end
+

test/runtests.jl

@@ -12,4 +12,5 @@ using BlockArrays, LinearAlgebra, Test
     include("test_blockproduct.jl")
     include("test_blockreduce.jl")
     include("test_blockdeque.jl")
+    include("test_cholesky.jl")
 end

test/test_cholesky.jl

@@ -0,0 +1,35 @@
+using BlockArrays, Test, LinearAlgebra
+
+
+
+@testset "Block cholesky" begin
+
+    # Generating random positive definite and symmetric matrices
+    A = BlockArray{Float64}(randn(9,9)+100I, fill(3,3), fill(3,3)); A = Symmetric(A)
+    B = BlockArray{Float64}(randn(55,55)+100I, 1:10, 1:10); B = Symmetric(B)
+    C = BlockArray{Float64}(randn(9,9)+100I, fill(3,3), fill(3,3)); C = Symmetric(C, :L)
+    D = BlockArray{Float64}(randn(55,55)+100I, 1:10, 1:10); D = Symmetric(D, :L)
+
+    A_T = Matrix(A)
+    B_T = Matrix(B)
+    C_T = Matrix(C)
+    D_T = Matrix(D)
+
+    #Tests on A
+    @test cholesky(A).U ≈ cholesky(A_T).U
+    @test cholesky(A).U'cholesky(A).U ≈ A
+
+    #Tests on B
+    @test cholesky(B).U ≈ cholesky(B_T).U
+    @test cholesky(B).U'cholesky(B).U ≈ B
+
+    #Tests on C
+    @test cholesky(C).L ≈ cholesky(C_T).L
+    @test cholesky(C).L*cholesky(C).L' ≈ C
+
+    #tests on D
+    @test cholesky(D).L ≈ cholesky(D_T).L
+    @test cholesky(D).L*cholesky(D).L' ≈ D
+
+end
+