From 7127e3f76ad70403eac7f98e79e2ae07cd147457 Mon Sep 17 00:00:00 2001
From: Theo Galy-Fajou <theo.galyfajou@gmail.com>
Date: Wed, 24 Nov 2021 15:17:46 +0100
Subject: [PATCH] Implemented suggestions

---
 src/linalg.jl  |  9 +++------
 test/linalg.jl | 11 +++++++++--
 2 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/src/linalg.jl b/src/linalg.jl
index 241d66a..82d56b8 100644
--- a/src/linalg.jl
+++ b/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)
diff --git a/test/linalg.jl b/test/linalg.jl
index aafd705..457fd0c 100644
--- a/test/linalg.jl
+++ b/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