Remove BlockArrays dep and handle blocked arrays manually

This and the transpose of the R factor are a huge slow-down for very large systems.
In the end, we copy the data by hand anyway, so do not use BlockArrays at all.

Author: pjaap

Project.toml

@@ -5,7 +5,6 @@ authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.ja
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
@@ -30,7 +29,6 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 [compat]
 ADTypes = "1.16.0"
 Aqua = "0.8"
-BlockArrays = "1.7.0"
 ChunkSplitters = "3.1.2"
 CommonSolve = "0.2"
 DiffResults = "1"

src/ExtendableFEM.jl

@@ -7,7 +7,6 @@ module ExtendableFEM
 
 using ADTypes: ADTypes, KnownJacobianSparsityDetector
 using ChunkSplitters: chunks
-using BlockArrays: BlockMatrix, BlockVector, Block, blocks, axes, mortar
 using CommonSolve: CommonSolve
 using DifferentiationInterface: DifferentiationInterface,
     AutoSparse, AutoForwardDiff, prepare_jacobian

src/solvers.jl

@@ -271,7 +271,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
         end
 
         ## restrictions only involve the blocks coressponding to the unknowns and not necessarily the full sol.entries
-        sol_freedofs = mortar([view(sol[u]) for u in unknowns])
+        sol_freedofs = vcat([view(sol[u]) for u in unknowns]...)
 
         if length(PD.restrictions) == 0
             if linsolve_needs_matrix
@@ -285,7 +285,7 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
 
             # block sizes for the block matrix
             block_sizes = [length(b_unrestricted), [ length(b) for b in restriction_rhss ]...]
-            total_size = sum(block_sizes)
+            block_ends = cumsum(block_sizes)
 
             # total number of additional LM dofs
             nLMs = @views sum(block_sizes[2:end])
@@ -315,14 +315,15 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
                         # we need the inverse column permutation
                         ipcol = invperm(pcol)
 
-                        # the new combined restriction block (add another transpose)
-                        combined_restriction_matrix = R[:, ipcol]'
+                        # the new combined restriction block (add another transpose and convert into CSC matrix)
+                        combined_restriction_matrix = sparse(R[:, ipcol]')
 
                         # the new combined restriction rhs
                         combined_restriction_rhs = Q'combined_restriction_rhs[prow]
 
                         # compress the  column space
                         qr_rank = rank(qr_result)
+                        SC.parameters[:verbosity] > 1 && @info "QR decomposition has rank $qr_rank"
                         @assert norm(combined_restriction_rhs[(qr_rank + 1):end]) ≤ 1.0e-12 * norm(combined_restriction_rhs) "the rhs of the restriction is not in the image"
 
                         combined_restriction_matrix = combined_restriction_matrix[:, 1:qr_rank]
@@ -337,63 +338,36 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
 
                         # update sizes
                         block_sizes = [length(b_unrestricted), length(combined_restriction_rhs)]
-                        total_size = sum(block_sizes)
+                        block_ends = cumsum(block_sizes)
 
                         SC.parameters[:verbosity] > 1 && @info "Created new CompressedRestriction"
                     end
 
-                    A_block = BlockMatrix(spzeros(Tv, total_size, total_size), block_sizes, block_sizes)
-                    A_block[Block(1, 1)] = A_unrestricted
+                    linsolve_A = spzeros(Tv, block_ends[end], block_ends[end])
+                    linsolve_A[1:block_ends[1], 1:block_ends[1]] = A_unrestricted
                 end
 
                 ## 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_rhss)
-                    rhs .-= B'sol_freedofs
+                    # rhs -= B'sol_freedofs
+                    mul!(rhs, B', sol_freedofs, -1.0, 0.0)
                 end
 
-                b_block = BlockVector(zeros(Tv, total_size), block_sizes)
-                b_block[Block(1)] = b_unrestricted
+                linsolve_b = zeros(Tv, block_ends[end])
+                linsolve_b[1:block_ends[1]] = b_unrestricted
 
                 for i in eachindex(restriction_matrices)
                     if linsolve_needs_matrix
-                        A_block[Block(1, i + 1)] = restriction_matrices[i]
-                        A_block[Block(i + 1, 1)] = transpose(restriction_matrices[i])
+                        linsolve_A[1:block_ends[1], (block_ends[i] + 1):block_ends[i + 1]] = restriction_matrices[i]
+                        linsolve_A[(block_ends[i] + 1):block_ends[i + 1], 1:block_ends[1]] = restriction_matrices[i]'
                     end
-                    b_block[Block(i + 1)] = restriction_rhss[i]
+                    linsolve_b[(block_ends[i] + 1):block_ends[i + 1]] = restriction_rhss[i]
 
                 end
-
-                # convert to dense vectors
-                linsolve_b = Vector(b_block)
-
-                if linsolve_needs_matrix
-
-                    # 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_row = [0, axes(A_block)[1].lasts...] # add leading zero
-                    lasts_col = [0, axes(A_block)[2].lasts...] # add leading zero
-
-                    (n_row, n_col) = size(blocks(A_block))
-
-
-                    SC.parameters[:verbosity] > 1 && @info "Assemble flat matrix out of blocks..."
-                    for i in 1:n_row, j in 1:n_col
-                        range_row = (lasts_row[i] + 1):lasts_row[i + 1]
-                        range_col = (lasts_col[j] + 1):lasts_col[j + 1]
-
-                        # write each block directly in the resulting matrix
-                        A_flat[range_row, range_col] = A_block[Block(i, j)]
-                    end
-                    SC.parameters[:verbosity] > 1 && @info "... matrix assembly complete!"
-
-                    linsolve_A = A_flat
-                end
+                SC.parameters[:verbosity] > 1 && @info "done setting blocked matrix/vector"
             end
         end
 
-
         if SC.parameters[:initialized]
             # set/update the linear system
             SC.linsolver.b = linsolve_b
@@ -442,7 +416,6 @@ function solve_linear_system!(A, b, sol, soltemp, residual, unknowns, freedofs,
 
         if length(PD.restrictions) > 0
             # extract all residuals for the restriction blocks
-            block_ends = cumsum(block_sizes)
             restriction_residuals = [norm(residual_flat[(block_ends[i] + 1):block_ends[i + 1]]) for i in 1:(length(block_sizes) - 1) ]
             push!(stats[:restriction_residuals], restriction_residuals)