def solve(self, b, u, axis, lu):
        """Solve Au=b

        Solve along axis if b and u are multidimensional arrays.

        Parameters
        ----------
        b, u : arrays of rhs and output
            Both can be multidimensional
        axis : int
            The axis we are solving over
        lu : LU-decomposition
            Can be either the output from splu, or a dia-matrix containing
            the L and U matrices. The latter is used in subclasses.
        """
        if axis > 0:
            u = np.moveaxis(u, axis, 0)
            if u is not b:
                b = np.moveaxis(b, axis, 0)

        s = slice(0, self.mat.shape[0])
        if b.ndim == 1:
            if b.dtype.char in 'fdg' or self.dtype in 'FDG':
                u[s] = lu.solve(b[s])

            else:
                u.real[s] = lu.solve(b[s].real)
                u.imag[s] = lu.solve(b[s].imag)

        else:
            N = b[s].shape[0]
            P = np.prod(b[s].shape[1:])
            br = b[s].reshape((N, P))
            if b.dtype.char in 'fdg' or self.dtype in 'FDG':
                u[s] = lu.solve(br).reshape(u[s].shape)
            else:
                u.real[s] = lu.solve(br.real).reshape(u[s].shape)
                u.imag[s] = lu.solve(br.imag).reshape(u[s].shape)

        if axis > 0:
            u = np.moveaxis(u, 0, axis)
            if u is not b:
                b = np.moveaxis(b, 0, axis)

        return u