'Solve the linear system Ax = rhs for x where A = self.
This is for symbolic matrices, for real or complex ones use
mpmath.lu_solve or mpmath.qr_solve.
See Also
lower_triangular_solve
upper_triangular_solve
gauss_jordan_solve
cholesky_solve
diagonal_solve
LDLsolve
QRsolve
pinv_solve
LUdecomposition'

def LUsolve(self, rhs, iszerofunc=_iszero):

    if (rhs.rows != self.rows):
        raise ShapeError('`self`    and    `rhs`    must    have    the    same    number    of    rows.')
    (A, perm) = self.LUdecomposition_Simple(iszerofunc=_iszero)
    n = self.rows
    b = rhs.permute_rows(perm).as_mutable()
    for i in range(n):
        for j in range(i):
            scale = A[(i, j)]
            b.zip_row_op(i, j, (lambda x, y: (x - (y * scale))))
    for i in range((n - 1), (-1), (-1)):
        for j in range((i + 1), n):
            scale = A[(i, j)]
            b.zip_row_op(i, j, (lambda x, y: (x - (y * scale))))
        scale = A[(i, i)]
        b.row_op(i, (lambda x, _: (x / scale)))
    return rhs.__class__(b)
