Use manual construction of sparse matrix

Author: pjaap

Project.toml

@@ -13,6 +13,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExtendableFEMBase = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -38,6 +39,7 @@ ExplicitImports = "1"
 ExtendableFEMBase = "1.3.0"
 ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
+FillArrays = "1.13.0"
 ForwardDiff = "0.10.35,1"
 GridVisualize = "1.8.1"
 IncompleteLU = "0.2.1"

src/ExtendableFEM.jl

@@ -5,8 +5,8 @@ $(read(joinpath(@__DIR__, "..", "README.md"), String))
 """
 module ExtendableFEM
 
-using BlockArrays: BlockMatrix, BlockVector, Block
 using ChunkSplitters: chunks
+using BlockArrays: BlockMatrix, BlockVector, Block, blocks, axes
 using CommonSolve: CommonSolve
 using DiffResults: DiffResults
 using DocStringExtensions: DocStringExtensions, TYPEDEF, TYPEDSIGNATURES
@@ -61,6 +61,7 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
 using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
     MTExtendableSparseMatrixCSC, findindex,
     rawupdateindex!
+using FillArrays: Zeros
 using ForwardDiff: ForwardDiff
 using GridVisualize: GridVisualize, GridVisualizer, gridplot!, reveal, save,
     scalarplot!, vectorplot!

src/common_restrictions/coupled_dofs_restriction.jl

@@ -32,7 +32,7 @@ function assemble!(R::CoupledDofsRestriction, SC; kwargs...)
     B = (R.coupling_matrix - LinearAlgebra.I)[:, unique_cols]
 
     R.parameters[:matrix] = B
-    return R.parameters[:rhs] = Zeros(length(unique_cols))
+    R.parameters[:rhs] = Zeros(length(unique_cols))
 
     return nothing
 end

src/solvers.jl

@@ -250,51 +250,72 @@ function solve_linear_system!(A, b, sol, soltemp, residual, linsolve, unknowns,
             b_unrestricted = residual.entries
         end
 
-        @timeit timer "LM restrictions" begin
-            ## add possible Lagrange restrictions
-            @timeit timer "prepare" begin
-                restriction_matrices = [length(freedofs) > 0 ? re.parameters[:matrix][freedofs, :] : re.parameters[:matrix] for re in PD.restrictions ]
-                restriction_rhs = [length(freedofs) > 0 ? re.parameters[:rhs][freedofs] : re.parameters[:rhs] for re in PD.restrictions ]
-
-                ## we need to add the (initial) solution to the rhs, since we work with the residual equation
-                for (B, rhs) in zip(restriction_matrices, restriction_rhs)
-                    rhs .+= B'sol.entries
-                end
-            end
+        if length(PD.restrictions) == 0
+            linsolve.A = A_unrestricted
+            linsolve.b = b_unrestricted
+        else
 
+            @timeit timer "LM restrictions" begin
+                ## add possible Lagrange restrictions
+                @timeit timer "prepare" begin
+                    restriction_matrices = [length(freedofs) > 0 ? re.parameters[:matrix][freedofs, :] : re.parameters[:matrix] for re in PD.restrictions ]
+                    restriction_rhs = [length(freedofs) > 0 ? re.parameters[:rhs][freedofs] : re.parameters[:rhs] for re in PD.restrictions ]
 
-            @timeit timer "compute blocks" begin
-                # block sizes for the block matrix
-                block_sizes = [size(A_unrestricted, 2), [ size(B, 2) for B in restriction_matrices ]...]
+                    ## we need to add the (initial) solution to the rhs, since we work with the residual equation
+                    for (B, rhs) in zip(restriction_matrices, restriction_rhs)
+                        rhs .+= B'sol.entries
+                    end
+                end
+
+                @timeit timer "set blocks" begin
+                    # block sizes for the block matrix
+                    block_sizes = [size(A_unrestricted, 2), [ size(B, 2) for B in restriction_matrices ]...]
 
-                total_size = sum(block_sizes)
-                Tv = eltype(A_unrestricted)
+                    total_size = sum(block_sizes)
+                    Tv = eltype(A_unrestricted)
 
-                ## create block matrix
-                A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
-                A_block[Block(1, 1)] = A_unrestricted
+                    ## create block matrix
+                    A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
+                    A_block[Block(1, 1)] = A_unrestricted
 
-                b_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                b_block[Block(1)] = b_unrestricted
+                    b_block = BlockVector(zeros(Tv, total_size), block_sizes)
+                    b_block[Block(1)] = b_unrestricted
 
-                u_unrestricted = linsolve.u
-                u_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                u_block[Block(1)] = u_unrestricted
+                    u_unrestricted = linsolve.u
+                    u_block = BlockVector(zeros(Tv, total_size), block_sizes)
+                    u_block[Block(1)] = u_unrestricted
 
-                for i in eachindex(PD.restrictions)
-                    A_block[Block(1, i + 1)] = restriction_matrices[i]
-                    A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
-                    b_block[Block(i + 1)] = restriction_rhs[i]
+                    for i in eachindex(PD.restrictions)
+                        A_block[Block(1, i + 1)] = restriction_matrices[i]
+                        A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
+                        b_block[Block(i + 1)] = restriction_rhs[i]
 
+                    end
                 end
-            end
 
-            @timeit timer "convert" begin
+                @timeit timer "convert" begin
+
+                    linsolve.b = Vector(b_block) # convert to dense vector
+                    linsolve.u = Vector(u_block) # convert to dense vector
+
+                    # linsolve.A = sparse(A_block) # convert to CSC Matrix is very slow https://github.com/JuliaArrays/BlockArrays.jl/issues/78
+                    # do it manually:
+                    A_flat = spzeros(size(A_block))
+
+                    lasts_i = [0, axes(A_block)[1].lasts...] # add zero
+                    lasts_j = [0, axes(A_block)[2].lasts...] # add zero
 
-                linsolve.A = sparse(A_block) # convert to CSC Matrix
-                linsolve.b = Vector(b_block) # convert to dense vector
-                linsolve.u = Vector(u_block) # convert to dense vector
+                    (ni, nj) = size(blocks(A_block))
 
+                    for i in 1:ni, j in 1:nj
+                        range_i = (lasts_i[i] + 1):lasts_i[i + 1]
+                        range_j = (lasts_j[j] + 1):lasts_j[j + 1]
+
+                        # write each block directly in the resulting matrix
+                        A_flat[range_i, range_j] = A_block[Block(i, j)]
+                    end
+                    linsolve.A = A_flat
+                end
             end
         end