|
| 1 | +module FerriteSparseMatrixCSR |
| 2 | + |
| 3 | +using Ferrite, SparseArrays, SparseMatricesCSR |
| 4 | +import Ferrite: AbstractSparsityPattern, CSRAssembler |
| 5 | +import Base: @propagate_inbounds |
| 6 | + |
| 7 | +# Could be generalized if https://github.com/JuliaSparse/SparseArrays.jl/pull/546 is merged |
| 8 | +function Ferrite.start_assemble(K::SparseMatrixCSR{<:Any, T}, f::Vector = T[]; fillzero::Bool = true, maxcelldofs_hint::Int = 0) where {T} |
| 9 | + fillzero && (Ferrite.fillzero!(K); Ferrite.fillzero!(f)) |
| 10 | + return CSRAssembler(K, f, zeros(Int, maxcelldofs_hint), zeros(Int, maxcelldofs_hint)) |
| 11 | +end |
| 12 | + |
| 13 | +@propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool) |
| 14 | + current_row = 1 |
| 15 | + ld = length(dofs) |
| 16 | + return @inbounds for Krow in sorteddofs |
| 17 | + maxlookups = sym ? current_row : ld |
| 18 | + Kerow = permutation[current_row] |
| 19 | + ci = 1 # col index pointer for the local matrix |
| 20 | + Ci = 1 # col index pointer for the global matrix |
| 21 | + nzr = nzrange(K, Krow) |
| 22 | + while Ci <= length(nzr) && ci <= maxlookups |
| 23 | + C = nzr[Ci] |
| 24 | + Kcol = K.colval[C] |
| 25 | + Kecol = permutation[ci] |
| 26 | + val = Ke[Kerow, Kecol] |
| 27 | + if Kcol == dofs[Kecol] |
| 28 | + # Match: add the value (if non-zero) and advance the pointers |
| 29 | + if !iszero(val) |
| 30 | + K.nzval[C] += val |
| 31 | + end |
| 32 | + ci += 1 |
| 33 | + Ci += 1 |
| 34 | + elseif Kcol < dofs[Kecol] |
| 35 | + # No match yet: advance the global matrix row pointer |
| 36 | + Ci += 1 |
| 37 | + else # Kcol > dofs[Kecol] |
| 38 | + # No match: no entry exist in the global matrix for this row. This is |
| 39 | + # allowed as long as the value which would have been inserted is zero. |
| 40 | + iszero(val) || Ferrite._missing_sparsity_pattern_error(Krow, Kcol) |
| 41 | + # Advance the local matrix row pointer |
| 42 | + ci += 1 |
| 43 | + end |
| 44 | + end |
| 45 | + # Make sure that remaining entries in this column of the local matrix are all zero |
| 46 | + for i in ci:maxlookups |
| 47 | + if !iszero(Ke[Kerow, permutation[i]]) |
| 48 | + Ferrite._missing_sparsity_pattern_error(Krow, sorteddofs[i]) |
| 49 | + end |
| 50 | + end |
| 51 | + current_row += 1 |
| 52 | + end |
| 53 | +end |
| 54 | + |
| 55 | +function Ferrite.zero_out_rows!(K::SparseMatrixCSR, ch::ConstraintHandler) |
| 56 | + @debug @assert issorted(ch.prescribed_dofs) |
| 57 | + for row in ch.prescribed_dofs |
| 58 | + r = nzrange(K, row) |
| 59 | + K.nzval[r] .= 0.0 |
| 60 | + end |
| 61 | + return |
| 62 | +end |
| 63 | + |
| 64 | +function Ferrite.zero_out_columns!(K::SparseMatrixCSR, ch::ConstraintHandler) |
| 65 | + colval = K.colval |
| 66 | + nzval = K.nzval |
| 67 | + return @inbounds for i in eachindex(colval, nzval) |
| 68 | + if haskey(ch.dofmapping, colval[i]) |
| 69 | + nzval[i] = 0 |
| 70 | + end |
| 71 | + end |
| 72 | +end |
| 73 | + |
| 74 | +function Ferrite.allocate_matrix(::Type{SparseMatrixCSR}, sp::AbstractSparsityPattern) |
| 75 | + return Ferrite.allocate_matrix(SparseMatrixCSR{1, Float64, Int64}, sp) |
| 76 | +end |
| 77 | + |
| 78 | +function Ferrite.allocate_matrix(::Type{SparseMatrixCSR{1, Tv, Ti}}, sp::AbstractSparsityPattern) where {Tv, Ti} |
| 79 | + return _allocate_matrix(SparseMatrixCSR{1, Tv, Ti}, sp, false) |
| 80 | +end |
| 81 | + |
| 82 | +function _allocate_matrix(::Type{SparseMatrixCSR{1, Tv, Ti}}, sp::AbstractSparsityPattern, sym::Bool) where {Tv, Ti} |
| 83 | + # 1. Setup rowptr |
| 84 | + rowptr = zeros(Ti, Ferrite.getnrows(sp) + 1) |
| 85 | + rowptr[1] = 1 |
| 86 | + for (row, colidxs) in enumerate(Ferrite.eachrow(sp)) |
| 87 | + for col in colidxs |
| 88 | + sym && row > col && continue |
| 89 | + rowptr[row + 1] += 1 |
| 90 | + end |
| 91 | + end |
| 92 | + cumsum!(rowptr, rowptr) |
| 93 | + nnz = rowptr[end] - 1 |
| 94 | + # 2. Allocate colval and nzval now that nnz is known |
| 95 | + colval = Vector{Ti}(undef, nnz) |
| 96 | + nzval = zeros(Tv, nnz) |
| 97 | + # 3. Populate colval. |
| 98 | + k = 1 |
| 99 | + for (row, colidxs) in zip(1:Ferrite.getnrows(sp), Ferrite.eachrow(sp)) # pairs(eachrow(sp)) |
| 100 | + for col in colidxs |
| 101 | + sym && row > col && continue |
| 102 | + colval[k] = col |
| 103 | + k += 1 |
| 104 | + end |
| 105 | + end |
| 106 | + S = SparseMatrixCSR{1}(Ferrite.getnrows(sp), Ferrite.getncols(sp), rowptr, colval, nzval) |
| 107 | + return S |
| 108 | +end |
| 109 | + |
| 110 | +end |
0 commit comments