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,
@@ -54,10 +54,12 @@ 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("blockbandedcholesky.jl")
 
 end # module

src/blockbandedcholesky.jl

@@ -0,0 +1,91 @@
+
+
+sparsecholesky(A::Symmetric{<:Real,<:BlockArray},
+::Val{false}=Val(false); check::Bool = true) = sparsecholesky!(cholcopy(A); check = check)
+
+
+
+function _blockbandedcholesky!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real
+    n = blocksize(A)[1]
+    k_end = 0
+
+    @inbounds begin
+        for i = 1:n
+            Pii = getblock(A,i,i)
+            for k = 1:i-1
+                muladd!(-one(T), getblock(A,k,i)', getblock(A,k,i), one(T), Pii)
+            end
+            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
+
+            k_start = n
+            while getblock(A,i,k_start) ≈ zeros(Float32, size(getblock(A,i,k_start))) && k_start > k_end
+                k_start -= 1
+            end
+            k_end = k_start
+
+            for j = k_start:-1:i+1
+                Pij = getblock(A,i,j)
+                for k = 1:i-1
+                    muladd!(-one(T), getblock(A,k,i)', getblock(A,k,j), one(T), Pij)
+                end
+                ldiv!(UpperTriangular(getblock(A,i,i))', Pij)
+            end
+        end
+    end
+
+    return UpperTriangular(A), 0
+end
+
+function _blockbandedcholesky!(A::BlockArray{T}, ::Type{LowerTriangular}) where T<:Real
+    n = blocksize(A)[1]
+    k_end = 0
+
+    @inbounds begin
+        for i = 1:n
+            Pii = getblock(A,i,i) 
+            for k = 1:i-1
+                muladd!(-one(T), getblock(A,i,k), getblock(A,i,k)', one(T), Pii)
+            end
+            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
+
+            k_start = n
+            while getblock(A,k_start,i) ≈ zeros(Float32, size(getblock(A,k_start,i))) && k_start > k_end
+                k_start -= 1
+            end
+            k_end = k_start
+
+            for j = k_start:-1:i+1
+                Pij = getblock(A,j,i)
+                for k = 1:i-1
+                    muladd!(-one(T), getblock(A,j,k), getblock(A,i,k)', one(T), Pij)
+                end
+                rdiv!(Pij, LowerTriangular(getblock(A,i,i))')
+            end
+        end
+    end
+
+    return LowerTriangular(A), 0
+end
+
+
+function sparsecholesky!(A::Symmetric{<:Real,<:BlockArray}, ::Val{false}=Val(false); check::Bool = true)
+    C, info = _blockbandedcholesky!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+    #check && LinearAlgebra.checkpositivedefinite(info)
+    return Cholesky(C.data, A.uplo, info)
+end

src/blockcholesky.jl

@@ -0,0 +1,82 @@
+
+
+##########################################
+# Cholesky Factorization on BlockMatrices#
+##########################################
+
+
+cholesky(A::Symmetric{<:Real,<:BlockArray},
+    ::Val{false}=Val(false); check::Bool = false) = cholesky!(cholcopy(A); check = check)
+
+
+function _block_chol!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real
+    n = blocksize(A)[1]
+
+    @inbounds begin
+        for i = 1:n
+            Pii = getblock(A,i,i) 
+            for k = 1:i-1
+                muladd!(-one(T), getblock(A,k,i)', getblock(A,k,i), one(T), Pii)
+            end
+            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
+                    muladd!(-one(T), getblock(A,k,i)', getblock(A,k,j), one(T), Pij)
+                end
+                ldiv!(UpperTriangular(getblock(A,i,i))', Pij)
+            end
+        end
+    end
+
+    return UpperTriangular(A), 0
+end
+
+
+function _block_chol!(A::BlockArray{T}, ::Type{LowerTriangular}) where T<:Real
+    n = blocksize(A)[1]
+
+    @inbounds begin
+        for i = 1:n
+            Pii = getblock(A,i,i) 
+            for k = 1:i-1
+                muladd!(-one(T), getblock(A,i,k), getblock(A,i,k)', one(T), Pii)
+            end
+            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)
+                for k = 1:i-1
+                    muladd!(-one(T), getblock(A,j,k), getblock(A,i,k)', one(T), Pij)
+                end
+                rdiv!(Pij, LowerTriangular(getblock(A,i,i))')
+            end
+        end
+    end
+
+    return LowerTriangular(A), 0
+end
+
+function cholesky!(A::Symmetric{<:Real,<:BlockArray}, ::Val{false}=Val(false); check::Bool = false)
+    C, info = _block_chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+    #check && LinearAlgebra.checkpositivedefinite(info)
+    return Cholesky(C.data, A.uplo, info)
+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,48 @@
+using BlockArrays, Test, LinearAlgebra
+
+
+
+@testset "Block cholesky" begin
+
+    # Generating random positive definite and symmetric matrices
+    A = BlockArray{Float32}(randn(9,9)+100I, fill(3,3), fill(3,3)); A = Symmetric(A)
+    B = BlockArray{Float32}(randn(55,55)+100I, 1:10, 1:10); B = Symmetric(B)
+    C = BlockArray{Float32}(randn(9,9)+100I, fill(3,3), fill(3,3)); C = Symmetric(C, :L)
+    D = BlockArray{Float32}(randn(55,55)+100I, 1:10, 1:10); D = Symmetric(D, :L)
+    E = BlockArray{Float32}(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{Float32}(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
+
+    #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
+
+    #Tests on non-PD matrices
+    @test cholesky(E2).info == 2
+    @test cholesky(E5).info == 5
+    @test cholesky(E8).info == 8
+
+end
+