Repository: hgrecco/pint

Issue:
numpy.linalg.solve returns wrong units
`numpy.linalg.solve` returns units that are inverted with respect to what is expected. E.g.

```python
A = np.random.rand(5, 5)*ureg.m
b = np.random.rand(5)*ureg.s
np.linalg.solve(A, b)
```

returns `<Unit('meter / second')>` for `x`, whereas `<Unit('second / meter')>` is expected, since the matrix equation

![equation](https://latex.codecogs.com/svg.latex?A%20%5Ccdot%20x%20%3D%20b)

expressed in units, corresponds to `meter * second / meter = second`.

(Note also, that the wrong behaviour is explicitly tested [here](https://github.com/hgrecco/pint/blob/b36d0a63ee24e7dada0cbf8b0499695ff57d17a8/pint/testsuite/test_numpy.py#L452).)

Discussion:

Patch:
diff --git a/CHANGES b/CHANGES
index b93b83b..cf50d1e 100644
--- a/CHANGES
+++ b/CHANGES
@@ -12,6 +12,7 @@ Pint Changelog
 - UnitsContainer returns false if other is str and cannnot be parsed
   (Issue #1179, thanks rfrowe)
 - Add Github Actions CI. (Issue #1236)
+- Fix numpy.linalg.solve unit output. (Issue #1246)
 
 0.16.1 (2020-09-22)
 -------------------
diff --git a/pint/numpy_func.py b/pint/numpy_func.py
index 1dce044..c335f3d 100644
--- a/pint/numpy_func.py
+++ b/pint/numpy_func.py
@@ -148,6 +148,7 @@ def get_op_output_unit(unit_op, first_input_units, all_args=None, size=None):
     - "sqrt": square root of `first_input_units`
     - "reciprocal": reciprocal of `first_input_units`
     - "size": `first_input_units` raised to the power of `size`
+    - "invdiv": inverse of `div`, product of all following units divided by first argument unit
 
     Parameters
     ----------
@@ -205,7 +206,15 @@ def get_op_output_unit(unit_op, first_input_units, all_args=None, size=None):
         if size is None:
             raise ValueError('size argument must be given when unit_op=="size"')
         result_unit = first_input_units ** size
-
+    elif unit_op == "invdiv":
+        # Start with first arg in numerator, all others in denominator
+        product = getattr(
+            all_args[0], "units", first_input_units._REGISTRY.parse_units("")
+        )
+        for x in all_args[1:]:
+            if hasattr(x, "units"):
+                product /= x.units
+        result_unit = product ** -1
     else:
         raise ValueError("Output unit method {} not understood".format(unit_op))
 
@@ -304,6 +313,7 @@ def implement_func(func_type, func_str, input_units=None, output_unit=None):
             "delta",
             "delta,div",
             "div",
+            "invdiv",
             "variance",
             "square",
             "sqrt",
@@ -878,7 +888,7 @@ for func_str in ["diff", "ediff1d"]:
 for func_str in ["gradient"]:
     implement_func("function", func_str, input_units=None, output_unit="delta,div")
 for func_str in ["linalg.solve"]:
-    implement_func("function", func_str, input_units=None, output_unit="div")
+    implement_func("function", func_str, input_units=None, output_unit="invdiv")
 for func_str in ["var", "nanvar"]:
     implement_func("function", func_str, input_units=None, output_unit="variance")
 


Test patch:
diff --git a/pint/testsuite/test_numpy.py b/pint/testsuite/test_numpy.py
index beea69d..44d4271 100644
--- a/pint/testsuite/test_numpy.py
+++ b/pint/testsuite/test_numpy.py
@@ -449,10 +449,13 @@ class TestNumpyMathematicalFunctions(TestNumpyMethods):
 
     @helpers.requires_array_function_protocol()
     def test_solve(self):
-        helpers.assert_quantity_almost_equal(
-            np.linalg.solve(self.q, [[3], [7]] * self.ureg.s),
-            self.Q_([[1], [1]], "m / s"),
-        )
+        A = self.q
+        b = [[3], [7]] * self.ureg.s
+        x = np.linalg.solve(A, b)
+
+        helpers.assert_quantity_almost_equal(x, self.Q_([[1], [1]], "s / m"))
+
+        helpers.assert_quantity_almost_equal(np.dot(A, x), b)
 
     # Arithmetic operations
     def test_addition_with_scalar(self):