Create blockcholesky.jl (#126)

* Create blockcholesky.jl

* Add files via upload

* overload cholesk!

* Including the test file in REPL will do the test

* Including test_cholesky.jl in runtests.jl

* Added lower Cholesky factorization

* Added PD check

* Changed the function name and types

* Added blockbanded case

* Rearrange

* Capable on both BlockArray and BlockBandedArray

* Changes

* Update Choleksy

* tests pass

* Support sparse block structure

* Update blockcholesky.jl

* Fix tests on 1.8, use ArrayLayouts.ldiv!

* Speed up BlockKron construction

* Update test_blockproduct.jl

* increase cov

Co-authored-by: Sheehan Olver <[email protected]>

Author: felixw98

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockArrays"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.16.17"
+version = "0.16.18"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,6 +18,7 @@ OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [targets]
-test = ["Base64", "OffsetArrays", "SparseArrays", "StaticArrays", "Test"]
+test = ["Base64", "OffsetArrays", "SparseArrays", "StaticArrays", "Test", "Random"]

src/BlockArrays.jl

@@ -22,12 +22,12 @@ import Base: @propagate_inbounds, Array, to_indices, to_index,
             unsafe_indices, first, last, size, length, unsafe_length,
             unsafe_convert,
             getindex, ndims, show,
-            step, 
+            step,
             broadcast, eltype, convert, similar,
             @_inline_meta, _maybetail, tail, @_propagate_inbounds_meta, reindex,
             RangeIndex, Int, Integer, Number,
             +, -, *, /, \, min, max, isless, in, copy, copyto!, axes, @deprecate,
-            BroadcastStyle, checkbounds, throw_boundserror, 
+            BroadcastStyle, checkbounds, throw_boundserror,
             ones, zeros, intersect, Slice, resize!
 using Base: ReshapedArray, dataids
 import Base: AbstractArray
@@ -36,13 +36,13 @@ import Base: AbstractArray
 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, RealHermSymComplexHerm
 import ArrayLayouts: _fill_lmul!, MatMulVecAdd, MatMulMatAdd, MatLmulVec, MatLdivVec,
-                        materialize!, MemoryLayout, sublayout, transposelayout, conjlayout, 
+                        materialize!, MemoryLayout, sublayout, transposelayout, conjlayout,
                         triangularlayout, triangulardata, _inv, _copyto!, axes_print_matrix_row,
                         colsupport, rowsupport, sub_materialize
 
-include("blockindices.jl")                        
+include("blockindices.jl")
 include("blockaxis.jl")
 include("abstractblockarray.jl")
 include("blockarray.jl")
@@ -51,6 +51,7 @@ include("views.jl")
 include("blocks.jl")
 
 include("blockbroadcast.jl")
+include("blockcholesky.jl")
 include("blocklinalg.jl")
 include("blockproduct.jl")
 include("show.jl")
@@ -63,4 +64,5 @@ include("blockarrayinterface.jl")
 @deprecate setblock!(A::AbstractBlockArray{T,N}, v, I::Vararg{Integer, N}) where {T,N} (A[Block(I...)] = v)
 
 
+
 end # module

src/blockcholesky.jl

@@ -0,0 +1,89 @@
+
+##########################################
+# Cholesky Factorization on BlockMatrices#
+##########################################
+
+
+function _diag_chol!(A::AbstractArray{T}, i::Int, uplo) where T<:Real
+    Pii = view(A,Block(i,i))
+    for k = Int(first(blockcolsupport(A,Block(i)))):i-1
+        muladd!(-one(T), view(A,Block(k,i))', view(A,Block(k,i)), one(T), Pii)
+    end
+    return LAPACK.potrf!('U', Pii)
+end
+
+
+function _nondiag_chol!(A::AbstractArray{T}, i::Int, n::Int, ::Type{UpperTriangular}) where T<:Real
+    for j = i+1:Int(last(blockrowsupport(A,Block(i))))
+        Pij = view(A,Block(i,j))
+        for k = Int(first(blockcolsupport(A,Block(j)))):i-1
+            muladd!(-one(T), view(A,Block(k,i))', view(A,Block(k,j)), one(T), Pij)
+        end
+        # use ArrayLayouts.ldiv! to take advantage of layout
+        ArrayLayouts.ldiv!(transpose(UpperTriangular(view(A,Block(i,i)))), Pij)
+    end
+end
+
+function _nondiag_chol!(A::AbstractArray{T}, i::Int, n::Int, ::Type{LowerTriangular}) where T<:Real
+    for j = i+1:Int(last(blockrowsupport(A,Block(i))))
+        Pij = view(A,Block(i,j))
+        for k = Int(first(blockcolsupport(A,Block(j)))):i-1
+            muladd!(-one(T), view(A,Block(k,i))', view(A,Block(k,j)), one(T), Pij)
+        end
+        ArrayLayouts.ldiv!(UpperTriangular(view(A,Block(i,i)))', Pij)
+    end
+end
+
+function _block_chol!(A::AbstractArray{T}, ::Type{UpperTriangular}) where T<:Real
+    n = blocksize(A)[1]
+
+    @inbounds begin
+        for i = 1:n
+            _, info = _diag_chol!(A, i, 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
+
+            _nondiag_chol!(A, i, n, UpperTriangular)
+        end
+    end
+
+    return UpperTriangular(A), 0
+end
+
+function _block_chol!(A::AbstractArray{T}, ::Type{LowerTriangular}) where T<:Real
+    n = blocksize(A)[1]
+    A = copy(transpose(A))
+
+    @inbounds begin
+        for i = 1:n
+            _, info = _diag_chol!(A, i, LowerTriangular)
+
+            if !iszero(info)
+                @assert info > 0
+                if i == 1
+                    return LowerTriangular(copy(transpose(A))), info
+                end
+                info += sum(size(A[Block(l,l)])[1] for l=1:i-1) 
+                return LowerTriangular(A), info
+            end
+    
+            _nondiag_chol!(A, i, n, LowerTriangular)
+        end
+    end
+
+    return LowerTriangular(transpose(A)), 0
+end
+
+function ArrayLayouts._cholesky!(layout, ::NTuple{2,BlockedUnitRange}, A::RealHermSymComplexHerm, ::ArrayLayouts.CNoPivot; check::Bool = true)
+    C, info = _block_chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
+    check && LinearAlgebra.checkpositivedefinite(info)
+    return Cholesky(C.data, A.uplo, info)
+end
+

src/blocklinalg.jl

@@ -138,7 +138,7 @@ function _copyto!(_, ::AbstractBlockLayout, dest::AbstractMatrix, src::AbstractM
         return dest
     end
 
-    @inbounds for J = blockaxes(src,2),K = blockaxes(src,1)
+    @inbounds for J = blockaxes(src,2),K = blockcolsupport(src,J)
         copyto!(view(dest,K,J), view(src,K,J))
     end
     dest

src/blockproduct.jl

@@ -89,6 +89,12 @@ kron_viewblock(args, k::Integer, j::Integer) = args[1][k,j]*BlockKron(tail(args)
 viewblock(K::BlockKron{<:Any,1}, k::Block{1}) = kron_viewblock(K.args, Int(k))
 viewblock(K::BlockKron{<:Any,2}, kj::Block{2}) = kron_viewblock(K.args, Int.(kj.n)...)
 
+# At the moment this only takes advantage of AbstractBlockLayout structure
+struct KronLayout <: AbstractBlockLayout end
+MemoryLayout(::Type{<:BlockKron}) = KronLayout()
+blockcolsupport(K::BlockKron, J) = Block.(colsupport(first(K.args), Int.(J)))
+blockrowsupport(K::BlockKron, J) = Block.(rowsupport(first(K.args), Int.(J)))
+
 # const SubKron{T,M1,M2,R1,R2} = SubArray{T,2,<:BlockKron{T,M1,M2},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
 
 

test/runtests.jl

@@ -1,6 +1,7 @@
 using BlockArrays, LinearAlgebra, Test
 
 
+
 include("test_blockindices.jl")
 include("test_blockarrays.jl")
 include("test_blockviews.jl")
@@ -12,3 +13,4 @@ include("test_blocklinalg.jl")
 include("test_blockproduct.jl")
 include("test_blockreduce.jl")
 include("test_blockdeque.jl")
+include("test_blockcholesky.jl")

test/test_blockcholesky.jl

@@ -0,0 +1,62 @@
+using BlockArrays, Test, LinearAlgebra, Random
+
+
+Random.seed!(0)
+
+@testset "Block cholesky" begin
+
+    # Generating random positive definite and symmetric matrices
+    A = BlockArray{Float32}(rand(9,9)+100I, fill(3,3), fill(3,3)); A = Symmetric(A)
+    B = BlockArray{Float32}(rand(55,55)+100I, 1:10, 1:10); B = Symmetric(B)
+    C = BlockArray{Float32}(rand(9,9)+100I, fill(3,3), fill(3,3)); C = Symmetric(C, :L)
+    D = BlockArray{Float32}(rand(55,55)+100I, 1:10, 1:10); D = Symmetric(D, :L)
+    E = BlockArray{Float32}(rand(9,9)+100I, fill(3,3), fill(3,3)); E = Symmetric(E)
+    D1 = copy(D); D1[1,1] = 0
+    D2 = copy(D); D2[2,2] = 0
+    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
+    if VERSION < v"1.8-"
+        @test_throws MethodError cholesky(nsym)
+    else
+        @test_throws DimensionMismatch cholesky(nsym)
+    end
+
+    #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_throws PosDefException cholesky(D1)
+    @test_throws PosDefException cholesky(D2)
+    @test_throws PosDefException cholesky(E2)
+    @test_throws PosDefException cholesky(E5)
+    @test_throws PosDefException cholesky(E8)
+
+    @test cholesky(D1; check=false).info == 1
+    @test cholesky(D2; check=false).info == 2
+    @test cholesky(E2; check=false).info == 2
+    @test cholesky(E5; check=false).info == 5
+    @test cholesky(E8; check=false).info == 8
+end
+

test/test_blockproduct.jl

@@ -124,18 +124,20 @@ using BlockArrays, Test
         K̃ = BlockKron(A,B)
         @test K == K̃ == kron(A,B)
         @test blocksize(K) == blocksize(K̃) == size(A)
+        @test BlockArrays.blockcolsupport(K, Block(1)) == BlockArrays.blockcolsupport(K̃, Block(1)) == Block.(BlockArrays.colsupport(A,1))
+        @test BlockArrays.blockrowsupport(K, Block(1)) == BlockArrays.blockrowsupport(K̃, Block(1)) == Block.(BlockArrays.rowsupport(A,1))
         @test K[Block(1,1)] == K̃[Block(1),Block(1)] == A[1,1]*B
         @test K[Block(2,3)] == K̃[Block(2),Block(3)] == A[2,3]*B
         C = randn(2,5)
         K̄ = BlockKron(A,B,C)
-        @test K̄ == blockkron(A,B,C) == kron(A,B,C)
+        @test K̄ ≈ blockkron(A,B,C) ≈ kron(A,B,C)
         @test K̄[Block(1,1)][Block(1,1)] ≈ A[1,1]*B[1,1]*C
         @test K̄[Block(2,3)][Block(3,4)] ≈ A[2,3]*B[3,4]*C
 
         @test blockkron(a,B) == kron(a,B)
         @test blockkron(A,b) == kron(A,b)
-        @test blockkron(A,b,c) == kron(A,b,c)
-        @test blockkron(A,b,C) == kron(A,b,C)
+        @test blockkron(A,b,c) ≈ kron(A,b,c)
+        @test blockkron(A,b,C) ≈ kron(A,b,C)
 
         @test_throws MethodError BlockKron()
         @test_throws MethodError BlockKron(a)