Skip to content

Create blockcholesky.jl #126

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 21 commits into from
Jun 16, 2022
Merged

Create blockcholesky.jl #126

Changes from 1 commit
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
9c47f78
Create blockcholesky.jl
felixw98 Jul 14, 2020
030a054
Add files via upload
felixw98 Jul 18, 2020
b977d8d
overload cholesk!
dlfivefifty Jul 19, 2020
0ff3e67
Including the test file in REPL will do the test
felixw98 Jul 20, 2020
dc54673
Including test_cholesky.jl in runtests.jl
felixw98 Jul 20, 2020
5605d8f
Added lower Cholesky factorization
felixw98 Jul 21, 2020
bb7355b
Added PD check
felixw98 Jul 26, 2020
2d782fe
Changed the function name and types
felixw98 Jul 27, 2020
72259da
Added blockbanded case
felixw98 Aug 23, 2020
883da65
Rearrange
felixw98 Aug 30, 2020
09a8a54
Capable on both BlockArray and BlockBandedArray
felixw98 Sep 1, 2020
b6f7d3a
Changes
felixw98 Sep 6, 2020
dc6e020
Merge branch 'master' into pr/126
dlfivefifty Jun 14, 2022
8793acd
Update Choleksy
dlfivefifty Jun 14, 2022
f994fc0
tests pass
dlfivefifty Jun 14, 2022
44b1490
Support sparse block structure
dlfivefifty Jun 14, 2022
2d43bda
Update blockcholesky.jl
dlfivefifty Jun 14, 2022
025a4ce
Fix tests on 1.8, use ArrayLayouts.ldiv!
dlfivefifty Jun 15, 2022
04af833
Speed up BlockKron construction
dlfivefifty Jun 16, 2022
f1059ef
Update test_blockproduct.jl
dlfivefifty Jun 16, 2022
99f54d5
increase cov
dlfivefifty Jun 16, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 27 additions & 0 deletions src/blockcholesky.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@

# Function for cholesky decomposition on block matries
# 'cholesky' is the build-in one in LAPACK

function Bcholesky(A::Symmetric{<:Any,<:BlockArray})
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In src/BlockArrays.jl, add

import LinearAlgebra: cholesky, cholesky!

And then you can overload cholesky and cholesky!

chol_U = BlockArray{Float64}(zeros(size(A)), fill(size(A)[1]÷blocksize(A)[1], blocksize(A)[1]), fill(size(A)[1]÷blocksize(A)[1], blocksize(A)[1]))

# Initializing the first role of blocks
chol_U[Block(1,1)] = cholesky!(A[Block(1,1)]).U
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When we change this to be properly in-place, we will want to replace this with cholesky!(Symmetric(getblock(parent(A),1,1))) which will modify A.

for j = 2:blocksize(A)[1]
chol_U[Block(1,j)] = chol_U[Block(1,1)]' \ A[Block(1,j)]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When this is in-place it will look something like:

P = parent(A)
ldiv!(UpperTriangular(getblock(P,1,1))', getblock(P,1,j))

end

# For the left blocks
for i = 2:blocksize(A)[1]
for j = i:blocksize(A)[1]
if j == i
chol_U[Block(i,j)] = cholesky!(A[Block(i,j)] - sum(chol_U[Block(k,j)]'chol_U[Block(k,j)] for k in 1:j-1)).U
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here we will want to use mul! to do this in-place. Something like:

Pij = getblock(P,i,j)
for k = 1:j-1
    mul!(Pij,getblock(P,k,j)',getblock(P,k,j),-1.0,1.0)
end
cholesky!(Pij)

else
chol_U[Block(i,j)] = chol_U[Block(i,i)]' \ (A[Block(i,j)] - sum(chol_U[Block(k,i)]'chol_U[Block(k,j)] for k in 1:j-1))
end
end
end

return chol_U

end