Create blockcholesky.jl

src/BlockArrays.jl

@@ -54,11 +54,11 @@ include("views.jl")
 include("blocks.jl")
 include("blockarrayinterface.jl")
 include("blockbroadcast.jl")
+include("blockcholesky.jl")
 include("blocklinalg.jl")
 include("blockproduct.jl")
 include("show.jl")
 include("blockreduce.jl")
 include("blockdeque.jl")
-include("blockcholesky.jl")
 
 end # module

src/blockcholesky.jl

@@ -1,37 +1,33 @@
 
 
-##########################
-# Cholesky Factorization #
-##########################
+##########################################
+# Cholesky Factorization on BlockMatrices#
+##########################################
 
-"""
-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)
+    ::Val{false}=Val(false); check::Bool = false) = 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
+        for i = 1: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))
-
+            Aii, info = LinearAlgebra._chol!(Pii, UpperTriangular)
+            if !iszero(info)
+                @assert info > 0
+                if i == 1
+                    return UpperTriangular(A), info
+                end
+                info += sum(size(A[Block(l,l)])[1] for l=1:i-1) 
+                return UpperTriangular(A), info
+            end
+
             for j = i+1:n
                 Pij = getblock(A,i,j)
                 for k = 1:i-1
@@ -41,26 +37,29 @@ function b_chol!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real
             end
         end
     end
-
-    return UpperTriangular(A)
+ 
+    return UpperTriangular(A), 0
 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
+        for i = 1: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))
+            Aii, info = LinearAlgebra._chol!(Pii, LowerTriangular)
+            if !iszero(info)
+                @assert info > 0
+                if i == 1
+                    return UpperTriangular(A), info
+                end
+                info += sum(size(A[Block(l,l)])[1] for l=1:i-1) 
+                return LowerTriangular(A), info
+            end
 
             for j = i+1:n
                 Pij = getblock(A,j,i)
@@ -72,12 +71,12 @@ function b_chol!(A::BlockArray{T}, ::Type{LowerTriangular}) where T<:Real
         end
     end
 
-    return LowerTriangular(A)
+    return LowerTriangular(A), 0
 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)
+function cholesky!(A::Symmetric{<:Real,<:BlockArray}, ::Val{false}=Val(false); check::Bool = false)
+    C, info = b_chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+    #check && LinearAlgebra.checkpositivedefinite(info)
+    return Cholesky(C.data, A.uplo, info)
 end
 

test/test_cholesky.jl

@@ -9,12 +9,20 @@ using BlockArrays, Test, LinearAlgebra
     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)
+    E = BlockArray{Float64}(randn(9,9)+100I, fill(3,3), fill(3,3)); E = Symmetric(E)
+    E2 = copy(E); E2[2,2] = 0
+    E5 = copy(E); E5[5,5] = 0
+    E8 = copy(E); E8[8,8] = 0
+    nsym = BlockArray{Float64}(randn(6,8), fill(2,3), fill(2,4))
 
     A_T = Matrix(A)
     B_T = Matrix(B)
     C_T = Matrix(C)
     D_T = Matrix(D)
 
+    #Test on nonsymmetric matrix
+    @test_throws MethodError cholesky(nsym)
+
     #Tests on A
     @test cholesky(A).U ≈ cholesky(A_T).U
     @test cholesky(A).U'cholesky(A).U ≈ A
@@ -27,9 +35,14 @@ using BlockArrays, Test, LinearAlgebra
     @test cholesky(C).L ≈ cholesky(C_T).L
     @test cholesky(C).L*cholesky(C).L' ≈ C
 
-    #tests on D
+    #Tests on D
     @test cholesky(D).L ≈ cholesky(D_T).L
     @test cholesky(D).L*cholesky(D).L' ≈ D
 
+    #Tests on non-PD matrices
+    @test cholesky(E2).info == 2
+    @test cholesky(E5).info == 5
+    @test cholesky(E8).info == 8
+
 end