Implemented suggestions

src/linalg.jl

@@ -153,12 +153,9 @@ end
 function LinearAlgebra.:\(B::BlockDiagonal, vm::AbstractVecOrMat)
  row_i = 1
  # BlockDiagonals with non-square blocks
- all(b -> size(b, 1) == size(b, 2), blocks(B)) || throw(
- ArgumentError(
- "Left division (\\) is not defined for BlockDiagonals with " *
- "non-square blocks (they are singular matrices).",
- ),
- )
+ if !all(is_square, blocks(B))
+ return Matrix(B) \ vm # Fallback on the generic LinearAlgebra method
+ end
  result = similar(vm)
  for block in blocks(B)
  nrow = size(block, 1)

test/linalg.jl

@@ -228,7 +228,14 @@ using Test
  end
  end # SVD
  @testset "Left division" begin
- @test b1 \ a ≈ inv(b1) * a
- @test_throws ArgumentError BlockDiagonal([A, B]) \ rand(rng, 2N + N1 + N2)
+ N1 = 20
+ N2 = 8
+ N3 = 5
+ A = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2)])
+ B = BlockDiagonal([rand(rng, N1, N2), rand(rng, N3, N1)])
+ x = rand(rng, N1 + N2)
+ y = rand(rng, N1 + N3)
+ @test A \ x ≈ inv(A) * x
+ @test B \ y ≈ Matrix(B) \ y
  end
 end