def tensorsolve(a, b, axes=None):
a, wrap = _makearray(a)<EOL>b = asarray(b)<EOL>an = a.ndim<EOL>if axes is not None:<EOL><INDENT>allaxes = list(range(<NUM_LIT:0>, an))<EOL>for k in axes:<EOL><INDENT>allaxes.remove(k)<EOL>allaxes.insert(an, k)<EOL><DEDENT>a = a.transpose(allaxes)<EOL><DEDENT>oldshape = a.shape[-(an-b.ndim):]<EOL>prod = <NUM_LIT:1><EOL>for k in oldshape:<EOL><INDENT>prod *= k<EOL><DEDENT>a = a.reshape(-<NUM_LIT:1>, prod)<EOL>b = b.ravel()<EOL>res = wrap(solve(a, b))<EOL>res.shape = oldshape<EOL>return res<EOL>
Solve the tensor equation ``a x = b`` for x.

It is assumed that all indices of `x` are summed over in the product,
together with the rightmost indices of `a`, as is done in, for example,
``tensordot(a, x, axes=len(b.shape))``.

Parameters
----------
a : array_like
    Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
    the shape of that sub-tensor of `a` consisting of the appropriate
    number of its rightmost indices, and must be such that
   ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
   'square').
b : array_like
    Right-hand tensor, which can be of any shape.
axes : tuple of ints, optional
    Axes in `a` to reorder to the right, before inversion.
    If None (default), no reordering is done.

Returns
-------
x : ndarray, shape Q

Raises
------
LinAlgError
    If `a` is singular or not 'square' (in the above sense).

See Also
--------
tensordot, tensorinv, einsum

Examples
--------
>>> a = np.eye(2*3*4)
>>> a.shape = (2*3, 4, 2, 3, 4)
>>> b = np.random.randn(2*3, 4)
>>> x = np.linalg.tensorsolve(a, b)
>>> x.shape
(2, 3, 4)
>>> np.allclose(np.tensordot(a, x, axes=3), b)
True