Repository: sympy/sympy

Issue:
LUsolve raises IndexError even though system has a unique solution
I encountered (essentially) this:
```python
>>> Matrix([[2, 1], [1, 0], [1, 0]]).LUsolve(Matrix([3, 1, 1]))
...
IndexError: Index out of range: a[2]
```
The reason the matrix was not in a canonical format is because it is generated programmatically. Manually, removing the superfluous row allows the system to be solved:
```python
>>> Matrix([[2, 1], [1, 0]]).LUsolve(Matrix([3, 1]))
⎡1⎤
⎢ ⎥
⎣1⎦
```
What do you think, should we allow LUsolve to handle this kind of "apparently" overdetermined systems?

Discussion:

Patch:
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 56329a9a6b..31dfc844d9 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -3589,14 +3589,24 @@ def LUsolve(self, rhs, iszerofunc=_iszero):
             raise ShapeError(
                 "`self` and `rhs` must have the same number of rows.")
 
+        m = self.rows
+        n = self.cols
+        if m < n:
+            raise NotImplementedError("Underdetermined systems not supported.")
         A, perm = self.LUdecomposition_Simple(iszerofunc=_iszero)
-        n = self.rows
         b = rhs.permute_rows(perm).as_mutable()
         # forward substitution, all diag entries are scaled to 1
-        for i in range(n):
-            for j in range(i):
+        for i in range(m):
+            for j in range(min(i, n)):
                 scale = A[i, j]
                 b.zip_row_op(i, j, lambda x, y: x - y * scale)
+        # consistency check for overdetermined systems
+        if m > n:
+            for i in range(n, m):
+                for j in range(b.cols):
+                    if not iszerofunc(b[i, j]):
+                        raise ValueError("The system is inconsistent.")
+            b = b[0:n, :]   # truncate zero rows if consistent
         # backward substitution
         for i in range(n - 1, -1, -1):
             for j in range(i + 1, n):


Test patch:
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index 35a4ba0a28..a82f26e4e2 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -640,6 +640,25 @@ def test_LUsolve():
     b = A*x
     soln = A.LUsolve(b)
     assert soln == x
+    A = Matrix([[2, 1], [1, 0], [1, 0]])   # issue 14548
+    b = Matrix([3, 1, 1])
+    assert A.LUsolve(b) == Matrix([1, 1])
+    b = Matrix([3, 1, 2])                  # inconsistent
+    raises(ValueError, lambda: A.LUsolve(b))
+    A = Matrix([[0, -1, 2],
+                [5, 10, 7],
+                [8,  3, 4],
+                [2, 3, 5],
+                [3, 6, 2],
+                [8, 3, 6]])
+    x = Matrix([2, 1, -4])
+    b = A*x
+    soln = A.LUsolve(b)
+    assert soln == x
+    A = Matrix([[0, -1, 2], [5, 10, 7]])  # underdetermined
+    x = Matrix([-1, 2, 0])
+    b = A*x
+    raises(NotImplementedError, lambda: A.LUsolve(b))
 
 
 def test_QRsolve():