diff --git a/Project.toml b/Project.toml
index af5e162e41..7cfdb6e609 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,15 +18,18 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [weakdeps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
+SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 
 [extensions]
 FerriteBlockArrays = "BlockArrays"
 FerriteMetis = "Metis"
+FerriteSparseMatrixCSR = "SparseMatricesCSR"
 
 [compat]
 BlockArrays = "0.16, 1"
 EnumX = "1"
 Metis = "1.3"
+SparseMatricesCSR = "0.6"
 NearestNeighbors = "0.4"
 OrderedCollections = "1"
 Preferences = "1"
@@ -49,8 +52,9 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [targets]
-test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "Pkg", "NBInclude", "ProgressMeter", "Random", "SHA", "Test", "TimerOutputs", "Logging"]
+test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "Pkg", "NBInclude", "ProgressMeter", "Random", "SHA", "SparseMatricesCSR", "Test", "TimerOutputs", "Logging"]
diff --git a/docs/src/devdocs/assembler.md b/docs/src/devdocs/assembler.md
new file mode 100644
index 0000000000..28b9df8ade
--- /dev/null
+++ b/docs/src/devdocs/assembler.md
@@ -0,0 +1,37 @@
+# [Assembler](@id devdocs-assembler)
+
+The assembler handles the insertion of the element matrices and element vectors into the system matrix and right hand side. While the CSC and CSR formats are the most common sparse matrix formats in practice, users might want to have optimized custom matrix formats for their specific use-case. The default assembler [`AssemblerSparsityPattern`](@ref) should be able to handle most cases in practice. To support a custom format users have to dispatch the following functions:
+
+```@docs
+Ferrite._assemble_inner!
+Ferrite.zero_out_rows!
+Ferrite.zero_out_columns!
+```
+
+and the `AbstractSparseMatrix` interface for their custom matrix type. Optional dispatches to speed up operations might be
+
+```@docs
+Ferrite.add_inhomogeneities!
+```
+
+## Custom Assembler
+
+In case the default assembler is insufficient, users can implement a custom assemblers. For this, they can create a custom type and dispatch the following functions.
+
+```@docs
+Ferrite.matrix_handle
+Ferrite.vector_handle
+start_assemble!
+finish_assemble!
+assemble!
+```
+
+For local elimination support the following functions might also need custom dispatches
+
+```@docs
+Ferrite._condense_local!
+```
+
+## Custom Matrix Type Sparsity Pattern
+
+TODO `create_sparsity_pattern`
diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md
index 9c16b83d7c..f292cdb55d 100644
--- a/docs/src/devdocs/index.md
+++ b/docs/src/devdocs/index.md
@@ -5,5 +5,5 @@ developing the library.
 
 ```@contents
 Depth = 1
-Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "performance.md"]
+Pages = ["reference_cells.md", "interpolations.md", "elements.md", "FEValues.md", "dofhandler.md", "assembler.md", "performance.md"]
 ```
diff --git a/ext/FerriteSparseMatrixCSR.jl b/ext/FerriteSparseMatrixCSR.jl
new file mode 100644
index 0000000000..07d24605ae
--- /dev/null
+++ b/ext/FerriteSparseMatrixCSR.jl
@@ -0,0 +1,79 @@
+module FerriteSparseMatrixCSR
+
+using Ferrite, SparseArrays, SparseMatricesCSR
+import Ferrite: AbstractSparsityPattern
+import Base: @propagate_inbounds
+
+@propagate_inbounds function Ferrite._assemble_inner!(K::SparseMatrixCSR, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool)
+    current_row = 1
+    ld = length(dofs)
+    @inbounds for Krow in sorteddofs
+        maxlookups = sym ? current_row : ld
+        Kerow = permutation[current_row]
+        ci = 1 # col index pointer for the local matrix
+        Ci = 1 # col index pointer for the global matrix
+        nzr = nzrange(K, Krow)
+        while Ci <= length(nzr) && ci <= maxlookups
+            C = nzr[Ci]
+            Kcol = K.colval[C]
+            Kecol = permutation[ci]
+            val = Ke[Kerow, Kecol]
+            if Kcol == dofs[Kecol]
+                # Match: add the value (if non-zero) and advance the pointers
+                if !iszero(val)
+                    K.nzval[C] += val
+                end
+                ci += 1
+                Ci += 1
+            elseif Kcol < dofs[Kecol]
+                # No match yet: advance the global matrix row pointer
+                Ci += 1
+            else # Kcol > dofs[Kecol]
+                # No match: no entry exist in the global matrix for this row. This is
+                # allowed as long as the value which would have been inserted is zero.
+                iszero(val) || Ferrite._missing_sparsity_pattern_error(Krow, Kcol)
+                # Advance the local matrix row pointer
+                ci += 1
+            end
+        end
+        # Make sure that remaining entries in this column of the local matrix are all zero
+        for i in ci:maxlookups
+            if !iszero(Ke[Kerow, permutation[i]])
+                Ferrite._missing_sparsity_pattern_error(Krow, sorteddofs[i])
+            end
+        end
+        current_row += 1
+    end
+end
+
+function Ferrite.zero_out_rows!(K::SparseMatrixCSR, ch::ConstraintHandler) # can be removed in 0.7 with #24711 merged
+    @debug @assert issorted(ch.prescribed_dofs)
+    for row in ch.prescribed_dofs
+        r = nzrange(K, row)
+        K.nzval[r] .= 0.0
+    end
+end
+
+function Ferrite.zero_out_columns!(K::SparseMatrixCSR, ch::ConstraintHandler)
+    colval = K.colval
+    nzval = K.nzval
+    @inbounds for i in eachindex(colval, nzval)
+        if haskey(ch.dofmapping, colval[i])
+            nzval[i] = 0
+        end
+    end
+end
+
+
+function Ferrite.allocate_matrix(::Type{SparseMatrixCSR}, sp::AbstractSparsityPattern)
+    return allocate_matrix(SparseMatrixCSR{1,Float64,Int}, sp)
+end
+
+function Ferrite.allocate_matrix(::Type{MatrixType}, sp::AbstractSparsityPattern) where {Bi, Tv, Ti, MatrixType <: SparseMatrixCSR{Bi, Tv, Ti}}
+    # Allocate CSC first ...
+    K = allocate_matrix(SparseMatrixCSC{Tv, Ti}, sp)
+    # ... and transform to SparseMatrixCSR
+    return SparseMatrixCSR{Bi}(transpose(sparse(transpose(K))))
+end
+
+end
diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl
index aaba4cfd7d..f08a136b48 100644
--- a/src/Dofs/ConstraintHandler.jl
+++ b/src/Dofs/ConstraintHandler.jl
@@ -495,7 +495,7 @@ end
 
 """
 
-    apply!(K::SparseMatrixCSC, rhs::AbstractVector, ch::ConstraintHandler)
+    apply!(K::AbstractSparseMatrix, rhs::AbstractVector, ch::ConstraintHandler)
 
 Adjust the matrix `K` and right hand side `rhs` to account for the Dirichlet boundary
 conditions specified in `ch` such that `K \\ rhs` gives the expected solution.
@@ -587,11 +587,11 @@ function _apply_v(v::AbstractVector, ch::ConstraintHandler, apply_zero::Bool)
     return v
 end
 
-function apply!(K::Union{SparseMatrixCSC,Symmetric}, ch::ConstraintHandler)
+function apply!(K::Union{AbstractSparseMatrix,Symmetric}, ch::ConstraintHandler)
     apply!(K, eltype(K)[], ch, true)
 end
 
-function apply_zero!(K::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler)
+function apply_zero!(K::Union{AbstractSparseMatrix,Symmetric}, f::AbstractVector, ch::ConstraintHandler)
     apply!(K, f, ch, true)
 end
 
@@ -600,7 +600,7 @@ end
 const APPLY_TRANSPOSE = ApplyStrategy.Transpose
 const APPLY_INPLACE = ApplyStrategy.Inplace
 
-function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false;
+function apply!(KK::Union{AbstractSparseMatrix,Symmetric{<:Any,<:AbstractSparseMatrix}}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false;
                 strategy::ApplyStrategy.T=ApplyStrategy.Transpose)
     @assert isclosed(ch)
     sym = isa(KK, Symmetric)
@@ -612,43 +612,14 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con
     m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver
 
     # Add inhomogeneities to f: (f - K * ch.inhomogeneities)
-    if !applyzero
-        @inbounds for i in 1:length(ch.inhomogeneities)
-            d = ch.prescribed_dofs[i]
-            v = ch.inhomogeneities[i]
-            if v != 0
-                for j in nzrange(K, d)
-                    r = K.rowval[j]
-                    sym && r > d && break # don't look below diagonal
-                    f[r] -= v * K.nzval[j]
-                end
-            end
-        end
-        if sym
-            # In the symmetric case, for a constrained dof `d`, we handle the contribution
-            # from `K[1:d, d]` in the loop above, but we are still missing the contribution
-            # from `K[(d+1):size(K,1), d]`. These values are not stored, but since the
-            # matrix is symmetric we can instead use `K[d, (d+1):size(K,1)]`. Looping over
-            # rows is slow, so loop over all columns again, and check if the row is a
-            # constrained row.
-            @inbounds for col in 1:size(K, 2)
-                for ri in nzrange(K, col)
-                    row = K.rowval[ri]
-                    row >= col && break
-                    if (i = get(ch.dofmapping, row, 0); i != 0)
-                        f[col] -= ch.inhomogeneities[i] * K.nzval[ri]
-                    end
-                end
-            end
-        end
-    end
+    !applyzero && add_inhomogeneities!(KK, f, ch.inhomogeneities, ch.prescribed_dofs, ch.dofmapping)
 
-    # Condense K (C' * K * C) and f (C' * f)
+    # Condense K := (C' * K * C) and f := (C' * f)
     _condense!(K, f, ch.dofcoefficients, ch.dofmapping, sym)
 
     # Remove constrained dofs from the matrix
-    zero_out_columns!(K, ch.prescribed_dofs)
-    zero_out_rows!(K, ch.dofmapping)
+    zero_out_columns!(K, ch)
+    zero_out_rows!(K, ch)
 
     # Add meandiag to constraint dofs
     @inbounds for i in 1:length(ch.inhomogeneities)
@@ -661,6 +632,50 @@ function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::Con
     end
 end
 
+"""
+    add_inhomogeneities!(K, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping::Dict{Int,Int})
+
+Compute "f -= K*inhomogeneities".
+By default this is a generic version via SpMSpV kernel.
+"""
+function add_inhomogeneities!(K, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping)
+    f .-= K*sparsevec(prescribed_dofs, inhomogeneities, size(K,2))
+end
+
+# Optimized version for SparseMatrixCSC
+add_inhomogeneities!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K, f, inhomogeneities, prescribed_dofs, dofmapping, false)
+add_inhomogeneities!(K::Symmetric{<:Any,<:SparseMatrixCSC}, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping) = add_inhomogeneities_csc!(K.data, f, inhomogeneities, prescribed_dofs, dofmapping, true)
+function add_inhomogeneities_csc!(K::SparseMatrixCSC, f::AbstractVector, inhomogeneities::AbstractVector, prescribed_dofs::AbstractVector{<:Integer}, dofmapping, sym::Bool)
+    @inbounds for i in 1:length(inhomogeneities)
+        d = prescribed_dofs[i]
+        v = inhomogeneities[i]
+        if v != 0
+            for j in nzrange(K, d)
+                r = K.rowval[j]
+                sym && r > d && break # don't look below diagonal
+                f[r] -= v * K.nzval[j]
+            end
+        end
+    end
+    if sym
+        # In the symmetric case, for a constrained dof `d`, we handle the contribution
+        # from `K[1:d, d]` in the loop above, but we are still missing the contribution
+        # from `K[(d+1):size(K,1), d]`. These values are not stored, but since the
+        # matrix is symmetric we can instead use `K[d, (d+1):size(K,1)]`. Looping over
+        # rows is slow, so loop over all columns again, and check if the row is a
+        # constrained row.
+        @inbounds for col in 1:size(K, 2)
+            for ri in nzrange(K, col)
+                row = K.rowval[ri]
+                row >= col && break
+                if (i = get(dofmapping, row, 0); i != 0)
+                    f[col] -= inhomogeneities[i] * K.nzval[ri]
+                end
+            end
+        end
+    end
+end
+
 # Fetch dof coefficients for a dof prescribed by an affine constraint. Return nothing if the
 # dof is not prescribed, or prescribed by DBC.
 @inline function coefficients_for_dof(dofmapping, dofcoeffs, dof)
@@ -669,6 +684,13 @@ end
     return dofcoeffs[idx]
 end
 
+
+function _condense!(K::AbstractSparseMatrix, f::AbstractVector, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, dofmapping::Dict{Int,Int}, sym::Bool=false) where T
+    # Return early if there are no non-trivial affine constraints
+    any(i -> !(i === nothing || isempty(i)), dofcoefficients) || return
+    error("condensation of ::$(typeof(K)) matrix not supported")
+end
+
 # Condenses K and f: C'*K*C, C'*f, in-place assuming the sparsity pattern is correct
 function _condense!(K::SparseMatrixCSC, f::AbstractVector, dofcoefficients::Vector{Union{Nothing, DofCoefficients{T}}}, dofmapping::Dict{Int,Int}, sym::Bool=false) where T
 
@@ -779,19 +801,19 @@ function create_constraint_matrix(ch::ConstraintHandler{dh,T}) where {dh,T}
 end
 
 # columns need to be stored entries, this is not checked
-function zero_out_columns!(K, dofs::Vector{Int}) # can be removed in 0.7 with #24711 merged
-    @debug @assert issorted(dofs)
-    for col in dofs
+function zero_out_columns!(K::SparseArrays.AbstractSparseMatrixCSC, ch::ConstraintHandler) # can be removed in 0.7 with #24711 merged
+    @debug @assert issorted(ch.prescribed_dofs)
+    for col in ch.prescribed_dofs
         r = nzrange(K, col)
         K.nzval[r] .= 0.0
     end
 end
 
-function zero_out_rows!(K, dofmapping)
+function zero_out_rows!(K::SparseArrays.AbstractSparseMatrixCSC, ch::ConstraintHandler)
     rowval = K.rowval
     nzval = K.nzval
     @inbounds for i in eachindex(rowval, nzval)
-        if haskey(dofmapping, rowval[i])
+        if haskey(ch.dofmapping, rowval[i])
             nzval[i] = 0
         end
     end
diff --git a/src/Ferrite.jl b/src/Ferrite.jl
index ca5ef421e5..86fb1c6c94 100644
--- a/src/Ferrite.jl
+++ b/src/Ferrite.jl
@@ -14,7 +14,7 @@ using NearestNeighbors:
 using OrderedCollections:
     OrderedSet
 using SparseArrays:
-    SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse
+    AbstractSparseMatrix, SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals, sparse, sparsevec
 using StaticArrays:
     StaticArrays, MArray, MMatrix, SArray, SMatrix, SVector
 using WriteVTK:
diff --git a/src/arrayutils.jl b/src/arrayutils.jl
index 2566ee338d..0e74e002c8 100644
--- a/src/arrayutils.jl
+++ b/src/arrayutils.jl
@@ -77,11 +77,11 @@ function addindex!(A::SparseMatrixCSC{Tv}, v::Tv, i::Int, j::Int) where Tv
     end
 end
 
-function fillzero!(A::SparseMatrixCSC{T}) where T
+function fillzero!(A::AbstractSparseMatrix{T}) where T
     fill!(nonzeros(A), zero(T))
     return A
 end
-function fillzero!(A::Symmetric{T,<:SparseMatrixCSC}) where T
+function fillzero!(A::Symmetric{T,<:AbstractSparseMatrix}) where T
     fillzero!(A.data)
     return A
 end
diff --git a/src/assembler.jl b/src/assembler.jl
index f6ad0448eb..ef3f5a3f8a 100644
--- a/src/assembler.jl
+++ b/src/assembler.jl
@@ -99,20 +99,14 @@ Return a reference to the underlying matrix/vector of the assembler.
 """
 matrix_handle, vector_handle
 
-struct AssemblerSparsityPattern{Tv,Ti} <: AbstractSparseAssembler
-    K::SparseMatrixCSC{Tv,Ti}
-    f::Vector{Tv}
-    permutation::Vector{Int}
-    sorteddofs::Vector{Int}
-end
-struct AssemblerSymmetricSparsityPattern{Tv,Ti} <: AbstractSparseAssembler
-    K::Symmetric{Tv,SparseMatrixCSC{Tv,Ti}}
+struct AssemblerSparsityPattern{Tv, MT <: Union{AbstractSparseMatrix{Tv}, Symmetric{Tv,<:AbstractSparseMatrix{Tv}}}} <: AbstractSparseAssembler
+    K::MT
     f::Vector{Tv}
     permutation::Vector{Int}
     sorteddofs::Vector{Int}
 end
 
-function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern,AssemblerSymmetricSparsityPattern})
+function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern})
     print(io, typeof(a), " for assembling into:\n - ")
     summary(io, a.K)
     if !isempty(a.f)
@@ -121,38 +115,28 @@ function Base.show(io::IO, ::MIME"text/plain", a::Union{AssemblerSparsityPattern
     end
 end
 
-matrix_handle(a::AssemblerSparsityPattern) = a.K
-matrix_handle(a::AssemblerSymmetricSparsityPattern) = a.K.data
-vector_handle(a::Union{AssemblerSparsityPattern, AssemblerSymmetricSparsityPattern}) = a.f
+matrix_handle(a::AssemblerSparsityPattern{<:Any, <:AbstractSparseMatrix}) = a.K
+matrix_handle(a::AssemblerSparsityPattern{<:Any, <:Symmetric}) = a.K.data
+vector_handle(a::AssemblerSparsityPattern) = a.f
 
 """
-    start_assemble(K::SparseMatrixCSC;            fillzero::Bool=true) -> AssemblerSparsityPattern
-    start_assemble(K::SparseMatrixCSC, f::Vector; fillzero::Bool=true) -> AssemblerSparsityPattern
+    start_assemble(K::AbstractSparseMatrix;            fillzero::Bool=true) -> AssemblerSparsityPattern
+    start_assemble(K::AbstractSparseMatrix, f::Vector; fillzero::Bool=true) -> AssemblerSparsityPattern
 
 Create a `AssemblerSparsityPattern` from the matrix `K` and optional vector `f`.
 
-    start_assemble(K::Symmetric{SparseMatrixCSC};                 fillzero::Bool=true) -> AssemblerSymmetricSparsityPattern
-    start_assemble(K::Symmetric{SparseMatrixCSC}, f::Vector=Td[]; fillzero::Bool=true) -> AssemblerSymmetricSparsityPattern
-
-Create a `AssemblerSymmetricSparsityPattern` from the matrix `K` and optional vector `f`.
-
-`AssemblerSparsityPattern` and `AssemblerSymmetricSparsityPattern` allocate workspace
-necessary for efficient matrix assembly. To assemble the contribution from an element, use
-[`assemble!`](@ref).
+`AssemblerSparsityPattern` allocate workspace necessary for efficient matrix assembly. To assemble
+ the contribution from an element, use [`assemble!`](@ref).
 
 The keyword argument `fillzero` can be set to `false` if `K` and `f` should not be zeroed
 out, but instead keep their current values.
 """
-start_assemble(K::Union{SparseMatrixCSC, Symmetric{<:Any,SparseMatrixCSC}}, f::Vector; fillzero::Bool)
+start_assemble(K::AbstractSparseMatrix, f::Vector; fillzero::Bool)
 
-function start_assemble(K::SparseMatrixCSC{T}, f::Vector=T[]; fillzero::Bool=true) where {T}
+function start_assemble(K::Union{AbstractSparseMatrix{T}, Symmetric{T,<: AbstractSparseMatrix{T}}}, f::Vector=T[]; fillzero::Bool=true) where {T}
     fillzero && (fillzero!(K); fillzero!(f))
     return AssemblerSparsityPattern(K, f, Int[], Int[])
 end
-function start_assemble(K::Symmetric{T,<:SparseMatrixCSC}, f::Vector=T[]; fillzero::Bool=true) where T
-    fillzero && (fillzero!(K); fillzero!(f))
-    return AssemblerSymmetricSparsityPattern(K, f, Int[], Int[])
-end
 
 """
     assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix)
@@ -172,15 +156,15 @@ end
 @propagate_inbounds function assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, fe::AbstractVector, Ke::AbstractMatrix)
     assemble!(A, dofs, Ke, fe)
 end
-@propagate_inbounds function assemble!(A::AssemblerSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
+@propagate_inbounds function assemble!(A::AssemblerSparsityPattern{<:Any, <:AbstractSparseMatrix}, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
     _assemble!(A, dofs, Ke, fe, false)
 end
-@propagate_inbounds function assemble!(A::AssemblerSymmetricSparsityPattern, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
+@propagate_inbounds function assemble!(A::AssemblerSparsityPattern{<:Any, <:Symmetric{<:Any,<:AbstractSparseMatrix}}, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector)
     _assemble!(A, dofs, Ke, fe, true)
 end
 
+# Main entry point for the CPU
 @propagate_inbounds function _assemble!(A::AbstractSparseAssembler, dofs::AbstractVector{Int}, Ke::AbstractMatrix, fe::AbstractVector, sym::Bool)
-    ld = length(dofs)
     @boundscheck checkbounds(Ke, keys(dofs), keys(dofs))
     if length(fe) != 0
         @boundscheck checkbounds(fe, keys(dofs))
@@ -192,12 +176,18 @@ end
     permutation = A.permutation
     sorteddofs = A.sorteddofs
     @boundscheck checkbounds(K, dofs, dofs)
+    ld = length(dofs)
     resize!(permutation, ld)
     resize!(sorteddofs, ld)
     copyto!(sorteddofs, dofs)
     sortperm2!(sorteddofs, permutation)
 
+    _assemble_inner!(K, Ke, dofs, sorteddofs, permutation, sym)
+end
+
+@propagate_inbounds function _assemble_inner!(K::SparseMatrixCSC, Ke::AbstractMatrix, dofs::AbstractVector, sorteddofs::AbstractVector, permutation::AbstractVector, sym::Bool)
     current_col = 1
+    ld = length(dofs)
     @inbounds for Kcol in sorteddofs
         maxlookups = sym ? current_col : ld
         Kecol = permutation[current_col]
diff --git a/test/runtests.jl b/test/runtests.jl
index 3a5e35e708..d09f24a385 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -58,6 +58,7 @@ include("test_apply_analytical.jl")
 include("PoolAllocator.jl")
 include("test_deprecations.jl")
 HAS_EXTENSIONS && include("blockarrays.jl")
+HAS_EXTENSIONS && include("test_assembler_extensions.jl")
 include("test_examples.jl")
 
 @test all(x -> isdefined(Ferrite, x), names(Ferrite))  # Test that all exported symbols are defined
diff --git a/test/test_assembler_extensions.jl b/test/test_assembler_extensions.jl
new file mode 100644
index 0000000000..feb2e2568c
--- /dev/null
+++ b/test/test_assembler_extensions.jl
@@ -0,0 +1,107 @@
+import SparseMatricesCSR: SparseMatrixCSR, sparsecsr
+using SparseArrays, LinearAlgebra
+
+@testset "SparseMatricesCSR extension" begin
+
+@testset "apply!(::SparseMatrixCSR,...)" begin
+    # Specifically this test that values below the diagonal of K2::Symmetric aren't touched
+    # and that the missing values are instead taken from above the diagonal.
+    grid = generate_grid(Line, (2,))
+    dh = DofHandler(grid)
+    add!(dh, :u, Lagrange{RefLine,1}())
+    close!(dh)
+    ch = ConstraintHandler(dh)
+    add!(ch, Dirichlet(:u, getfacetset(grid, "left"), x -> 1))
+    close!(ch)
+    K0 = sparse(rand(3, 3))
+    K0 = K0'*K0
+    K1 = SparseMatrixCSR(transpose(copy(K0)))
+    K2 = Symmetric(SparseMatrixCSR(transpose(copy(K0))))
+    @test K0 == K1
+    @test K1 == K2
+    sol = [1.0, 2.0, 3.0]
+    f0 = K0*sol
+    f1 = K1*sol
+    f2 = K2*sol
+    apply!(K0, f0, ch)
+    apply!(K1, f1, ch)
+    apply!(K2, f2, ch)
+    @test K0 == K1
+    @test K1 == K2
+    @test f0 == f1
+    @test f1 == f2
+    # Error for affine constraints
+    ch = ConstraintHandler(dh)
+    add!(ch, AffineConstraint(1, [3 => 1.0], 1.0))
+    close!(ch)
+    @test_throws ErrorException("condensation of ::SparseMatrixCSR{1, Float64, Int64} matrix not supported") apply!(K2, f2, ch)
+end
+
+@testset "assembly integration" begin
+    # Setup simple problem
+    grid = generate_grid(Line, (2,))
+    dh = DofHandler(grid)
+    add!(dh, :u, Lagrange{RefLine, 1}())
+    close!(dh)
+
+    # Check if the matrix is correctly allocated
+    K = allocate_matrix(SparseMatrixCSR, dh)
+    I = [1,1,2,2,2,3,3]
+    J = [1,2,1,2,3,2,3]
+    V = zeros(7)
+    @test K == sparsecsr(I,J,V)
+    f = zeros(3)
+
+    # Check that incuding the ch doesnot mess up the pattern
+    ch = ConstraintHandler(dh)
+    add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 1))
+    close!(ch)
+    @test K == allocate_matrix(SparseMatrixCSR, dh, ch)
+
+    # Check if assembly works
+    assembler = start_assemble(K, f)
+    ke = [-1.0 1.0; 2.0 -1.0]
+    fe = [1.0,2.0]
+    assemble!(assembler, [1,2], ke,fe)
+    assemble!(assembler, [3,2], ke,fe)
+    I = [1,1,2,2,2,3,3]
+    J = [1,2,1,2,3,2,3]
+    V = [-1.0,1.0,2.0,-2.0,2.0,1.0,-1.0]
+    @test K ≈ sparsecsr(I,J,V)
+    @test f ≈ [1.0,4.0,1.0]
+
+    # Check if constraint handler integration works
+    apply!(K,f,ch)
+    I = [1,1,2,2,2,3,3]
+    J = [1,2,1,2,3,2,3]
+    V = [4/3,0.0,0.0,-2.0,2.0,1.0,-1.0]
+    @test K ≈ sparsecsr(I,J,V)
+    @test f ≈ [4/3,2.0,1.0]
+
+    # Check if coupling works
+    grid = generate_grid(Quadrilateral, (2,2))
+    ip = Lagrange{RefQuadrilateral,1}()
+    dh = DofHandler(grid)
+    add!(dh, :u, ip)
+    add!(dh, :v, ip)
+    close!(dh)
+
+    Ke_zeros = zeros(ndofs_per_cell(dh,1), ndofs_per_cell(dh,1))
+    Ke_rand = rand(ndofs_per_cell(dh,1), ndofs_per_cell(dh,1))
+    dofs = celldofs(dh,1)
+
+    for c1 ∈ [true, false], c2 ∈ [true, false], c3 ∈ [true, false], c4 ∈ [true, false]
+        coupling = [c1; c2;; c3; c4]
+        K = allocate_matrix(SparseMatrixCSR, dh; coupling)
+        a = start_assemble(K)
+        assemble!(a, dofs, Ke_zeros)
+        if all(coupling)
+            assemble!(a, dofs, Ke_rand)
+            @test Ke_rand ≈ K[dofs,dofs]
+        else
+            @test_throws ErrorException assemble!(a, dofs, Ke_rand)
+        end
+    end
+end
+
+end