Skip to content

Commit

Permalink
Resolved merge conflicts
Browse files Browse the repository at this point in the history
  • Loading branch information
@bonevbs
bonevbs committed Apr 9, 2021
2 parents 2564a52 + a3e6f8e commit 3b3a11b
Show file tree
Hide file tree
Showing 9 changed files with 663 additions and 615 deletions.
4 changes: 2 additions & 2 deletions src/HssMatrices.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ module HssMatrices
# linearmap.jl
export LinearMap, HermitianLinearMap
# hssmatrix.jl
export HssLeaf, HssNode, HssMatrix, hss, isleaf, isbranch, ishss, hssrank, full, checkdims, prune_leaves!, blkdiagm
export HssMatrix, hss, isleaf, isbranch, ishss, hssrank, full, checkdims, prune_leaves!, blkdiagm
# prrqr.jl
export prrqr, prrqr!
# compression.jl
Expand Down Expand Up @@ -97,7 +97,7 @@ module HssMatrices
opts.noversampling ≥ 1 || throw(ArgumentError("noversampling"))
end

function setopts(opts::HssOptions=HssOptions(Float64); args...)
function setopts(opts=HssOptions(Float64); args...)
opts = copy(opts; args...)
chkopts!(opts)
global atol = opts.atol
Expand Down
224 changes: 122 additions & 102 deletions src/compression.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,24 +41,28 @@ function compress(A::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, opts::HssOpt
Brow = Array{T}(undef, m, 0)
Bcol = Array{T}(undef, 0, n)
if isleaf(rcl) && isleaf(ccl)
hssA, _, _ = _compress!(A, Brow, Bcol, rcl.data, ccl.data, opts.atol, opts.rtol)
hssA, _, _ = _compress!(A, Brow, Bcol, rcl.data, ccl.data, opts.atol, opts.rtol; rootnode=true)
elseif isbranch(rcl) && isbranch(ccl)
hssA, _, _ = _compress!(A, Brow, Bcol, rcl, ccl, opts.atol, opts.rtol)
hssA, _, _ = _compress!(A, Brow, Bcol, rcl, ccl, opts.atol, opts.rtol; rootnode=true)
else
throw(ArgumentError("row and column clusters are not compatible"))
end
return hssA
end

# leaf node function for compression
function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rows::UnitRange{Int}, cols::UnitRange{Int}, atol::Float64, rtol::Float64) where T
U, Brow = _compress_block!(Brow, atol, rtol)
V, Bcol = _compress_block!(Bcol', atol, rtol) #TODO: write code that is better at dealing with Julia's lazy transpose
return HssLeaf(A[rows, cols], U, V), Brow, copy(Bcol')
function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rows::UnitRange{Int}, cols::UnitRange{Int}, atol::Float64, rtol::Float64; rootnode=false) where T
if rootnode
return HssMatrix(A[rows, cols]), Brow, Bcol
else
U, Brow = _compress_block!(Brow, atol, rtol)
V, Bcol = _compress_block!(Bcol', atol, rtol) #TODO: write code that is better at dealing with Julia's lazy transpose
return HssMatrix(A[rows, cols], U, V), Brow, copy(Bcol')
end
end

# branch node function for compression
function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, atol::Float64, rtol::Float64) where T
function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, atol::Float64, rtol::Float64; rootnode=false) where T
m1 = length(rcl.left.data); m2 = length(rcl.right.data)
n1 = length(ccl.left.data); n2 = length(ccl.right.data)

Expand Down Expand Up @@ -96,20 +100,23 @@ function _compress!(A::Matrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, rcl::Cluster
# clean up stuff from the front and form the composed HSS block row/col for compression
Brow = [Brow1[:, n2+1:end]; Brow2[:, rn1+1:end]]
Bcol = [Bcol1[m2+1:end, :] Bcol2[rm1+1:end, :]]

# do the actual compression and disentangle blocks of the translation operators
R, Brow = _compress_block!(Brow, atol, rtol)
R1 = R[1:rm1, :]
R2 = R[rm1+1:end, :]

X = copy(Bcol')
W, Bcol = _compress_block!(copy(Bcol'), atol, rtol);
Bcol = copy(Bcol')
W1 = W[1:rn1, :]
W2 = W[rn1+1:end, :]

# call recursively the
hssA = HssNode(A11, A22, B12, B21, R1, W1, R2, W2)

if !rootnode
# do the actual compression and disentangle blocks of the translation operators
R, Brow = _compress_block!(Brow, atol, rtol)
R1 = R[1:rm1, :]
R2 = R[rm1+1:end, :]

X = copy(Bcol')
W, Bcol = _compress_block!(copy(Bcol'), atol, rtol);
Bcol = copy(Bcol')
W1 = W[1:rn1, :]
W2 = W[rn1+1:end, :]

hssA = HssMatrix(A11, A22, B12, B21, R1, W1, R2, W2, rootnode)
else
hssA = HssMatrix(A11, A22, B12, B21, rootnode)
end
return hssA, Brow, Bcol
end

Expand Down Expand Up @@ -183,7 +190,7 @@ function recompress!(hssA::HssMatrix{T}, opts::HssOptions=HssOptions(T); args...
return hssA
end

function _recompress!(hssA::HssNode{T}, Brow::Matrix{T}, Bcol::Matrix{T}, atol, rtol) where T
function _recompress!(hssA::HssMatrix{T}, Brow::Matrix{T}, Bcol::Matrix{T}, atol, rtol) where T
# compress B12
Brow1 = [hssA.B12 hssA.R1*Brow]
Bcol2 = [hssA.B12' hssA.W2*Bcol]
Expand Down Expand Up @@ -270,7 +277,7 @@ function randcompress(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTr
Ωrow = randn(m, k+r)
Scol = A*Ωcol # this should invoke the magic of the linearoperator.jl type
Srow = A'*Ωrow
hssA = _randcompress!(hssA, A, Scol, Srow, Ωcol, Ωrow, 0, 0, opts.atol, opts.rtol; rootnode=true)
hssA, _, _, _, _, _, _, _, _ = _randcompress!(hssA, A, Scol, Srow, Ωcol, Ωrow, 0, 0, opts.atol, opts.rtol; rootnode=true)
return hssA
end

Expand All @@ -290,7 +297,7 @@ Randomized HSS compression.
julia> hssA = randcompress_adaptive(A, rcl, ccl, kest=20)
```
"""
function randcompress_adaptive(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree, opts::HssOptions=HssOptions(T); kest=10, args...) where T
function randcompress_adaptive(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree, opts=HssOptions(T); kest=10, args...) where T
opts = copy(opts; args...)
chkopts!(opts)
m, n = size(A)
Expand All @@ -310,7 +317,7 @@ function randcompress_adaptive(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::

while failed && k < maxrank
# TODO: In further versions we might want to re-use the information gained during previous attempts
hssA = _randcompress!(hssA, A, Scol, Srow, Ωcol, Ωrow, 0, 0, opts.atol, opts.rtol; rootnode=true)
hssA, _, _, _, _, _, _, _, _ = _randcompress!(hssA, A, Scol, Srow, Ωcol, Ωrow, 0, 0, opts.atol, opts.rtol; rootnode=true)

Ωcol_test = randn(n, bs)
Ωrow_test = randn(m, bs)
Expand All @@ -336,105 +343,118 @@ function randcompress_adaptive(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::
end

# this function compresses given the sampling matrix of rank k
function _randcompress!(hssA::HssLeaf, A, Scol::Matrix, Srow::Matrix, Ωcol::Matrix, Ωrow::Matrix, ro::Int, co::Int, atol::Float64, rtol::Float64; rootnode=false)
Scol .= Scol .- hssA.D * Ωcol
Srow .= Srow .- hssA.D' * Ωrow

if rootnode return hssA end

# take care of column-space
Xcol, Jcol = _interpolate(Scol', atol, rtol)
hssA.U = Xcol'
Scol = Scol[Jcol, :]
U = Xcol'
Jcol .= ro .+ Jcol

# same for the row-space
Xrow, Jrow = _interpolate(Srow', atol, rtol)
hssA.V = Xrow'
Srow = Srow[Jrow, :]
V = Xrow'
Jrow .= co .+ Jrow

return hssA, Scol, Srow, Ωcol, Ωrow, Jcol, Jrow, U, V
end
function _randcompress!(hssA::HssNode, A, Scol::Matrix, Srow::Matrix, Ωcol::Matrix, Ωrow::Matrix, ro::Int, co::Int, atol::Float64, rtol::Float64; rootnode=false)
m1, n1 = hssA.sz1; m2, n2 = hssA.sz2
hssA.A11, Scol1, Srow1, Ωcol1, Ωrow1, Jcol1, Jrow1, U1, V1 = _randcompress!(hssA.A11, A, Scol[1:m1, :], Srow[1:n1, :], Ωcol[1:n1, :], Ωrow[1:m1, :], ro, co, atol, rtol)
hssA.A22, Scol2, Srow2, Ωcol2, Ωrow2, Jcol2, Jrow2, U2, V2 = _randcompress!(hssA.A22, A, Scol[m1+1:end, :], Srow[n1+1:end, :], Ωcol[n1+1:end, :], Ωrow[m1+1:end, :], ro+m1, co+n1, atol, rtol)

# update the sampling matrix based on the extracted generators
Ωcol2 = V2' * Ωcol2
Ωcol1 = V1' * Ωcol1
Ωrow2 = U2' * Ωrow2
Ωrow1 = U1' * Ωrow1
# step 1 in the algorithm: look only at extracted rows/cols
Jcol = [Jcol1; Jcol2]
Jrow = [Jrow1; Jrow2]
Ωcol = [Ωcol1; Ωcol2]
Ωrow = [Ωrow1; Ωrow2]
# extract the correct generator blocks
hssA.B12 = A[Jcol1, Jrow2]
hssA.B21 = A[Jcol2, Jrow1]
# subtract the diagonal block
Scol = [Scol1 .- hssA.B12 * Ωcol2; Scol2 .- hssA.B21 * Ωcol1 ];
Srow = [Srow1 .- hssA.B21' * Ωrow2; Srow2 .- hssA.B12' * Ωrow1 ];

if rootnode
rk1, wk1 = gensize(hssA.A11)
rk2, wk2 = gensize(hssA.A22)
hssA.R1 = Matrix{eltype(hssA)}(undef, rk1, 0)
hssA.R2 = Matrix{eltype(hssA)}(undef, rk2, 0)
hssA.W1 = Matrix{eltype(hssA)}(undef, wk1, 0)
hssA.W2 = Matrix{eltype(hssA)}(undef, wk2, 0)

return hssA
function _randcompress!(hssA::HssMatrix, A::AbstractMatOrLinOp, Scol::Matrix, Srow::Matrix, Ωcol::Matrix, Ωrow::Matrix, ro::Int, co::Int, atol::Float64, rtol::Float64; rootnode=false)
if isleaf(hssA)
Scol .= Scol .- hssA.D * Ωcol
Srow .= Srow .- hssA.D' * Ωrow

#if rootnode return hssA, Scol, Srow, Ωcol, Ωrow, Jcol, Jrow, hssA.U, hssA.V end

# take care of column-space
Xcol, Jcol = _interpolate(Scol', atol, rtol)
hssA.U = Xcol'
Scol = Scol[Jcol, :]
U = Xcol'
Jcol .= ro .+ Jcol

# same for the row-space
Xrow, Jrow = _interpolate(Srow', atol, rtol)
hssA.V = Xrow'
Srow = Srow[Jrow, :]
V = Xrow'
Jrow .= co .+ Jrow
else
# take care of the columns
Xcol, Jcolloc = _interpolate(Scol', atol, rtol)
hssA.R1 = Xcol[:, 1:size(Scol1, 1)]'
hssA.R2 = Xcol[:, size(Scol1, 1)+1:end]'
Scol = Scol[Jcolloc, :]
Jcol = Jcol[Jcolloc]
U = [hssA.R1; hssA.R2]
m1, n1 = hssA.sz1; m2, n2 = hssA.sz2
hssA.A11, Scol1, Srow1, Ωcol1, Ωrow1, Jcol1, Jrow1, U1, V1 = _randcompress!(hssA.A11, A, Scol[1:m1, :], Srow[1:n1, :], Ωcol[1:n1, :], Ωrow[1:m1, :], ro, co, atol, rtol)
hssA.A22, Scol2, Srow2, Ωcol2, Ωrow2, Jcol2, Jrow2, U2, V2 = _randcompress!(hssA.A22, A, Scol[m1+1:end, :], Srow[n1+1:end, :], Ωcol[n1+1:end, :], Ωrow[m1+1:end, :], ro+m1, co+n1, atol, rtol)

# take care of the rows
Xrow, Jrowloc = _interpolate(Srow', atol, rtol)
hssA.W1 = Xrow[:, 1:size(Srow1, 1)]'
hssA.W2 = Xrow[:, size(Srow1, 1)+1:end]'
Srow = Srow[Jrowloc, :]
Jrow = Jrow[Jrowloc]
V = [hssA.W1; hssA.W2]

return hssA, Scol, Srow, Ωcol, Ωrow, Jcol, Jrow, U, V
# update the sampling matrix based on the extracted generators
Ωcol2 = V2' * Ωcol2
Ωcol1 = V1' * Ωcol1
Ωrow2 = U2' * Ωrow2
Ωrow1 = U1' * Ωrow1
# step 1 in the algorithm: look only at extracted rows/cols
Jcol = [Jcol1; Jcol2]
Jrow = [Jrow1; Jrow2]
Ωcol = [Ωcol1; Ωcol2]
Ωrow = [Ωrow1; Ωrow2]
# extract the correct generator blocks
hssA.B12 = A[Jcol1, Jrow2]
hssA.B21 = A[Jcol2, Jrow1]
# subtract the diagonal block
Scol = [Scol1 .- hssA.B12 * Ωcol2; Scol2 .- hssA.B21 * Ωcol1 ];
Srow = [Srow1 .- hssA.B21' * Ωrow2; Srow2 .- hssA.B12' * Ωrow1 ];

if rootnode
rk1, wk1 = gensize(hssA.A11)
rk2, wk2 = gensize(hssA.A22)
hssA.R1 = Matrix{eltype(hssA)}(undef, rk1, 0)
hssA.R2 = Matrix{eltype(hssA)}(undef, rk2, 0)
hssA.W1 = Matrix{eltype(hssA)}(undef, wk1, 0)
hssA.W2 = Matrix{eltype(hssA)}(undef, wk2, 0)
hssA.rootnode = true

U = similar(U1, rk1+rk2, 0)
V = similar(V1, wk1+wk2, 0)
else
# take care of the columns
Xcol, Jcolloc = _interpolate(Scol', atol, rtol)
hssA.R1 = Xcol[:, 1:size(Scol1, 1)]'
hssA.R2 = Xcol[:, size(Scol1, 1)+1:end]'
Scol = Scol[Jcolloc, :]
Jcol = Jcol[Jcolloc]
U = [hssA.R1; hssA.R2]

# take care of the rows
Xrow, Jrowloc = _interpolate(Srow', atol, rtol)
hssA.W1 = Xrow[:, 1:size(Srow1, 1)]'
hssA.W2 = Xrow[:, size(Srow1, 1)+1:end]'
Srow = Srow[Jrowloc, :]
Jrow = Jrow[Jrowloc]
V = [hssA.W1; hssA.W2]
end
end
return hssA, Scol, Srow, Ωcol, Ωrow, Jcol, Jrow, U, V
end

# interpolative decomposition
# still gotta figure out which qr decomposition to use
function _interpolate(A::AbstractMatrix{T}, atol::Float64, rtol::Float64) where T
size(A,2) == 0 && return Matrix{T}(undef, 0,0), []
size(A,2) == 0 && return Matrix{eltype(A)}(undef, 0, 0), Vector{Int}()
#_, R, p = prrqr(A, tol; reltol=reltol)
#rk = min(size(R)...)
_, R, p = qr(A, Val(true))
tol = min( atol, rtol * abs(R[1,1]) )
rk = sum(abs.(diag(R)) .> tol)
J = p[1:rk];
J = p[1:rk]
X = R[1:rk, 1:rk]\R[1:rk,invperm(p)]
return X, J
end

# extracts the block-diagonal of a matrix as HSS matrix of rank 0
function hss_blkdiag(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree) where T
function hss_blkdiag(A::AbstractMatOrLinOp{T}, rcl::ClusterTree, ccl::ClusterTree; rootnode=true) where T
m = length(rcl.data); n = length(ccl.data)
if isleaf(rcl) # only check row cluster as we have already checked cluster equality
D = A[rcl.data, ccl.data]
return HssLeaf(convert(Matrix{T}, D))
D = convert(Matrix{T}, A[rcl.data, ccl.data])
if rootnode
return HssMatrix(D)
else
return HssMatrix(D, Matrix{T}(undef,m,0), Matrix{T}(undef,n,0))
end
elseif isbranch(rcl)
A11 = hss_blkdiag(A, rcl.left, ccl.left)
A22 = hss_blkdiag(A, rcl.right, ccl.right)
A11 = hss_blkdiag(A, rcl.left, ccl.left; rootnode=false)
A22 = hss_blkdiag(A, rcl.right, ccl.right; rootnode=false)
B12 = Matrix{T}(undef,0,0)
B21 = Matrix{T}(undef,0,0)
return HssNode(A11, A22, B12, B21)
if rootnode
return HssMatrix(A11, A22, B12, B21)
else
R1 = Matrix{T}(undef,0,0)
W1 = Matrix{T}(undef,0,0)
R2 = Matrix{T}(undef,0,0)
W2 = Matrix{T}(undef,0,0)
return HssMatrix(A11, A22, B12, B21, R1, W1, R2, W2)
end
else
throw(ErrorException("Encountered node with only one child. Make sure that each node in the binary tree has either two children or is a leaf."))
end
Expand Down
8 changes: 4 additions & 4 deletions src/constructors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,16 +4,16 @@

function lowrank2hss(U::Matrix{T}, V::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree) where T
(k = size(U,2)) == size(V,2) || throw(ArgumentError("second dimension of U and V must agree"))
return _lowrank2hss(U, V, rcl, ccl, k)
return _lowrank2hss(U, V, rcl, ccl, k; rootnode=true)
end

function _lowrank2hss(U::Matrix{T}, V::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, k::Int) where T
function _lowrank2hss(U::Matrix{T}, V::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, k::Int; rootnode=false) where T
if isleaf(rcl) && isleaf(ccl)
HssLeaf(U[rcl.data,:]*V[ccl.data,:]', U[rcl.data,:], V[ccl.data,:])
HssMatrix(U[rcl.data,:]*V[ccl.data,:]', U[rcl.data,:], V[ccl.data,:], rootnode)
elseif isbranch(rcl) && isbranch(ccl)
A11 = _lowrank2hss(U, V, rcl.left, ccl.left, k)
A22 = _lowrank2hss(U, V, rcl.right, ccl.right, k)
HssNode(A11, A22, Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k))
HssMatrix(A11, A22, Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), rootnode)
else
throw(ArgumentError("row and column clusters are not compatible"))
end
Expand Down
Loading

2 comments on commit 3b3a11b

@bonevbs
Copy link
Owner Author

@bonevbs bonevbs commented on 3b3a11b Apr 9, 2021

Choose a reason for hiding this comment

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

@JuliaRegistrator register

Release notes:

  • re-wrote the underlying datastructure. HssLeaf and HssNode are replaced by HssMatrix.
  • achieved type-stability in random compression

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

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

Registration pull request created: JuliaRegistries/General/33911

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.2 -m "<description of version>" 3b3a11b92bc1921917a44cccef485af1f208bc88
git push origin v0.1.2

Please sign in to comment.