Skip to content
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

Jump to bottom

More efficient way of creating condensed sparsity pattern #436

Merged
merged 2 commits into from
Oct 5, 2022
Merged

More efficient way of creating condensed sparsity pattern #436

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ad853f1
More efficient way of creating condensed sparse
lijas Mar 8, 2022
e17d09c
Fix w.r.t review
lijas Oct 5, 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
73 changes: 51 additions & 22 deletions src/Dofs/ConstraintHandler.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ Collection of constraints.
"""
struct ConstraintHandler{DH<:AbstractDofHandler,T}
dbcs::Vector{Dirichlet}
acs::Vector{AffineConstraint}
acs::Vector{AffineConstraint{T}}
prescribed_dofs::Vector{Int}
free_dofs::Vector{Int}
inhomogeneities::Vector{T}
Expand All @@ -74,7 +74,7 @@ end

function ConstraintHandler(dh::AbstractDofHandler)
@assert isclosed(dh)
ConstraintHandler(Dirichlet[], AffineConstraint[], Int[], Int[], Float64[], Dict{Int,Int}(), BCValues{Float64}[], dh, ScalarWrapper(false))
ConstraintHandler(Dirichlet[], AffineConstraint{Float64}[], Int[], Int[], Float64[], Dict{Int,Int}(), BCValues{Float64}[], dh, ScalarWrapper(false))
end

"""
Expand Down Expand Up @@ -588,12 +588,18 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con
end

# Similar to Ferrite._condense!(K, ch), but only add the non-zero entries to K (that arises from the condensation process)
function _condense_sparsity_pattern!(K::SparseMatrixCSC, acs::Vector{AffineConstraint})
function _condense_sparsity_pattern!(K::SparseMatrixCSC{T}, acs::Vector{AffineConstraint{T}}) where T
ndofs = size(K, 1)

(length(acs) == 0) && return
# Store linear constraint index for each constrained dof
distribute = Dict{Int,Int}(acs[c].constrained_dof => c for c in 1:length(acs))

#Adding new entries to K is extremely slow, so create a new sparsity triplet for the condensed sparsity pattern
N = length(acs)*2 # TODO: Better size estimate for additional condensed sparsity pattern.
I = Int[]; resize!(I, N)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Could use sizehint! here instead and then use push! instead of setindex!.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yes, and perhaps in _create_sparsity_pattern() aswell? How much does push! Grow the vector if it grows beyond the size hint?

J = Int[]; resize!(J, N)

cnt = 0
for col in 1:ndofs
# Since we will possibly be pushing new entries to K, the field K.rowval will grow.
# Therefor we must extract this before iterating over K
Expand All @@ -606,7 +612,10 @@ function _condense_sparsity_pattern!(K::SparseMatrixCSC, acs::Vector{AffineConst
if drow != 0
ac = acs[drow]
for (d, _) in ac.entries
add_entry!(K, d, col)
if !_addindex_sparsematrix!(K, 0.0, d, col)
cnt += 1
_add_or_grow(cnt, I, J, d, col)
end
end
end
end
Expand All @@ -616,24 +625,41 @@ function _condense_sparsity_pattern!(K::SparseMatrixCSC, acs::Vector{AffineConst
if drow == 0
ac = acs[dcol]
for (d, _) in ac.entries
add_entry!(K, row, d)
if !_addindex_sparsematrix!(K, 0.0, row, d)
cnt += 1
_add_or_grow(cnt, I, J, row, d)
end
end
else
ac1 = acs[dcol]
for (d1, _) in ac1.entries
ac2 = acs[distribute[row]]
for (d2, _) in ac2.entries
add_entry!(K, d1, d2)
if !_addindex_sparsematrix!(K, 0.0, d1, d2)
cnt += 1
_add_or_grow(cnt, I, J, d1, d2)
end
end
end
end
end
end
end

resize!(I, cnt)
resize!(J, cnt)

# Fill the sparse matrix with a non-zero value so that :+ operation does not remove entries with value zero.
V = fill(1.0, length(I))
K2 = sparse(I, J, V, ndofs, ndofs)

K .+= K2

return nothing
end

# Condenses K and f: C'*K*C, C'*f, in-place assuming the sparsity pattern is correct
function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineConstraint})
function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineConstraint{T}}) where T

ndofs = size(K, 1)
condense_f = !(length(f) == 0)
Expand All @@ -652,7 +678,7 @@ function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineCon
ac = acs[drow]
for (d, v) in ac.entries
Kval = K.nzval[a]
_addindex_sparsematrix!(K, v * Kval, d, col)
_addindex_sparsematrix!(K, v * Kval, d, col) || _sparsity_error()
end

# Perform f - K*g. However, this has already been done in outside this functions so we skip this.
Expand All @@ -669,15 +695,15 @@ function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineCon
ac = acs[dcol]
for (d,v) in ac.entries
Kval = K.nzval[a]
_addindex_sparsematrix!(K, v * Kval, row, d)
_addindex_sparsematrix!(K, v * Kval, row, d) || _sparsity_error()
end
else
ac1 = acs[dcol]
for (d1,v1) in ac1.entries
ac2 = acs[drow]
for (d2,v2) in ac2.entries
Kval = K.nzval[a]
_addindex_sparsematrix!(K, v1 * v2 * Kval, d1, d2)
_addindex_sparsematrix!(K, v1 * v2 * Kval, d1, d2) || _sparsity_error()
end
end
end
Expand All @@ -695,8 +721,20 @@ function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineCon
end
end

_sparsity_error() = error("Sparsity pattern missing entries for the condensation pattern. Make sure to call `create_sparsity_pattern(dh::DofHandler, ch::ConstraintHandler) when using linear constraints.`")

function _add_or_grow(cnt::Int, I::Vector{Int}, J::Vector{Int}, dofi::Int, dofj::Int)
if cnt > length(J)
resize!(I, trunc(Int, length(I) * 1.5))
resize!(J, trunc(Int, length(J) * 1.5))
end
I[cnt] = dofi
J[cnt] = dofj
end

# Copied from SparseArrays._setindex_scalar!(...)
# Custom SparseArrays._setindex_scalar!() that throws error if entry K(_i,_j) does not exist
# Returns true if it successfully added the new value, returns false otherwise.
function _addindex_sparsematrix!(A::SparseMatrixCSC{Tv,Ti}, v::Tv, i::Ti, j::Ti) where {Tv, Ti}
if !((1 <= i <= size(A, 1)) & (1 <= j <= size(A, 2)))
throw(BoundsError(A, (i,j)))
Expand All @@ -707,18 +745,9 @@ function _addindex_sparsematrix!(A::SparseMatrixCSC{Tv,Ti}, v::Tv, i::Ti, j::Ti)
if searchk <= coljlastk && rowvals(A)[searchk] == i
# Column j contains entry A[i,j]. Update and return
nonzeros(A)[searchk] += v
return A
end
error("Sparsity pattern missing entries for the condensation pattern. Make sure to call `create_sparsity_pattern(dh::DofHandler, ch::ConstraintHandler) when using linear constraints.`")
end

# A[i,j] += 0.0 does not add entries to sparse matrices, so we need to first add 1.0, and then remove it
# TODO: Maybe this can be done for vectors i and j instead of doing it individually?
function add_entry!(A::SparseMatrixCSC, i::Int, j::Int)
if iszero(A[i,j]) # Check first if zero to not remove already non-zero entries
A[i,j] = 1
A[i,j] = 0
return true
end
return false
end

"""
Expand Down
15 changes: 9 additions & 6 deletions src/Dofs/DofHandler.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -356,6 +356,7 @@ See the [Sparsity Pattern](@ref) section of the manual.
create_symmetric_sparsity_pattern(dh::DofHandler) = Symmetric(_create_sparsity_pattern(dh, nothing, true), :U)

function _create_sparsity_pattern(dh::DofHandler, ch#=::Union{ConstraintHandler, Nothing}=#, sym::Bool)
@assert isclosed(dh)
ncells = getncells(dh.grid)
n = ndofs_per_cell(dh)
N = sym ? div(n*(n+1), 2) * ncells : n^2 * ncells
Expand Down Expand Up @@ -392,17 +393,19 @@ function _create_sparsity_pattern(dh::DofHandler, ch#=::Union{ConstraintHandler,

resize!(I, cnt)
resize!(J, cnt)
V = zeros(length(I))
K = sparse(I, J, V)

# Add entries to K corresponding to condensation due the linear constraints
# Note, this requires the K matrix, which is why we can't push!() to the I,J,V
# triplet directly.
# If ConstraintHandler is given, create the condensation pattern due to affine constraints
if ch !== nothing
@assert isclosed(ch)

V = ones(length(I))
K = sparse(I, J, V, ndofs(dh), ndofs(dh))
_condense_sparsity_pattern!(K, ch.acs)
fill!(K.nzval, 0.0)
else
V = zeros(length(I))
K = sparse(I, J, V, ndofs(dh), ndofs(dh))
end

return K
end

Expand Down
17 changes: 10 additions & 7 deletions src/Dofs/MixedDofHandler.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -371,6 +371,8 @@ end

# TODO if not too slow it can replace the "Grid-version"
function _create_sparsity_pattern(dh::MixedDofHandler, ch#=::Union{ConstraintHandler,Nothing}=#, sym::Bool)
@assert isclosed(dh)

ncells = getncells(dh.grid)
N::Int = 0
for element_id = 1:ncells # TODO check for correctness
Expand Down Expand Up @@ -410,19 +412,20 @@ function _create_sparsity_pattern(dh::MixedDofHandler, ch#=::Union{ConstraintHan
end
resize!(I, cnt)
resize!(J, cnt)
V = zeros(length(I))
K = sparse(I, J, V)

# Add entries to K corresponding to condensation due the linear constraints
# Note, this requires the K matrix, which is why we can't push!() to the I,J,V
# triplet directly.
# If ConstraintHandler is given, create the condensation pattern due to affine constraints
if ch !== nothing
@assert isclosed(ch)

V = ones(length(I))
K = sparse(I, J, V, ndofs(dh), ndofs(dh))
_condense_sparsity_pattern!(K, ch.acs)
fill!(K.nzval, 0.0)
else
V = zeros(length(I))
K = sparse(I, J, V, ndofs(dh), ndofs(dh))
end

return K

end

create_sparsity_pattern(dh::MixedDofHandler) = _create_sparsity_pattern(dh, nothing, false)
Expand Down