def matrix_representation(op):
    """Return a matrix representation of a linear operator.

    Parameters
    ----------
    op : `Operator`
        The linear operator of which one wants a matrix representation.
        If the domain or range is a `ProductSpace`, it must be a power-space.

    Returns
    -------
    matrix : `numpy.ndarray`
        The matrix representation of the operator.
        The shape will be ``op.domain.shape + op.range.shape`` and the dtype
        is the promoted (greatest) dtype of the domain and range.

    Examples
    --------
    Approximate a matrix on its own:

    >>> mat = np.array([[1, 2, 3],
    ...                 [4, 5, 6],
    ...                 [7, 8, 9]])
    >>> op = odl.MatrixOperator(mat)
    >>> matrix_representation(op)
    array([[1, 2, 3],
           [4, 5, 6],
           [7, 8, 9]])

    It also works with `ProductSpace`'s and higher dimensional `TensorSpace`'s.
    In this case, the returned "matrix" will also be higher dimensional:

    >>> space = odl.uniform_discr([0, 0], [2, 2], (2, 2))
    >>> grad = odl.Gradient(space)
    >>> tensor = odl.matrix_representation(grad)
    >>> tensor.shape == (2, 2, 2, 2, 2)
    True

    Since the "matrix" is now higher dimensional, we need to use e.g.
    `numpy.tensordot` if we want to compute with the matrix representation:

    >>> x = space.element(lambda x: x[0] ** 2 + 2 * x[1] ** 2)
    >>> grad(x)
    ProductSpace(uniform_discr([ 0.,  0.], [ 2.,  2.], (2, 2)), 2).element([
    <BLANKLINE>
            [[ 2.  ,  2.  ],
             [-2.75, -6.75]],
    <BLANKLINE>
            [[ 4.  , -4.75],
             [ 4.  , -6.75]]
    ])
    >>> np.tensordot(tensor, x, axes=grad.domain.ndim)
    array([[[ 2.  ,  2.  ],
            [-2.75, -6.75]],
    <BLANKLINE>
           [[ 4.  , -4.75],
            [ 4.  , -6.75]]])

    Notes
    ----------
    The algorithm works by letting the operator act on all unit vectors, and
    stacking the output as a matrix.
    """

    if not op.is_linear:
        raise ValueError('the operator is not linear')

    if not (isinstance(op.domain, TensorSpace) or
            (isinstance(op.domain, ProductSpace) and
             op.domain.is_power_space and
             all(isinstance(spc, TensorSpace) for spc in op.domain))):
        raise TypeError('operator domain {!r} is neither `TensorSpace` '
                        'nor `ProductSpace` with only equal `TensorSpace` '
                        'components'.format(op.domain))

    if not (isinstance(op.range, TensorSpace) or
            (isinstance(op.range, ProductSpace) and
             op.range.is_power_space and
             all(isinstance(spc, TensorSpace) for spc in op.range))):
        raise TypeError('operator range {!r} is neither `TensorSpace` '
                        'nor `ProductSpace` with only equal `TensorSpace` '
                        'components'.format(op.range))

    # Generate the matrix
    dtype = np.promote_types(op.domain.dtype, op.range.dtype)
    matrix = np.zeros(op.range.shape + op.domain.shape, dtype=dtype)
    tmp_ran = op.range.element()  # Store for reuse in loop
    tmp_dom = op.domain.zero()  # Store for reuse in loop

    for j in nd_iterator(op.domain.shape):
        tmp_dom[j] = 1.0

        op(tmp_dom, out=tmp_ran)
        matrix[(Ellipsis,) + j] = tmp_ran.asarray()

        tmp_dom[j] = 0.0

    return matrix
Return a matrix representation of a linear operator.

    Parameters
    ----------
    op : `Operator`
        The linear operator of which one wants a matrix representation.
        If the domain or range is a `ProductSpace`, it must be a power-space.

    Returns
    -------
    matrix : `numpy.ndarray`
        The matrix representation of the operator.
        The shape will be ``op.domain.shape + op.range.shape`` and the dtype
        is the promoted (greatest) dtype of the domain and range.

    Examples
    --------
    Approximate a matrix on its own:

    >>> mat = np.array([[1, 2, 3],
    ...                 [4, 5, 6],
    ...                 [7, 8, 9]])
    >>> op = odl.MatrixOperator(mat)
    >>> matrix_representation(op)
    array([[1, 2, 3],
           [4, 5, 6],
           [7, 8, 9]])

    It also works with `ProductSpace`'s and higher dimensional `TensorSpace`'s.
    In this case, the returned "matrix" will also be higher dimensional:

    >>> space = odl.uniform_discr([0, 0], [2, 2], (2, 2))
    >>> grad = odl.Gradient(space)
    >>> tensor = odl.matrix_representation(grad)
    >>> tensor.shape == (2, 2, 2, 2, 2)
    True

    Since the "matrix" is now higher dimensional, we need to use e.g.
    `numpy.tensordot` if we want to compute with the matrix representation:

    >>> x = space.element(lambda x: x[0] ** 2 + 2 * x[1] ** 2)
    >>> grad(x)
    ProductSpace(uniform_discr([ 0.,  0.], [ 2.,  2.], (2, 2)), 2).element([
    <BLANKLINE>
            [[ 2.  ,  2.  ],
             [-2.75, -6.75]],
    <BLANKLINE>
            [[ 4.  , -4.75],
             [ 4.  , -6.75]]
    ])
    >>> np.tensordot(tensor, x, axes=grad.domain.ndim)
    array([[[ 2.  ,  2.  ],
            [-2.75, -6.75]],
    <BLANKLINE>
           [[ 4.  , -4.75],
            [ 4.  , -6.75]]])

    Notes
    ----------
    The algorithm works by letting the operator act on all unit vectors, and
    stacking the output as a matrix.