Update get_op_output_unit with new type invdiv.
It outputs the product of the following units over the first one in the args list.

Update tests with values & np.dot(A, x) == b.
Where x = np.linalg.solve(A, b)

Closes #1246

1224: Updated linter r=hgrecco a=hgrecco

- [x] Closes #1223 (insert issue number)
- [x] Executed ``pre-commit run --all-files`` with no errors
- [ ] The change is fully covered by automated unit tests
- [ ] Documented in docs/ as appropriate
- [ ] Added an entry to the CHANGES file


Co-authored-by: Hernan <hernan.grecco@gmail.com>

1120: revise the unit computation for np.prod r=hgrecco a=keewis

#1087 left out the case where `axis` and `where` are specified and also used the size of the result to compute the output unit if only `axis` was specified. This changes that to use `a.shape[axis]` instead and also implements the support for both `axis` and `where` by broadcasting `where` against the array, applying `np.sum` along `axis` and using the only one unique value (`0` doesn't count) as an exponent. In case there's more than that, it will try to cast to `dimensionless`.

I'm not quite sure if using `np.broadcast_arrays` is the best way to get the exponents, though.

Edit: **Todo**: make the error message easier to understand

- [x] Closes #867
- [x] Executed ``black -t py36 . && isort -rc . && flake8`` with no errors
- [x] The change is fully covered by automated unit tests
- [ ] Documented in docs/ as appropriate
- [ ] Added an entry to the CHANGES file


Co-authored-by: Keewis <keewis@posteo.de>