Repository: sympy/sympy

Issue:
Error in inverse of BlockMatrix
The inverse computation of a square `BlockMatrix` with more than 3x3 blocks sometimes provides a incorrect solution. For example:

```python
import sympy as sym
from IPython.display import display

n  = 10
I  = sym.Identity(n)
O  = sym.ZeroMatrix(n, n)

A  = sym.BlockMatrix([[  I,  O,  O,  O ],
                      [  O,  I,  O,  O ],
                      [  O,  O,  I,  O ],
                      [  I,  O,  O,  I ]])
Ainv = sym.block_collapse(A.inv())
Am1A = sym.block_collapse(Ainv @ A)
display(Am1A)
```
```
>>> Matrix([
[  I, 0, 0, 0],
[  0, I, 0, 0],
[  0, 0, I, 0],
[2*I, 0, 0, I]])
```


```python
Id = sym.BlockMatrix([[  I,  O,  O,  O ],
                      [  O,  I,  O,  O ],
                      [  O,  O,  I,  O ],
                      [  O,  O,  O,  I ]])
print('inv(A)*A == Id: %s' % (Am1A == Id))
```
```
>>> False
```

Notice that the correct solution should be the following `B` matrix:
```python
B  = sym.BlockMatrix([[  I,  O,  O,  O ],
                      [  O,  I,  O,  O ],
                      [  O,  O,  I,  O ],
                      [ -I,  O,  O,  I ]])
BA = sym.block_collapse(B @ A)
print('B*A  == Id: %s' % (BA == Id))
```
```
>>> True
```

*Remark*: the same issue is present with any of the following methods:
```python
Ainv = sym.block_collapse(A**-1)
Ainv = sym.block_collapse(A.inverse())
Ainv = sym.block_collapse(sym.Inverse(A))
```

Discussion:

Patch:
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 193e429cc8..d7809926ed 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -777,7 +777,7 @@ def bc_dist(expr):
         new_B = [
             [factor * B[i, j] for j in range(B.cols)] for i in range(B.rows)]
         return BlockMatrix(new_B)
-    return unpacked
+    return expr
 
 
 def bc_matmul(expr):


Test patch:
diff --git a/sympy/matrices/expressions/tests/test_blockmatrix.py b/sympy/matrices/expressions/tests/test_blockmatrix.py
index b196fd7c8a..f52295b6b5 100644
--- a/sympy/matrices/expressions/tests/test_blockmatrix.py
+++ b/sympy/matrices/expressions/tests/test_blockmatrix.py
@@ -415,3 +415,18 @@ def test_block_lu_decomposition():
     #LU decomposition
     L, U = X.LUdecomposition()
     assert block_collapse(L*U) == X
+
+def test_issue_21866():
+    n  = 10
+    I  = Identity(n)
+    O  = ZeroMatrix(n, n)
+    A  = BlockMatrix([[  I,  O,  O,  O ],
+                      [  O,  I,  O,  O ],
+                      [  O,  O,  I,  O ],
+                      [  I,  O,  O,  I ]])
+    Ainv = block_collapse(A.inv())
+    AinvT = BlockMatrix([[  I,  O,  O,  O ],
+                      [  O,  I,  O,  O ],
+                      [  O,  O,  I,  O ],
+                      [  -I,  O,  O,  I ]])
+    assert Ainv == AinvT