Repository: sympy/sympy

Issue:
Series return an incorrect result
The latest sympy master:

```
In [1]: f=exp(x**3)*cos(x**6)

In [2]: f.series(x, 0, 14)
Out[2]: 
          6    9    12         
     3   x    x    x      ⎛ 14⎞
1 + x  + ── + ── + ─── + O⎝x  ⎠
         2    6     24         

In [3]: f.series(x, 0, 19)
Out[3]: 
          6    9       12       15        18         
     3   x    x    11⋅x     59⋅x     179⋅x      ⎛ 19⎞
1 + x  + ── + ── - ────── - ────── - ─────── + O⎝x  ⎠
         2    6      24      120       720           
```


Discussion:

Patch:
diff --git a/sympy/core/add.py b/sympy/core/add.py
index eca10ae1c6..4d22b7555a 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -878,7 +878,7 @@ def as_real_imag(self, deep=True, **hints):
         return (self.func(*re_part), self.func(*im_part))
 
     def _eval_as_leading_term(self, x):
-        from sympy import expand_mul, Order
+        from sympy import expand_mul, factor_terms
 
         old = self
 
@@ -888,36 +888,26 @@ def _eval_as_leading_term(self, x):
 
         infinite = [t for t in expr.args if t.is_infinite]
 
-        leading_terms = [t.as_leading_term(x) for t in expr.args]
-
-        min, new_expr = Order(0), 0
-
-        try:
-            for term in leading_terms:
-                order = Order(term, x)
-                if not min or order not in min:
-                    min = order
-                    new_expr = term
-                elif order == min:
-                    new_expr += term
-
-        except TypeError:
-            return expr
-
-        new_expr=new_expr.together()
-        if new_expr.is_Add:
-            new_expr = new_expr.simplify()
-
-        if not new_expr:
-            # simple leading term analysis gave us cancelled terms but we have to send
+        expr = expr.func(*[t.as_leading_term(x) for t in expr.args]).removeO()
+        if not expr:
+            # simple leading term analysis gave us 0 but we have to send
             # back a term, so compute the leading term (via series)
             return old.compute_leading_term(x)
-
-        elif new_expr is S.NaN:
+        elif expr is S.NaN:
             return old.func._from_args(infinite)
-
+        elif not expr.is_Add:
+            return expr
         else:
-            return new_expr
+            plain = expr.func(*[s for s, _ in expr.extract_leading_order(x)])
+            rv = factor_terms(plain, fraction=False)
+            rv_simplify = rv.simplify()
+            # if it simplifies to an x-free expression, return that;
+            # tests don't fail if we don't but it seems nicer to do this
+            if x not in rv_simplify.free_symbols:
+                if rv_simplify.is_zero and plain.is_zero is not True:
+                    return (expr - plain)._eval_as_leading_term(x)
+                return rv_simplify
+            return rv
 
     def _eval_adjoint(self):
         return self.func(*[t.adjoint() for t in self.args])
diff --git a/sympy/core/function.py b/sympy/core/function.py
index 9ab9f5c6aa..af1207faf1 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -692,9 +692,6 @@ def _eval_nseries(self, x, n, logx):
             if e == e1:
                 #for example when e = sin(x+1) or e = sin(cos(x))
                 #let's try the general algorithm
-                if len(e.args) == 1:
-                    # issue 14411
-                    e = e.func(e.args[0].cancel())
                 term = e.subs(x, S.Zero)
                 if term.is_finite is False or term is S.NaN:
                     raise PoleError("Cannot expand %s around 0" % (self))
@@ -724,7 +721,7 @@ def _eval_nseries(self, x, n, logx):
         nterms = n + 2
         cf = Order(arg.as_leading_term(x), x).getn()
         if cf != 0:
-            nterms = int(nterms / cf)
+            nterms = (n/cf).ceiling()
         for i in range(nterms):
             g = self.taylor_term(i, arg, g)
             g = g.nseries(x, n=n, logx=logx)
diff --git a/sympy/core/power.py b/sympy/core/power.py
index a3de65b720..b2e7c76b5d 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1456,13 +1456,11 @@ def _eval_nseries(self, x, n, logx):
 
                 # express "rest" as: rest = 1 + k*x**l + ... + O(x**n)
                 rest = expand_mul((b - prefactor)/prefactor)
-                rest = rest.simplify() #test_issue_6364
 
                 if rest.is_Order:
                     return 1/prefactor + rest/prefactor + O(x**n, x)
 
                 k, l = rest.leadterm(x)
-
                 if l.is_Rational and l > 0:
                     pass
                 elif l.is_number and l > 0:
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 9935305f7e..1821ff06e3 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -884,7 +884,7 @@ def _eval_is_extended_nonnegative(self):
     def _eval_nseries(self, x, n, logx):
         # NOTE Please see the comment at the beginning of this file, labelled
         #      IMPORTANT.
-        from sympy import cancel, Order, logcombine
+        from sympy import cancel, Order
         if not logx:
             logx = log(x)
         if self.args[0] == x:
@@ -905,8 +905,6 @@ def _eval_nseries(self, x, n, logx):
             s = self.args[0].nseries(x, n=n, logx=logx)
         a, b = s.leadterm(x)
         p = cancel(s/(a*x**b) - 1)
-        if p.has(exp):
-            p = logcombine(p)
         g = None
         l = []
         for i in range(n + 2):
diff --git a/sympy/matrices/decompositions.py b/sympy/matrices/decompositions.py
index 3ae82f0345..43726ba1ca 100644
--- a/sympy/matrices/decompositions.py
+++ b/sympy/matrices/decompositions.py
@@ -11,7 +11,7 @@
 from .determinant import _find_reasonable_pivot_naive
 
 
-def _rank_decomposition(M, iszerofunc=_iszero, simplify=False, dotprodsimp=None):
+def _rank_decomposition(M, iszerofunc=_iszero, simplify=False):
     r"""Returns a pair of matrices (`C`, `F`) with matching rank
     such that `A = C F`.
 
@@ -26,11 +26,6 @@ def _rank_decomposition(M, iszerofunc=_iszero, simplify=False, dotprodsimp=None)
         A function used to simplify elements when looking for a
         pivot. By default SymPy's ``simplify`` is used.
 
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
     Returns
     =======
 
@@ -102,14 +97,14 @@ def _rank_decomposition(M, iszerofunc=_iszero, simplify=False, dotprodsimp=None)
     rref
     """
 
-    F, pivot_cols = M.rref(simplify=simplify, iszerofunc=iszerofunc,
-            pivots=True, dotprodsimp=dotprodsimp)
+    (F, pivot_cols) = M.rref(
+        simplify=simplify, iszerofunc=iszerofunc, pivots=True)
     rank = len(pivot_cols)
 
     C = M.extract(range(M.rows), pivot_cols)
     F = F[:rank, :]
 
-    return C, F
+    return (C, F)
 
 
 def _liupc(M):
diff --git a/sympy/matrices/dense.py b/sympy/matrices/dense.py
index 90a0098768..c619c8c2a6 100644
--- a/sympy/matrices/dense.py
+++ b/sympy/matrices/dense.py
@@ -172,10 +172,69 @@ def _eval_matrix_mul_elementwise(self, other):
         return classof(self, other)._new(self.rows, self.cols, mat, copy=False)
 
     def _eval_inverse(self, **kwargs):
-        return self.inv(method=kwargs.get('method', 'GE'),
-                        iszerofunc=kwargs.get('iszerofunc', _iszero),
-                        dotprodsimp=kwargs.get('dotprodsimp', None),
-                        try_block_diag=kwargs.get('try_block_diag', False))
+        """Return the matrix inverse using the method indicated (default
+        is Gauss elimination).
+
+        kwargs
+        ======
+
+        method : ('GE', 'LU', or 'ADJ')
+        iszerofunc
+        try_block_diag
+
+        Notes
+        =====
+
+        According to the ``method`` keyword, it calls the appropriate method:
+
+          GE .... inverse_GE(); default
+          LU .... inverse_LU()
+          ADJ ... inverse_ADJ()
+
+        According to the ``try_block_diag`` keyword, it will try to form block
+        diagonal matrices using the method get_diag_blocks(), invert these
+        individually, and then reconstruct the full inverse matrix.
+
+        Note, the GE and LU methods may require the matrix to be simplified
+        before it is inverted in order to properly detect zeros during
+        pivoting. In difficult cases a custom zero detection function can
+        be provided by setting the ``iszerosfunc`` argument to a function that
+        should return True if its argument is zero. The ADJ routine computes
+        the determinant and uses that to detect singular matrices in addition
+        to testing for zeros on the diagonal.
+
+        See Also
+        ========
+
+        inverse_LU
+        inverse_GE
+        inverse_ADJ
+        """
+        from sympy.matrices import diag
+
+        method = kwargs.get('method', 'GE')
+        iszerofunc = kwargs.get('iszerofunc', _iszero)
+        dotprodsimp = kwargs.get('dotprodsimp', None)
+        if kwargs.get('try_block_diag', False):
+            blocks = self.get_diag_blocks()
+            r = []
+            for block in blocks:
+                r.append(block.inv(method=method, iszerofunc=iszerofunc,
+                        dotprodsimp=dotprodsimp))
+            return diag(*r)
+
+        M = self.as_mutable()
+        if method == "GE":
+            rv = M.inverse_GE(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
+        elif method == "LU":
+            rv = M.inverse_LU(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
+        elif method == "ADJ":
+            rv = M.inverse_ADJ(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
+        else:
+            # make sure to add an invertibility check (as in inverse_LU)
+            # if a new method is added.
+            raise ValueError("Inversion method unrecognized")
+        return self._new(rv)
 
     def _eval_scalar_mul(self, other):
         mat = [other*a for a in self._mat]
diff --git a/sympy/matrices/determinant.py b/sympy/matrices/determinant.py
index 88086ca332..7849675dc3 100644
--- a/sympy/matrices/determinant.py
+++ b/sympy/matrices/determinant.py
@@ -1,877 +1,877 @@
-from __future__ import division, print_function
-
-from types import FunctionType
-
-from sympy.core.numbers import Float, Integer
-from sympy.core.singleton import S
-from sympy.core.symbol import _uniquely_named_symbol
-from sympy.polys import PurePoly, cancel
-from sympy.simplify.simplify import (simplify as _simplify,
-    dotprodsimp as _dotprodsimp)
-
-from .common import MatrixError, NonSquareMatrixError
-from .utilities import _iszero, _is_zero_after_expand_mul
-
-
-def _find_reasonable_pivot(col, iszerofunc=_iszero, simpfunc=_simplify):
-    """ Find the lowest index of an item in ``col`` that is
-    suitable for a pivot.  If ``col`` consists only of
-    Floats, the pivot with the largest norm is returned.
-    Otherwise, the first element where ``iszerofunc`` returns
-    False is used.  If ``iszerofunc`` doesn't return false,
-    items are simplified and retested until a suitable
-    pivot is found.
-
-    Returns a 4-tuple
-        (pivot_offset, pivot_val, assumed_nonzero, newly_determined)
-    where pivot_offset is the index of the pivot, pivot_val is
-    the (possibly simplified) value of the pivot, assumed_nonzero
-    is True if an assumption that the pivot was non-zero
-    was made without being proved, and newly_determined are
-    elements that were simplified during the process of pivot
-    finding."""
-
-    newly_determined = []
-    col = list(col)
-    # a column that contains a mix of floats and integers
-    # but at least one float is considered a numerical
-    # column, and so we do partial pivoting
-    if all(isinstance(x, (Float, Integer)) for x in col) and any(
-            isinstance(x, Float) for x in col):
-        col_abs = [abs(x) for x in col]
-        max_value = max(col_abs)
-        if iszerofunc(max_value):
-            # just because iszerofunc returned True, doesn't
-            # mean the value is numerically zero.  Make sure
-            # to replace all entries with numerical zeros
-            if max_value != 0:
-                newly_determined = [(i, 0) for i, x in enumerate(col) if x != 0]
-            return (None, None, False, newly_determined)
-        index = col_abs.index(max_value)
-        return (index, col[index], False, newly_determined)
-
-    # PASS 1 (iszerofunc directly)
-    possible_zeros = []
-    for i, x in enumerate(col):
-        is_zero = iszerofunc(x)
-        # is someone wrote a custom iszerofunc, it may return
-        # BooleanFalse or BooleanTrue instead of True or False,
-        # so use == for comparison instead of `is`
-        if is_zero == False:
-            # we found something that is definitely not zero
-            return (i, x, False, newly_determined)
-        possible_zeros.append(is_zero)
-
-    # by this point, we've found no certain non-zeros
-    if all(possible_zeros):
-        # if everything is definitely zero, we have
-        # no pivot
-        return (None, None, False, newly_determined)
-
-    # PASS 2 (iszerofunc after simplify)
-    # we haven't found any for-sure non-zeros, so
-    # go through the elements iszerofunc couldn't
-    # make a determination about and opportunistically
-    # simplify to see if we find something
-    for i, x in enumerate(col):
-        if possible_zeros[i] is not None:
-            continue
-        simped = simpfunc(x)
-        is_zero = iszerofunc(simped)
-        if is_zero == True or is_zero == False:
-            newly_determined.append((i, simped))
-        if is_zero == False:
-            return (i, simped, False, newly_determined)
-        possible_zeros[i] = is_zero
-
-    # after simplifying, some things that were recognized
-    # as zeros might be zeros
-    if all(possible_zeros):
-        # if everything is definitely zero, we have
-        # no pivot
-        return (None, None, False, newly_determined)
-
-    # PASS 3 (.equals(0))
-    # some expressions fail to simplify to zero, but
-    # ``.equals(0)`` evaluates to True.  As a last-ditch
-    # attempt, apply ``.equals`` to these expressions
-    for i, x in enumerate(col):
-        if possible_zeros[i] is not None:
-            continue
-        if x.equals(S.Zero):
-            # ``.iszero`` may return False with
-            # an implicit assumption (e.g., ``x.equals(0)``
-            # when ``x`` is a symbol), so only treat it
-            # as proved when ``.equals(0)`` returns True
-            possible_zeros[i] = True
-            newly_determined.append((i, S.Zero))
-
-    if all(possible_zeros):
-        return (None, None, False, newly_determined)
-
-    # at this point there is nothing that could definitely
-    # be a pivot.  To maintain compatibility with existing
-    # behavior, we'll assume that an illdetermined thing is
-    # non-zero.  We should probably raise a warning in this case
-    i = possible_zeros.index(None)
-    return (i, col[i], True, newly_determined)
-
-
-def _find_reasonable_pivot_naive(col, iszerofunc=_iszero, simpfunc=None):
-    """
-    Helper that computes the pivot value and location from a
-    sequence of contiguous matrix column elements. As a side effect
-    of the pivot search, this function may simplify some of the elements
-    of the input column. A list of these simplified entries and their
-    indices are also returned.
-    This function mimics the behavior of _find_reasonable_pivot(),
-    but does less work trying to determine if an indeterminate candidate
-    pivot simplifies to zero. This more naive approach can be much faster,
-    with the trade-off that it may erroneously return a pivot that is zero.
-
-    ``col`` is a sequence of contiguous column entries to be searched for
-    a suitable pivot.
-    ``iszerofunc`` is a callable that returns a Boolean that indicates
-    if its input is zero, or None if no such determination can be made.
-    ``simpfunc`` is a callable that simplifies its input. It must return
-    its input if it does not simplify its input. Passing in
-    ``simpfunc=None`` indicates that the pivot search should not attempt
-    to simplify any candidate pivots.
-
-    Returns a 4-tuple:
-    (pivot_offset, pivot_val, assumed_nonzero, newly_determined)
-    ``pivot_offset`` is the sequence index of the pivot.
-    ``pivot_val`` is the value of the pivot.
-    pivot_val and col[pivot_index] are equivalent, but will be different
-    when col[pivot_index] was simplified during the pivot search.
-    ``assumed_nonzero`` is a boolean indicating if the pivot cannot be
-    guaranteed to be zero. If assumed_nonzero is true, then the pivot
-    may or may not be non-zero. If assumed_nonzero is false, then
-    the pivot is non-zero.
-    ``newly_determined`` is a list of index-value pairs of pivot candidates
-    that were simplified during the pivot search.
-    """
-
-    # indeterminates holds the index-value pairs of each pivot candidate
-    # that is neither zero or non-zero, as determined by iszerofunc().
-    # If iszerofunc() indicates that a candidate pivot is guaranteed
-    # non-zero, or that every candidate pivot is zero then the contents
-    # of indeterminates are unused.
-    # Otherwise, the only viable candidate pivots are symbolic.
-    # In this case, indeterminates will have at least one entry,
-    # and all but the first entry are ignored when simpfunc is None.
-    indeterminates = []
-    for i, col_val in enumerate(col):
-        col_val_is_zero = iszerofunc(col_val)
-        if col_val_is_zero == False:
-            # This pivot candidate is non-zero.
-            return i, col_val, False, []
-        elif col_val_is_zero is None:
-            # The candidate pivot's comparison with zero
-            # is indeterminate.
-            indeterminates.append((i, col_val))
-
-    if len(indeterminates) == 0:
-        # All candidate pivots are guaranteed to be zero, i.e. there is
-        # no pivot.
-        return None, None, False, []
-
-    if simpfunc is None:
-        # Caller did not pass in a simplification function that might
-        # determine if an indeterminate pivot candidate is guaranteed
-        # to be nonzero, so assume the first indeterminate candidate
-        # is non-zero.
-        return indeterminates[0][0], indeterminates[0][1], True, []
-
-    # newly_determined holds index-value pairs of candidate pivots
-    # that were simplified during the search for a non-zero pivot.
-    newly_determined = []
-    for i, col_val in indeterminates:
-        tmp_col_val = simpfunc(col_val)
-        if id(col_val) != id(tmp_col_val):
-            # simpfunc() simplified this candidate pivot.
-            newly_determined.append((i, tmp_col_val))
-            if iszerofunc(tmp_col_val) == False:
-                # Candidate pivot simplified to a guaranteed non-zero value.
-                return i, tmp_col_val, False, newly_determined
-
-    return indeterminates[0][0], indeterminates[0][1], True, newly_determined
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _berkowitz_toeplitz_matrix(M, dotprodsimp=None):
-    """Return (A,T) where T the Toeplitz matrix used in the Berkowitz algorithm
-    corresponding to ``M`` and A is the first principal submatrix.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-    """
-
-    # the 0 x 0 case is trivial
-    if M.rows == 0 and M.cols == 0:
-        return M._new(1,1, [M.one])
-
-    #
-    # Partition M = [ a_11  R ]
-    #                  [ C     A ]
-    #
-
-    a, R = M[0,0],   M[0, 1:]
-    C, A = M[1:, 0], M[1:,1:]
-
-    #
-    # The Toeplitz matrix looks like
-    #
-    #  [ 1                                     ]
-    #  [ -a         1                          ]
-    #  [ -RC       -a        1                 ]
-    #  [ -RAC     -RC       -a       1         ]
-    #  [ -RA**2C -RAC      -RC      -a       1 ]
-    #  etc.
-
-    # Compute the diagonal entries.
-    # Because multiplying matrix times vector is so much
-    # more efficient than matrix times matrix, recursively
-    # compute -R * A**n * C.
-    diags = [C]
-    for i in range(M.rows - 2):
-        diags.append(A.multiply(diags[i], dotprodsimp=dotprodsimp))
-    diags = [(-R).multiply(d, dotprodsimp=dotprodsimp)[0, 0] for d in diags]
-    diags = [M.one, -a] + diags
-
-    def entry(i,j):
-        if j > i:
-            return M.zero
-        return diags[i - j]
-
-    toeplitz = M._new(M.cols + 1, M.rows, entry)
-    return (A, toeplitz)
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _berkowitz_vector(M, dotprodsimp=None):
-    """ Run the Berkowitz algorithm and return a vector whose entries
-        are the coefficients of the characteristic polynomial of ``M``.
-
-        Given N x N matrix, efficiently compute
-        coefficients of characteristic polynomials of ``M``
-        without division in the ground domain.
-
-        This method is particularly useful for computing determinant,
-        principal minors and characteristic polynomial when ``M``
-        has complicated coefficients e.g. polynomials. Semi-direct
-        usage of this algorithm is also important in computing
-        efficiently sub-resultant PRS.
-
-        Assuming that M is a square matrix of dimension N x N and
-        I is N x N identity matrix, then the Berkowitz vector is
-        an N x 1 vector whose entries are coefficients of the
-        polynomial
-
-                        charpoly(M) = det(t*I - M)
-
-        As a consequence, all polynomials generated by Berkowitz
-        algorithm are monic.
-
-        Parameters
-        ==========
-
-        dotprodsimp : bool, optional
-            Specifies whether intermediate term algebraic simplification is
-            used during matrix multiplications to control expression blowup
-            and thus speed up calculation.
-
-        For more information on the implemented algorithm refer to:
-
-        [1] S.J. Berkowitz, On computing the determinant in small
-            parallel time using a small number of processors, ACM,
-            Information Processing Letters 18, 1984, pp. 147-150
-
-        [2] M. Keber, Division-Free computation of sub-resultants
-            using Bezout matrices, Tech. Report MPI-I-2006-1-006,
-            Saarbrucken, 2006
-    """
-
-    # handle the trivial cases
-    if M.rows == 0 and M.cols == 0:
-        return M._new(1, 1, [M.one])
-    elif M.rows == 1 and M.cols == 1:
-        return M._new(2, 1, [M.one, -M[0,0]])
-
-    submat, toeplitz = _berkowitz_toeplitz_matrix(M, dotprodsimp=dotprodsimp)
-
-    return toeplitz.multiply(_berkowitz_vector(submat, dotprodsimp=dotprodsimp),
-            dotprodsimp=dotprodsimp)
-
-
-def _adjugate(M, method="berkowitz", dotprodsimp=None):
-    """Returns the adjugate, or classical adjoint, of
-    a matrix.  That is, the transpose of the matrix of cofactors.
-
-    https://en.wikipedia.org/wiki/Adjugate
-
-    Parameters
-    ==========
-
-    method : string, optional
-        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
-        "lu".
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        to control expression blowup during matrix multiplication. If this
-        is true then the simplify function is not used.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2], [3, 4]])
-    >>> M.adjugate()
-    Matrix([
-    [ 4, -2],
-    [-3,  1]])
-
-    See Also
-    ========
-
-    cofactor_matrix
-    sympy.matrices.common.MatrixCommon.transpose
-    """
-
-    return M.cofactor_matrix(method=method, dotprodsimp=dotprodsimp).transpose()
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _charpoly(M, x='lambda', simplify=_simplify, dotprodsimp=None):
-    """Computes characteristic polynomial det(x*I - M) where I is
-    the identity matrix.
-
-    A PurePoly is returned, so using different variables for ``x`` does
-    not affect the comparison or the polynomials:
-
-    Parameters
-    ==========
-
-    x : string, optional
-        Name for the "lambda" variable, defaults to "lambda".
-
-    simplify : function, optional
-        Simplification function to use on the characteristic polynomial
-        calculated. Defaults to ``simplify``, if ``dotprodsimp`` is ``True``
-        then this is ignored.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        to control expression blowup during matrix multiplication. If this
-        is true then the simplify function is not used.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> from sympy.abc import x, y
-    >>> M = Matrix([[1, 3], [2, 0]])
-    >>> M.charpoly()
-    PurePoly(lambda**2 - lambda - 6, lambda, domain='ZZ')
-    >>> M.charpoly(x) == M.charpoly(y)
-    True
-    >>> M.charpoly(x) == M.charpoly(y)
-    True
-
-    Specifying ``x`` is optional; a symbol named ``lambda`` is used by
-    default (which looks good when pretty-printed in unicode):
-
-    >>> M.charpoly().as_expr()
-    lambda**2 - lambda - 6
-
-    And if ``x`` clashes with an existing symbol, underscores will
-    be prepended to the name to make it unique:
-
-    >>> M = Matrix([[1, 2], [x, 0]])
-    >>> M.charpoly(x).as_expr()
-    _x**2 - _x - 2*x
-
-    Whether you pass a symbol or not, the generator can be obtained
-    with the gen attribute since it may not be the same as the symbol
-    that was passed:
-
-    >>> M.charpoly(x).gen
-    _x
-    >>> M.charpoly(x).gen == x
-    False
-
-    Notes
-    =====
-
-    The Samuelson-Berkowitz algorithm is used to compute
-    the characteristic polynomial efficiently and without any
-    division operations.  Thus the characteristic polynomial over any
-    commutative ring without zero divisors can be computed.
-
-    See Also
-    ========
-
-    det
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    if dotprodsimp:
-        simplify = lambda e: e
-
-    berk_vector = _berkowitz_vector(M, dotprodsimp=dotprodsimp)
-    x = _uniquely_named_symbol(x, berk_vector)
-
-    return PurePoly([simplify(a) for a in berk_vector], x)
-
-
-def _cofactor(M, i, j, method="berkowitz", dotprodsimp=None):
-    """Calculate the cofactor of an element.
-
-    Parameters
-    ==========
-
-    method : string, optional
-        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
-        "lu".
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        to control expression blowup during matrix multiplication. If this
-        is true then the simplify function is not used.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2], [3, 4]])
-    >>> M.cofactor(0, 1)
-    -3
-
-    See Also
-    ========
-
-    cofactor_matrix
-    minor
-    minor_submatrix
-    """
-
-    if not M.is_square or M.rows < 1:
-        raise NonSquareMatrixError()
-
-    return (-1)**((i + j) % 2) * M.minor(i, j, method, dotprodsimp=dotprodsimp)
-
-
-def _cofactor_matrix(M, method="berkowitz", dotprodsimp=None):
-    """Return a matrix containing the cofactor of each element.
-
-    Parameters
-    ==========
-
-    method : string, optional
-        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
-        "lu".
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        to control expression blowup during matrix multiplication. If this
-        is true then the simplify function is not used.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2], [3, 4]])
-    >>> M.cofactor_matrix()
-    Matrix([
-    [ 4, -3],
-    [-2,  1]])
-
-    See Also
-    ========
-
-    cofactor
-    minor
-    minor_submatrix
-    """
-
-    if not M.is_square or M.rows < 1:
-        raise NonSquareMatrixError()
-
-    return M._new(M.rows, M.cols,
-            lambda i, j: M.cofactor(i, j, method, dotprodsimp=dotprodsimp))
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _det(M, method="bareiss", iszerofunc=None, dotprodsimp=None):
-    """Computes the determinant of a matrix if ``M`` is a concrete matrix object
-    otherwise return an expressions ``Determinant(M)`` if ``M`` is a
-    ``MatrixSymbol`` or other expression.
-
-    Parameters
-    ==========
-
-    method : string, optional
-        Specifies the algorithm used for computing the matrix determinant.
-
-        If the matrix is at most 3x3, a hard-coded formula is used and the
-        specified method is ignored. Otherwise, it defaults to
-        ``'bareiss'``.
-
-        If it is set to ``'bareiss'``, Bareiss' fraction-free algorithm will
-        be used.
-
-        If it is set to ``'berkowitz'``, Berkowitz' algorithm will be used.
-
-        Otherwise, if it is set to ``'lu'``, LU decomposition will be used.
-
-        .. note::
-            For backward compatibility, legacy keys like "bareis" and
-            "det_lu" can still be used to indicate the corresponding
-            methods.
-            And the keys are also case-insensitive for now. However, it is
-            suggested to use the precise keys for specifying the method.
-
-    iszerofunc : FunctionType or None, optional
-        If it is set to ``None``, it will be defaulted to ``_iszero`` if the
-        method is set to ``'bareiss'``, and ``_is_zero_after_expand_mul`` if
-        the method is set to ``'lu'``.
-
-        It can also accept any user-specified zero testing function, if it
-        is formatted as a function which accepts a single symbolic argument
-        and returns ``True`` if it is tested as zero and ``False`` if it
-        tested as non-zero, and also ``None`` if it is undecidable.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Returns
-    =======
-
-    det : Basic
-        Result of determinant.
-
-    Raises
-    ======
-
-    ValueError
-        If unrecognized keys are given for ``method`` or ``iszerofunc``.
-
-    NonSquareMatrixError
-        If attempted to calculate determinant from a non-square matrix.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix, MatrixSymbol, eye, det
-    >>> M = Matrix([[1, 2], [3, 4]])
-    >>> M.det()
-    -2
-    """
-
-    # sanitize `method`
-    method = method.lower()
-
-    if method == "bareis":
-        method = "bareiss"
-    elif method == "det_lu":
-        method = "lu"
-
-    if method not in ("bareiss", "berkowitz", "lu"):
-        raise ValueError("Determinant method '%s' unrecognized" % method)
-
-    if iszerofunc is None:
-        if method == "bareiss":
-            iszerofunc = _is_zero_after_expand_mul
-        elif method == "lu":
-            iszerofunc = _iszero
-
-    elif not isinstance(iszerofunc, FunctionType):
-        raise ValueError("Zero testing method '%s' unrecognized" % iszerofunc)
-
-    n = M.rows
-
-    if n == M.cols: # square check is done in individual method functions
-        if n == 0:
-            return M.one
-        elif n == 1:
-            return M[0,0]
-        elif n == 2:
-            m = M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]
-            return _dotprodsimp(m) if dotprodsimp else m
-        elif n == 3:
-            m =  (M[0, 0] * M[1, 1] * M[2, 2]
-                + M[0, 1] * M[1, 2] * M[2, 0]
-                + M[0, 2] * M[1, 0] * M[2, 1]
-                - M[0, 2] * M[1, 1] * M[2, 0]
-                - M[0, 0] * M[1, 2] * M[2, 1]
-                - M[0, 1] * M[1, 0] * M[2, 2])
-            return _dotprodsimp(m) if dotprodsimp else m
-
-    if method == "bareiss":
-        return M._eval_det_bareiss(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    elif method == "berkowitz":
-        return M._eval_det_berkowitz(dotprodsimp=dotprodsimp)
-    elif method == "lu":
-        return M._eval_det_lu(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    else:
-        raise MatrixError('unknown method for calculating determinant')
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _det_bareiss(M, iszerofunc=_is_zero_after_expand_mul, dotprodsimp=None):
-    """Compute matrix determinant using Bareiss' fraction-free
-    algorithm which is an extension of the well known Gaussian
-    elimination method. This approach is best suited for dense
-    symbolic matrices and will result in a determinant with
-    minimal number of fractions. It means that less term
-    rewriting is needed on resulting formulae.
-
-    Parameters
-    ==========
-
-    iszerofunc : function, optional
-        The function to use to determine zeros when doing an LU decomposition.
-        Defaults to ``lambda x: x.is_zero``.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    TODO: Implement algorithm for sparse matrices (SFF),
-    http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps.
-    """
-
-    # Recursively implemented Bareiss' algorithm as per Deanna Richelle Leggett's
-    # thesis http://www.math.usm.edu/perry/Research/Thesis_DRL.pdf
-    def bareiss(mat, cumm=1):
-        if mat.rows == 0:
-            return mat.one
-        elif mat.rows == 1:
-            return mat[0, 0]
-
-        # find a pivot and extract the remaining matrix
-        # With the default iszerofunc, _find_reasonable_pivot slows down
-        # the computation by the factor of 2.5 in one test.
-        # Relevant issues: #10279 and #13877.
-        pivot_pos, pivot_val, _, _ = _find_reasonable_pivot(mat[:, 0], iszerofunc=iszerofunc)
-        if pivot_pos is None:
-            return mat.zero
-
-        # if we have a valid pivot, we'll do a "row swap", so keep the
-        # sign of the det
-        sign = (-1) ** (pivot_pos % 2)
-
-        # we want every row but the pivot row and every column
-        rows = list(i for i in range(mat.rows) if i != pivot_pos)
-        cols = list(range(mat.cols))
-        tmp_mat = mat.extract(rows, cols)
-
-        def entry(i, j):
-            ret = (pivot_val*tmp_mat[i, j + 1] - mat[pivot_pos, j + 1]*tmp_mat[i, 0]) / cumm
-            if dotprodsimp:
-                return _dotprodsimp(ret)
-            elif not ret.is_Atom:
-                return cancel(ret)
-            return ret
-
-        return sign*bareiss(M._new(mat.rows - 1, mat.cols - 1, entry), pivot_val)
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    if M.rows == 0:
-        return M.one
-        # sympy/matrices/tests/test_matrices.py contains a test that
-        # suggests that the determinant of a 0 x 0 matrix is one, by
-        # convention.
-
-    return bareiss(M)
-
-
-def _det_berkowitz(M, dotprodsimp=None):
-    """ Use the Berkowitz algorithm to compute the determinant.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    if M.rows == 0:
-        return M.one
-        # sympy/matrices/tests/test_matrices.py contains a test that
-        # suggests that the determinant of a 0 x 0 matrix is one, by
-        # convention.
-
-    berk_vector = _berkowitz_vector(M, dotprodsimp=dotprodsimp)
-    return (-1)**(len(berk_vector) - 1) * berk_vector[-1]
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _det_LU(M, iszerofunc=_iszero, simpfunc=None, dotprodsimp=None):
-    """ Computes the determinant of a matrix from its LU decomposition.
-    This function uses the LU decomposition computed by
-    LUDecomposition_Simple().
-
-    The keyword arguments iszerofunc and simpfunc are passed to
-    LUDecomposition_Simple().
-    iszerofunc is a callable that returns a boolean indicating if its
-    input is zero, or None if it cannot make the determination.
-    simpfunc is a callable that simplifies its input.
-    The default is simpfunc=None, which indicate that the pivot search
-    algorithm should not attempt to simplify any candidate pivots.
-    If simpfunc fails to simplify its input, then it must return its input
-    instead of a copy.
-
-    Parameters
-    ==========
-
-    iszerofunc : function, optional
-        The function to use to determine zeros when doing an LU decomposition.
-        Defaults to ``lambda x: x.is_zero``.
-
-    simpfunc : function, optional
-        The simplification function to use when looking for zeros for pivots.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    if M.rows == 0:
-        return M.one
-        # sympy/matrices/tests/test_matrices.py contains a test that
-        # suggests that the determinant of a 0 x 0 matrix is one, by
-        # convention.
-
-    lu, row_swaps = M.LUdecomposition_Simple(iszerofunc=iszerofunc,
-            simpfunc=simpfunc, dotprodsimp=dotprodsimp)
-    # P*A = L*U => det(A) = det(L)*det(U)/det(P) = det(P)*det(U).
-    # Lower triangular factor L encoded in lu has unit diagonal => det(L) = 1.
-    # P is a permutation matrix => det(P) in {-1, 1} => 1/det(P) = det(P).
-    # LUdecomposition_Simple() returns a list of row exchange index pairs, rather
-    # than a permutation matrix, but det(P) = (-1)**len(row_swaps).
-
-    # Avoid forming the potentially time consuming  product of U's diagonal entries
-    # if the product is zero.
-    # Bottom right entry of U is 0 => det(A) = 0.
-    # It may be impossible to determine if this entry of U is zero when it is symbolic.
-    if iszerofunc(lu[lu.rows-1, lu.rows-1]):
-        return M.zero
-
-    # Compute det(P)
-    det = -M.one if len(row_swaps)%2 else M.one
-
-    # Compute det(U) by calculating the product of U's diagonal entries.
-    # The upper triangular portion of lu is the upper triangular portion of the
-    # U factor in the LU decomposition.
-    for k in range(lu.rows):
-        det *= lu[k, k]
-
-    # return det(P)*det(U)
-    return det
-
-
-def _minor(M, i, j, method="berkowitz", dotprodsimp=None):
-    """Return the (i,j) minor of ``M``.  That is,
-    return the determinant of the matrix obtained by deleting
-    the `i`th row and `j`th column from ``M``.
-
-    Parameters
-    ==========
-
-    i, j : int
-        The row and column to exclude to obtain the submatrix.
-
-    method : string, optional
-        Method to use to find the determinant of the submatrix, can be
-        "bareiss", "berkowitz" or "lu".
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
-    >>> M.minor(1, 1)
-    -12
-
-    See Also
-    ========
-
-    minor_submatrix
-    cofactor
-    det
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    return M.minor_submatrix(i, j).det(method=method, dotprodsimp=dotprodsimp)
-
-
-def _minor_submatrix(M, i, j):
-    """Return the submatrix obtained by removing the `i`th row
-    and `j`th column from ``M`` (works with Pythonic negative indices).
-
-    Parameters
-    ==========
-
-    i, j : int
-        The row and column to exclude to obtain the submatrix.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
-    >>> M.minor_submatrix(1, 1)
-    Matrix([
-    [1, 3],
-    [7, 9]])
-
-    See Also
-    ========
-
-    minor
-    cofactor
-    """
-
-    if i < 0:
-        i += M.rows
-    if j < 0:
-        j += M.cols
-
-    if not 0 <= i < M.rows or not 0 <= j < M.cols:
-        raise ValueError("`i` and `j` must satisfy 0 <= i < ``M.rows`` "
-                            "(%d)" % M.rows + "and 0 <= j < ``M.cols`` (%d)." % M.cols)
-
-    rows = [a for a in range(M.rows) if a != i]
-    cols = [a for a in range(M.cols) if a != j]
-
-    return M.extract(rows, cols)
+from __future__ import division, print_function
+
+from types import FunctionType
+
+from sympy.core.numbers import Float, Integer
+from sympy.core.singleton import S
+from sympy.core.symbol import _uniquely_named_symbol
+from sympy.polys import PurePoly, cancel
+from sympy.simplify.simplify import (simplify as _simplify,
+    dotprodsimp as _dotprodsimp)
+
+from .common import MatrixError, NonSquareMatrixError
+from .utilities import _iszero, _is_zero_after_expand_mul
+
+
+def _find_reasonable_pivot(col, iszerofunc=_iszero, simpfunc=_simplify):
+    """ Find the lowest index of an item in ``col`` that is
+    suitable for a pivot.  If ``col`` consists only of
+    Floats, the pivot with the largest norm is returned.
+    Otherwise, the first element where ``iszerofunc`` returns
+    False is used.  If ``iszerofunc`` doesn't return false,
+    items are simplified and retested until a suitable
+    pivot is found.
+
+    Returns a 4-tuple
+        (pivot_offset, pivot_val, assumed_nonzero, newly_determined)
+    where pivot_offset is the index of the pivot, pivot_val is
+    the (possibly simplified) value of the pivot, assumed_nonzero
+    is True if an assumption that the pivot was non-zero
+    was made without being proved, and newly_determined are
+    elements that were simplified during the process of pivot
+    finding."""
+
+    newly_determined = []
+    col = list(col)
+    # a column that contains a mix of floats and integers
+    # but at least one float is considered a numerical
+    # column, and so we do partial pivoting
+    if all(isinstance(x, (Float, Integer)) for x in col) and any(
+            isinstance(x, Float) for x in col):
+        col_abs = [abs(x) for x in col]
+        max_value = max(col_abs)
+        if iszerofunc(max_value):
+            # just because iszerofunc returned True, doesn't
+            # mean the value is numerically zero.  Make sure
+            # to replace all entries with numerical zeros
+            if max_value != 0:
+                newly_determined = [(i, 0) for i, x in enumerate(col) if x != 0]
+            return (None, None, False, newly_determined)
+        index = col_abs.index(max_value)
+        return (index, col[index], False, newly_determined)
+
+    # PASS 1 (iszerofunc directly)
+    possible_zeros = []
+    for i, x in enumerate(col):
+        is_zero = iszerofunc(x)
+        # is someone wrote a custom iszerofunc, it may return
+        # BooleanFalse or BooleanTrue instead of True or False,
+        # so use == for comparison instead of `is`
+        if is_zero == False:
+            # we found something that is definitely not zero
+            return (i, x, False, newly_determined)
+        possible_zeros.append(is_zero)
+
+    # by this point, we've found no certain non-zeros
+    if all(possible_zeros):
+        # if everything is definitely zero, we have
+        # no pivot
+        return (None, None, False, newly_determined)
+
+    # PASS 2 (iszerofunc after simplify)
+    # we haven't found any for-sure non-zeros, so
+    # go through the elements iszerofunc couldn't
+    # make a determination about and opportunistically
+    # simplify to see if we find something
+    for i, x in enumerate(col):
+        if possible_zeros[i] is not None:
+            continue
+        simped = simpfunc(x)
+        is_zero = iszerofunc(simped)
+        if is_zero == True or is_zero == False:
+            newly_determined.append((i, simped))
+        if is_zero == False:
+            return (i, simped, False, newly_determined)
+        possible_zeros[i] = is_zero
+
+    # after simplifying, some things that were recognized
+    # as zeros might be zeros
+    if all(possible_zeros):
+        # if everything is definitely zero, we have
+        # no pivot
+        return (None, None, False, newly_determined)
+
+    # PASS 3 (.equals(0))
+    # some expressions fail to simplify to zero, but
+    # ``.equals(0)`` evaluates to True.  As a last-ditch
+    # attempt, apply ``.equals`` to these expressions
+    for i, x in enumerate(col):
+        if possible_zeros[i] is not None:
+            continue
+        if x.equals(S.Zero):
+            # ``.iszero`` may return False with
+            # an implicit assumption (e.g., ``x.equals(0)``
+            # when ``x`` is a symbol), so only treat it
+            # as proved when ``.equals(0)`` returns True
+            possible_zeros[i] = True
+            newly_determined.append((i, S.Zero))
+
+    if all(possible_zeros):
+        return (None, None, False, newly_determined)
+
+    # at this point there is nothing that could definitely
+    # be a pivot.  To maintain compatibility with existing
+    # behavior, we'll assume that an illdetermined thing is
+    # non-zero.  We should probably raise a warning in this case
+    i = possible_zeros.index(None)
+    return (i, col[i], True, newly_determined)
+
+
+def _find_reasonable_pivot_naive(col, iszerofunc=_iszero, simpfunc=None):
+    """
+    Helper that computes the pivot value and location from a
+    sequence of contiguous matrix column elements. As a side effect
+    of the pivot search, this function may simplify some of the elements
+    of the input column. A list of these simplified entries and their
+    indices are also returned.
+    This function mimics the behavior of _find_reasonable_pivot(),
+    but does less work trying to determine if an indeterminate candidate
+    pivot simplifies to zero. This more naive approach can be much faster,
+    with the trade-off that it may erroneously return a pivot that is zero.
+
+    ``col`` is a sequence of contiguous column entries to be searched for
+    a suitable pivot.
+    ``iszerofunc`` is a callable that returns a Boolean that indicates
+    if its input is zero, or None if no such determination can be made.
+    ``simpfunc`` is a callable that simplifies its input. It must return
+    its input if it does not simplify its input. Passing in
+    ``simpfunc=None`` indicates that the pivot search should not attempt
+    to simplify any candidate pivots.
+
+    Returns a 4-tuple:
+    (pivot_offset, pivot_val, assumed_nonzero, newly_determined)
+    ``pivot_offset`` is the sequence index of the pivot.
+    ``pivot_val`` is the value of the pivot.
+    pivot_val and col[pivot_index] are equivalent, but will be different
+    when col[pivot_index] was simplified during the pivot search.
+    ``assumed_nonzero`` is a boolean indicating if the pivot cannot be
+    guaranteed to be zero. If assumed_nonzero is true, then the pivot
+    may or may not be non-zero. If assumed_nonzero is false, then
+    the pivot is non-zero.
+    ``newly_determined`` is a list of index-value pairs of pivot candidates
+    that were simplified during the pivot search.
+    """
+
+    # indeterminates holds the index-value pairs of each pivot candidate
+    # that is neither zero or non-zero, as determined by iszerofunc().
+    # If iszerofunc() indicates that a candidate pivot is guaranteed
+    # non-zero, or that every candidate pivot is zero then the contents
+    # of indeterminates are unused.
+    # Otherwise, the only viable candidate pivots are symbolic.
+    # In this case, indeterminates will have at least one entry,
+    # and all but the first entry are ignored when simpfunc is None.
+    indeterminates = []
+    for i, col_val in enumerate(col):
+        col_val_is_zero = iszerofunc(col_val)
+        if col_val_is_zero == False:
+            # This pivot candidate is non-zero.
+            return i, col_val, False, []
+        elif col_val_is_zero is None:
+            # The candidate pivot's comparison with zero
+            # is indeterminate.
+            indeterminates.append((i, col_val))
+
+    if len(indeterminates) == 0:
+        # All candidate pivots are guaranteed to be zero, i.e. there is
+        # no pivot.
+        return None, None, False, []
+
+    if simpfunc is None:
+        # Caller did not pass in a simplification function that might
+        # determine if an indeterminate pivot candidate is guaranteed
+        # to be nonzero, so assume the first indeterminate candidate
+        # is non-zero.
+        return indeterminates[0][0], indeterminates[0][1], True, []
+
+    # newly_determined holds index-value pairs of candidate pivots
+    # that were simplified during the search for a non-zero pivot.
+    newly_determined = []
+    for i, col_val in indeterminates:
+        tmp_col_val = simpfunc(col_val)
+        if id(col_val) != id(tmp_col_val):
+            # simpfunc() simplified this candidate pivot.
+            newly_determined.append((i, tmp_col_val))
+            if iszerofunc(tmp_col_val) == False:
+                # Candidate pivot simplified to a guaranteed non-zero value.
+                return i, tmp_col_val, False, newly_determined
+
+    return indeterminates[0][0], indeterminates[0][1], True, newly_determined
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _berkowitz_toeplitz_matrix(M, dotprodsimp=None):
+    """Return (A,T) where T the Toeplitz matrix used in the Berkowitz algorithm
+    corresponding to ``M`` and A is the first principal submatrix.
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+    """
+
+    # the 0 x 0 case is trivial
+    if M.rows == 0 and M.cols == 0:
+        return M._new(1,1, [M.one])
+
+    #
+    # Partition M = [ a_11  R ]
+    #                  [ C     A ]
+    #
+
+    a, R = M[0,0],   M[0, 1:]
+    C, A = M[1:, 0], M[1:,1:]
+
+    #
+    # The Toeplitz matrix looks like
+    #
+    #  [ 1                                     ]
+    #  [ -a         1                          ]
+    #  [ -RC       -a        1                 ]
+    #  [ -RAC     -RC       -a       1         ]
+    #  [ -RA**2C -RAC      -RC      -a       1 ]
+    #  etc.
+
+    # Compute the diagonal entries.
+    # Because multiplying matrix times vector is so much
+    # more efficient than matrix times matrix, recursively
+    # compute -R * A**n * C.
+    diags = [C]
+    for i in range(M.rows - 2):
+        diags.append(A.multiply(diags[i], dotprodsimp=dotprodsimp))
+    diags = [(-R).multiply(d, dotprodsimp=dotprodsimp)[0, 0] for d in diags]
+    diags = [M.one, -a] + diags
+
+    def entry(i,j):
+        if j > i:
+            return M.zero
+        return diags[i - j]
+
+    toeplitz = M._new(M.cols + 1, M.rows, entry)
+    return (A, toeplitz)
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _berkowitz_vector(M, dotprodsimp=None):
+    """ Run the Berkowitz algorithm and return a vector whose entries
+        are the coefficients of the characteristic polynomial of ``M``.
+
+        Given N x N matrix, efficiently compute
+        coefficients of characteristic polynomials of ``M``
+        without division in the ground domain.
+
+        This method is particularly useful for computing determinant,
+        principal minors and characteristic polynomial when ``M``
+        has complicated coefficients e.g. polynomials. Semi-direct
+        usage of this algorithm is also important in computing
+        efficiently sub-resultant PRS.
+
+        Assuming that M is a square matrix of dimension N x N and
+        I is N x N identity matrix, then the Berkowitz vector is
+        an N x 1 vector whose entries are coefficients of the
+        polynomial
+
+                        charpoly(M) = det(t*I - M)
+
+        As a consequence, all polynomials generated by Berkowitz
+        algorithm are monic.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is
+            used during matrix multiplications to control expression blowup
+            and thus speed up calculation.
+
+        For more information on the implemented algorithm refer to:
+
+        [1] S.J. Berkowitz, On computing the determinant in small
+            parallel time using a small number of processors, ACM,
+            Information Processing Letters 18, 1984, pp. 147-150
+
+        [2] M. Keber, Division-Free computation of sub-resultants
+            using Bezout matrices, Tech. Report MPI-I-2006-1-006,
+            Saarbrucken, 2006
+    """
+
+    # handle the trivial cases
+    if M.rows == 0 and M.cols == 0:
+        return M._new(1, 1, [M.one])
+    elif M.rows == 1 and M.cols == 1:
+        return M._new(2, 1, [M.one, -M[0,0]])
+
+    submat, toeplitz = _berkowitz_toeplitz_matrix(M, dotprodsimp=dotprodsimp)
+
+    return toeplitz.multiply(_berkowitz_vector(submat, dotprodsimp=dotprodsimp),
+            dotprodsimp=dotprodsimp)
+
+
+def _adjugate(M, method="berkowitz", dotprodsimp=None):
+    """Returns the adjugate, or classical adjoint, of
+    a matrix.  That is, the transpose of the matrix of cofactors.
+
+    https://en.wikipedia.org/wiki/Adjugate
+
+    Parameters
+    ==========
+
+    method : string, optional
+        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
+        "lu".
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        to control expression blowup during matrix multiplication. If this
+        is true then the simplify function is not used.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2], [3, 4]])
+    >>> M.adjugate()
+    Matrix([
+    [ 4, -2],
+    [-3,  1]])
+
+    See Also
+    ========
+
+    cofactor_matrix
+    sympy.matrices.common.MatrixCommon.transpose
+    """
+
+    return M.cofactor_matrix(method=method, dotprodsimp=dotprodsimp).transpose()
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _charpoly(M, x='lambda', simplify=_simplify, dotprodsimp=None):
+    """Computes characteristic polynomial det(x*I - M) where I is
+    the identity matrix.
+
+    A PurePoly is returned, so using different variables for ``x`` does
+    not affect the comparison or the polynomials:
+
+    Parameters
+    ==========
+
+    x : string, optional
+        Name for the "lambda" variable, defaults to "lambda".
+
+    simplify : function, optional
+        Simplification function to use on the characteristic polynomial
+        calculated. Defaults to ``simplify``, if ``dotprodsimp`` is ``True``
+        then this is ignored.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        to control expression blowup during matrix multiplication. If this
+        is true then the simplify function is not used.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> from sympy.abc import x, y
+    >>> M = Matrix([[1, 3], [2, 0]])
+    >>> M.charpoly()
+    PurePoly(lambda**2 - lambda - 6, lambda, domain='ZZ')
+    >>> M.charpoly(x) == M.charpoly(y)
+    True
+    >>> M.charpoly(x) == M.charpoly(y)
+    True
+
+    Specifying ``x`` is optional; a symbol named ``lambda`` is used by
+    default (which looks good when pretty-printed in unicode):
+
+    >>> M.charpoly().as_expr()
+    lambda**2 - lambda - 6
+
+    And if ``x`` clashes with an existing symbol, underscores will
+    be prepended to the name to make it unique:
+
+    >>> M = Matrix([[1, 2], [x, 0]])
+    >>> M.charpoly(x).as_expr()
+    _x**2 - _x - 2*x
+
+    Whether you pass a symbol or not, the generator can be obtained
+    with the gen attribute since it may not be the same as the symbol
+    that was passed:
+
+    >>> M.charpoly(x).gen
+    _x
+    >>> M.charpoly(x).gen == x
+    False
+
+    Notes
+    =====
+
+    The Samuelson-Berkowitz algorithm is used to compute
+    the characteristic polynomial efficiently and without any
+    division operations.  Thus the characteristic polynomial over any
+    commutative ring without zero divisors can be computed.
+
+    See Also
+    ========
+
+    det
+    """
+
+    if not M.is_square:
+        raise NonSquareMatrixError()
+
+    if dotprodsimp:
+        simplify = lambda e: e
+
+    berk_vector = _berkowitz_vector(M, dotprodsimp=dotprodsimp)
+    x = _uniquely_named_symbol(x, berk_vector)
+
+    return PurePoly([simplify(a) for a in berk_vector], x)
+
+
+def _cofactor(M, i, j, method="berkowitz", dotprodsimp=None):
+    """Calculate the cofactor of an element.
+
+    Parameters
+    ==========
+
+    method : string, optional
+        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
+        "lu".
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        to control expression blowup during matrix multiplication. If this
+        is true then the simplify function is not used.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2], [3, 4]])
+    >>> M.cofactor(0, 1)
+    -3
+
+    See Also
+    ========
+
+    cofactor_matrix
+    minor
+    minor_submatrix
+    """
+
+    if not M.is_square or M.rows < 1:
+        raise NonSquareMatrixError()
+
+    return (-1)**((i + j) % 2) * M.minor(i, j, method, dotprodsimp=dotprodsimp)
+
+
+def _cofactor_matrix(M, method="berkowitz", dotprodsimp=None):
+    """Return a matrix containing the cofactor of each element.
+
+    Parameters
+    ==========
+
+    method : string, optional
+        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
+        "lu".
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        to control expression blowup during matrix multiplication. If this
+        is true then the simplify function is not used.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2], [3, 4]])
+    >>> M.cofactor_matrix()
+    Matrix([
+    [ 4, -3],
+    [-2,  1]])
+
+    See Also
+    ========
+
+    cofactor
+    minor
+    minor_submatrix
+    """
+
+    if not M.is_square or M.rows < 1:
+        raise NonSquareMatrixError()
+
+    return M._new(M.rows, M.cols,
+            lambda i, j: M.cofactor(i, j, method, dotprodsimp=dotprodsimp))
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _det(M, method="bareiss", iszerofunc=None, dotprodsimp=None):
+    """Computes the determinant of a matrix if ``M`` is a concrete matrix object
+    otherwise return an expressions ``Determinant(M)`` if ``M`` is a
+    ``MatrixSymbol`` or other expression.
+
+    Parameters
+    ==========
+
+    method : string, optional
+        Specifies the algorithm used for computing the matrix determinant.
+
+        If the matrix is at most 3x3, a hard-coded formula is used and the
+        specified method is ignored. Otherwise, it defaults to
+        ``'bareiss'``.
+
+        If it is set to ``'bareiss'``, Bareiss' fraction-free algorithm will
+        be used.
+
+        If it is set to ``'berkowitz'``, Berkowitz' algorithm will be used.
+
+        Otherwise, if it is set to ``'lu'``, LU decomposition will be used.
+
+        .. note::
+            For backward compatibility, legacy keys like "bareis" and
+            "det_lu" can still be used to indicate the corresponding
+            methods.
+            And the keys are also case-insensitive for now. However, it is
+            suggested to use the precise keys for specifying the method.
+
+    iszerofunc : FunctionType or None, optional
+        If it is set to ``None``, it will be defaulted to ``_iszero`` if the
+        method is set to ``'bareiss'``, and ``_is_zero_after_expand_mul`` if
+        the method is set to ``'lu'``.
+
+        It can also accept any user-specified zero testing function, if it
+        is formatted as a function which accepts a single symbolic argument
+        and returns ``True`` if it is tested as zero and ``False`` if it
+        tested as non-zero, and also ``None`` if it is undecidable.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Returns
+    =======
+
+    det : Basic
+        Result of determinant.
+
+    Raises
+    ======
+
+    ValueError
+        If unrecognized keys are given for ``method`` or ``iszerofunc``.
+
+    NonSquareMatrixError
+        If attempted to calculate determinant from a non-square matrix.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix, MatrixSymbol, eye, det
+    >>> M = Matrix([[1, 2], [3, 4]])
+    >>> M.det()
+    -2
+    """
+
+    # sanitize `method`
+    method = method.lower()
+
+    if method == "bareis":
+        method = "bareiss"
+    elif method == "det_lu":
+        method = "lu"
+
+    if method not in ("bareiss", "berkowitz", "lu"):
+        raise ValueError("Determinant method '%s' unrecognized" % method)
+
+    if iszerofunc is None:
+        if method == "bareiss":
+            iszerofunc = _is_zero_after_expand_mul
+        elif method == "lu":
+            iszerofunc = _iszero
+
+    elif not isinstance(iszerofunc, FunctionType):
+        raise ValueError("Zero testing method '%s' unrecognized" % iszerofunc)
+
+    n = M.rows
+
+    if n == M.cols: # square check is done in individual method functions
+        if n == 0:
+            return M.one
+        elif n == 1:
+            return M[0,0]
+        elif n == 2:
+            m = M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]
+            return _dotprodsimp(m) if dotprodsimp else m
+        elif n == 3:
+            m =  (M[0, 0] * M[1, 1] * M[2, 2]
+                + M[0, 1] * M[1, 2] * M[2, 0]
+                + M[0, 2] * M[1, 0] * M[2, 1]
+                - M[0, 2] * M[1, 1] * M[2, 0]
+                - M[0, 0] * M[1, 2] * M[2, 1]
+                - M[0, 1] * M[1, 0] * M[2, 2])
+            return _dotprodsimp(m) if dotprodsimp else m
+
+    if method == "bareiss":
+        return M._eval_det_bareiss(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
+    elif method == "berkowitz":
+        return M._eval_det_berkowitz(dotprodsimp=dotprodsimp)
+    elif method == "lu":
+        return M._eval_det_lu(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
+    else:
+        raise MatrixError('unknown method for calculating determinant')
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _det_bareiss(M, iszerofunc=_is_zero_after_expand_mul, dotprodsimp=None):
+    """Compute matrix determinant using Bareiss' fraction-free
+    algorithm which is an extension of the well known Gaussian
+    elimination method. This approach is best suited for dense
+    symbolic matrices and will result in a determinant with
+    minimal number of fractions. It means that less term
+    rewriting is needed on resulting formulae.
+
+    Parameters
+    ==========
+
+    iszerofunc : function, optional
+        The function to use to determine zeros when doing an LU decomposition.
+        Defaults to ``lambda x: x.is_zero``.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    TODO: Implement algorithm for sparse matrices (SFF),
+    http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps.
+    """
+
+    # Recursively implemented Bareiss' algorithm as per Deanna Richelle Leggett's
+    # thesis http://www.math.usm.edu/perry/Research/Thesis_DRL.pdf
+    def bareiss(mat, cumm=1):
+        if mat.rows == 0:
+            return mat.one
+        elif mat.rows == 1:
+            return mat[0, 0]
+
+        # find a pivot and extract the remaining matrix
+        # With the default iszerofunc, _find_reasonable_pivot slows down
+        # the computation by the factor of 2.5 in one test.
+        # Relevant issues: #10279 and #13877.
+        pivot_pos, pivot_val, _, _ = _find_reasonable_pivot(mat[:, 0], iszerofunc=iszerofunc)
+        if pivot_pos is None:
+            return mat.zero
+
+        # if we have a valid pivot, we'll do a "row swap", so keep the
+        # sign of the det
+        sign = (-1) ** (pivot_pos % 2)
+
+        # we want every row but the pivot row and every column
+        rows = list(i for i in range(mat.rows) if i != pivot_pos)
+        cols = list(range(mat.cols))
+        tmp_mat = mat.extract(rows, cols)
+
+        def entry(i, j):
+            ret = (pivot_val*tmp_mat[i, j + 1] - mat[pivot_pos, j + 1]*tmp_mat[i, 0]) / cumm
+            if dotprodsimp:
+                return _dotprodsimp(ret)
+            elif not ret.is_Atom:
+                return cancel(ret)
+            return ret
+
+        return sign*bareiss(M._new(mat.rows - 1, mat.cols - 1, entry), pivot_val)
+
+    if not M.is_square:
+        raise NonSquareMatrixError()
+
+    if M.rows == 0:
+        return M.one
+        # sympy/matrices/tests/test_matrices.py contains a test that
+        # suggests that the determinant of a 0 x 0 matrix is one, by
+        # convention.
+
+    return bareiss(M)
+
+
+def _det_berkowitz(M, dotprodsimp=None):
+    """ Use the Berkowitz algorithm to compute the determinant.
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+    """
+
+    if not M.is_square:
+        raise NonSquareMatrixError()
+
+    if M.rows == 0:
+        return M.one
+        # sympy/matrices/tests/test_matrices.py contains a test that
+        # suggests that the determinant of a 0 x 0 matrix is one, by
+        # convention.
+
+    berk_vector = _berkowitz_vector(M, dotprodsimp=dotprodsimp)
+    return (-1)**(len(berk_vector) - 1) * berk_vector[-1]
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _det_LU(M, iszerofunc=_iszero, simpfunc=None, dotprodsimp=None):
+    """ Computes the determinant of a matrix from its LU decomposition.
+    This function uses the LU decomposition computed by
+    LUDecomposition_Simple().
+
+    The keyword arguments iszerofunc and simpfunc are passed to
+    LUDecomposition_Simple().
+    iszerofunc is a callable that returns a boolean indicating if its
+    input is zero, or None if it cannot make the determination.
+    simpfunc is a callable that simplifies its input.
+    The default is simpfunc=None, which indicate that the pivot search
+    algorithm should not attempt to simplify any candidate pivots.
+    If simpfunc fails to simplify its input, then it must return its input
+    instead of a copy.
+
+    Parameters
+    ==========
+
+    iszerofunc : function, optional
+        The function to use to determine zeros when doing an LU decomposition.
+        Defaults to ``lambda x: x.is_zero``.
+
+    simpfunc : function, optional
+        The simplification function to use when looking for zeros for pivots.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+    """
+
+    if not M.is_square:
+        raise NonSquareMatrixError()
+
+    if M.rows == 0:
+        return M.one
+        # sympy/matrices/tests/test_matrices.py contains a test that
+        # suggests that the determinant of a 0 x 0 matrix is one, by
+        # convention.
+
+    lu, row_swaps = M.LUdecomposition_Simple(iszerofunc=iszerofunc,
+            simpfunc=simpfunc, dotprodsimp=dotprodsimp)
+    # P*A = L*U => det(A) = det(L)*det(U)/det(P) = det(P)*det(U).
+    # Lower triangular factor L encoded in lu has unit diagonal => det(L) = 1.
+    # P is a permutation matrix => det(P) in {-1, 1} => 1/det(P) = det(P).
+    # LUdecomposition_Simple() returns a list of row exchange index pairs, rather
+    # than a permutation matrix, but det(P) = (-1)**len(row_swaps).
+
+    # Avoid forming the potentially time consuming  product of U's diagonal entries
+    # if the product is zero.
+    # Bottom right entry of U is 0 => det(A) = 0.
+    # It may be impossible to determine if this entry of U is zero when it is symbolic.
+    if iszerofunc(lu[lu.rows-1, lu.rows-1]):
+        return M.zero
+
+    # Compute det(P)
+    det = -M.one if len(row_swaps)%2 else M.one
+
+    # Compute det(U) by calculating the product of U's diagonal entries.
+    # The upper triangular portion of lu is the upper triangular portion of the
+    # U factor in the LU decomposition.
+    for k in range(lu.rows):
+        det *= lu[k, k]
+
+    # return det(P)*det(U)
+    return det
+
+
+def _minor(M, i, j, method="berkowitz", dotprodsimp=None):
+    """Return the (i,j) minor of ``M``.  That is,
+    return the determinant of the matrix obtained by deleting
+    the `i`th row and `j`th column from ``M``.
+
+    Parameters
+    ==========
+
+    i, j : int
+        The row and column to exclude to obtain the submatrix.
+
+    method : string, optional
+        Method to use to find the determinant of the submatrix, can be
+        "bareiss", "berkowitz" or "lu".
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    >>> M.minor(1, 1)
+    -12
+
+    See Also
+    ========
+
+    minor_submatrix
+    cofactor
+    det
+    """
+
+    if not M.is_square:
+        raise NonSquareMatrixError()
+
+    return M.minor_submatrix(i, j).det(method=method, dotprodsimp=dotprodsimp)
+
+
+def _minor_submatrix(M, i, j):
+    """Return the submatrix obtained by removing the `i`th row
+    and `j`th column from ``M`` (works with Pythonic negative indices).
+
+    Parameters
+    ==========
+
+    i, j : int
+        The row and column to exclude to obtain the submatrix.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    >>> M.minor_submatrix(1, 1)
+    Matrix([
+    [1, 3],
+    [7, 9]])
+
+    See Also
+    ========
+
+    minor
+    cofactor
+    """
+
+    if i < 0:
+        i += M.rows
+    if j < 0:
+        j += M.cols
+
+    if not 0 <= i < M.rows or not 0 <= j < M.cols:
+        raise ValueError("`i` and `j` must satisfy 0 <= i < ``M.rows`` "
+                            "(%d)" % M.rows + "and 0 <= j < ``M.cols`` (%d)." % M.cols)
+
+    rows = [a for a in range(M.rows) if a != i]
+    cols = [a for a in range(M.cols) if a != j]
+
+    return M.extract(rows, cols)
diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py
index 487a1bfa3c..cd28826bda 100644
--- a/sympy/matrices/eigen.py
+++ b/sympy/matrices/eigen.py
@@ -1,998 +1,998 @@
-from __future__ import division, print_function
-
-from types import FunctionType
-
-from mpmath.libmp.libmpf import prec_to_dps
-
-from sympy.core.compatibility import default_sort_key
-from sympy.core.logic import fuzzy_and, fuzzy_or
-from sympy.core.numbers import Float
-from sympy.core.symbol import Dummy
-from sympy.functions.elementary.miscellaneous import sqrt
-from sympy.polys import roots
-from sympy.simplify import nsimplify, simplify as _simplify
-from sympy.utilities.exceptions import SymPyDeprecationWarning
-
-from .common import (MatrixError, NonSquareMatrixError,
-    NonPositiveDefiniteMatrixError)
-
-from .utilities import _iszero
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _eigenvals(M, error_when_incomplete=True, dotprodsimp=None, **flags):
-    r"""Return eigenvalues using the Berkowitz agorithm to compute
-    the characteristic polynomial.
-
-    Parameters
-    ==========
-
-    error_when_incomplete : bool, optional
-        If it is set to ``True``, it will raise an error if not all
-        eigenvalues are computed. This is caused by ``roots`` not returning
-        a full list of eigenvalues.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    simplify : bool or function, optional
-        If it is set to ``True``, it attempts to return the most
-        simplified form of expressions returned by applying default
-        simplification method in every routine.
-
-        If it is set to ``False``, it will skip simplification in this
-        particular routine to save computation resources.
-
-        If a function is passed to, it will attempt to apply
-        the particular function as simplification method.
-
-    rational : bool, optional
-        If it is set to ``True``, every floating point numbers would be
-        replaced with rationals before computation. It can solve some
-        issues of ``roots`` routine not working well with floats.
-
-    multiple : bool, optional
-        If it is set to ``True``, the result will be in the form of a
-        list.
-
-        If it is set to ``False``, the result will be in the form of a
-        dictionary.
-
-    Returns
-    =======
-
-    eigs : list or dict
-        Eigenvalues of a matrix. The return format would be specified by
-        the key ``multiple``.
-
-    Raises
-    ======
-
-    MatrixError
-        If not enough roots had got computed.
-
-    NonSquareMatrixError
-        If attempted to compute eigenvalues from a non-square matrix.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1])
-    >>> M.eigenvals()
-    {-1: 1, 0: 1, 2: 1}
-
-    See Also
-    ========
-
-    MatrixDeterminant.charpoly
-    eigenvects
-
-    Notes
-    =====
-
-    Eigenvalues of a matrix `A` can be computed by solving a matrix
-    equation `\det(A - \lambda I) = 0`
-    """
-
-    simplify = flags.get('simplify', False) # Collect simplify flag before popped up, to reuse later in the routine.
-    multiple = flags.get('multiple', False) # Collect multiple flag to decide whether return as a dict or list.
-    rational = flags.pop('rational', True)
-
-    if not M:
-        return {}
-
-    if rational:
-        M = M.applyfunc(
-            lambda x: nsimplify(x, rational=True) if x.has(Float) else x)
-
-    if M.is_upper or M.is_lower:
-        if not M.is_square:
-            raise NonSquareMatrixError()
-
-        diagonal_entries = [M[i, i] for i in range(M.rows)]
-
-        if multiple:
-            eigs = diagonal_entries
-
-        else:
-            eigs = {}
-
-            for diagonal_entry in diagonal_entries:
-                if diagonal_entry not in eigs:
-                    eigs[diagonal_entry] = 0
-
-                eigs[diagonal_entry] += 1
-
-    else:
-        flags.pop('simplify', None)  # pop unsupported flag
-
-        if isinstance(simplify, FunctionType):
-            eigs = roots(M.charpoly(x=Dummy('x'), simplify=simplify,
-                    dotprodsimp=dotprodsimp), **flags)
-        else:
-            eigs = roots(M.charpoly(x=Dummy('x'), dotprodsimp=dotprodsimp), **flags)
-
-    # make sure the algebraic multiplicity sums to the
-    # size of the matrix
-    if error_when_incomplete and (sum(eigs.values()) if
-            isinstance(eigs, dict) else len(eigs)) != M.cols:
-        raise MatrixError("Could not compute eigenvalues for {}".format(M))
-
-    # Since 'simplify' flag is unsupported in roots()
-    # simplify() function will be applied once at the end of the routine.
-    if not simplify:
-        return eigs
-    if not isinstance(simplify, FunctionType):
-        simplify = _simplify
-
-    # With 'multiple' flag set true, simplify() will be mapped for the list
-    # Otherwise, simplify() will be mapped for the keys of the dictionary
-    if not multiple:
-        return {simplify(key): value for key, value in eigs.items()}
-    else:
-        return [simplify(value) for value in eigs]
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero,
-        dotprodsimp=None, **flags):
-    """Return list of triples (eigenval, multiplicity, eigenspace).
-
-    Parameters
-    ==========
-
-    error_when_incomplete : bool, optional
-        Raise an error when not all eigenvalues are computed. This is
-        caused by ``roots`` not returning a full list of eigenvalues.
-
-    iszerofunc : function, optional
-        Specifies a zero testing function to be used in ``rref``.
-
-        Default value is ``_iszero``, which uses SymPy's naive and fast
-        default assumption handler.
-
-        It can also accept any user-specified zero testing function, if it
-        is formatted as a function which accepts a single symbolic argument
-        and returns ``True`` if it is tested as zero and ``False`` if it
-        is tested as non-zero, and ``None`` if it is undecidable.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    simplify : bool or function, optional
-        If ``True``, ``as_content_primitive()`` will be used to tidy up
-        normalization artifacts.
-
-        It will also be used by the ``nullspace`` routine.
-
-    chop : bool or positive number, optional
-        If the matrix contains any Floats, they will be changed to Rationals
-        for computation purposes, but the answers will be returned after
-        being evaluated with evalf. The ``chop`` flag is passed to ``evalf``.
-        When ``chop=True`` a default precision will be used; a number will
-        be interpreted as the desired level of precision.
-
-    Returns
-    =======
-    ret : [(eigenval, multiplicity, eigenspace), ...]
-        A ragged list containing tuples of data obtained by ``eigenvals``
-        and ``nullspace``.
-
-        ``eigenspace`` is a list containing the ``eigenvector`` for each
-        eigenvalue.
-
-        ``eigenvector`` is a vector in the form of a ``Matrix``. e.g.
-        a vector of length 3 is returned as ``Matrix([a_1, a_2, a_3])``.
-
-    Raises
-    ======
-
-    NotImplementedError
-        If failed to compute nullspace.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1])
-    >>> M.eigenvects()
-    [(-1, 1, [Matrix([
-    [-1],
-    [ 1],
-    [ 0]])]), (0, 1, [Matrix([
-    [ 0],
-    [-1],
-    [ 1]])]), (2, 1, [Matrix([
-    [2/3],
-    [1/3],
-    [  1]])])]
-
-    See Also
-    ========
-
-    eigenvals
-    MatrixSubspaces.nullspace
-    """
-
-    def eigenspace(eigenval):
-        """Get a basis for the eigenspace for a particular eigenvalue"""
-
-        m   = M - M.eye(M.rows) * eigenval
-        ret = m.nullspace(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-        # the nullspace for a real eigenvalue should be
-        # non-trivial.  If we didn't find an eigenvector, try once
-        # more a little harder
-        if len(ret) == 0 and simplify:
-            ret = m.nullspace(iszerofunc=iszerofunc, simplify=True, dotprodsimp=dotprodsimp)
-        if len(ret) == 0:
-            raise NotImplementedError(
-                    "Can't evaluate eigenvector for eigenvalue %s" % eigenval)
-
-        return ret
-
-    simplify = flags.get('simplify', True)
-
-    if not isinstance(simplify, FunctionType):
-        simpfunc = _simplify if simplify else lambda x: x
-
-    primitive = flags.get('simplify', False)
-    chop      = flags.pop('chop', False)
-
-    flags.pop('multiple', None)  # remove this if it's there
-
-    has_floats = M.has(Float) # roots doesn't like Floats, so replace them with Rationals
-
-    if has_floats:
-        M = M.applyfunc(lambda x: nsimplify(x, rational=True))
-
-    eigenvals = M.eigenvals(rational=False,
-            error_when_incomplete=error_when_incomplete,
-            dotprodsimp=dotprodsimp, **flags)
-
-    ret = [(val, mult, eigenspace(val)) for val, mult in
-                sorted(eigenvals.items(), key=default_sort_key)]
-
-    if primitive:
-        # if the primitive flag is set, get rid of any common
-        # integer denominators
-        def denom_clean(l):
-            from sympy import gcd
-            return [(v / gcd(list(v))).applyfunc(simpfunc) for v in l]
-
-        ret = [(val, mult, denom_clean(es)) for val, mult, es in ret]
-
-    if has_floats:
-        # if we had floats to start with, turn the eigenvectors to floats
-        ret = [(val.evalf(chop=chop), mult, [v.evalf(chop=chop) for v in es])
-                for val, mult, es in ret]
-
-    return ret
-
-
-def _is_diagonalizable_with_eigen(M, reals_only=False, dotprodsimp=None):
-    """See _is_diagonalizable. This function returns the bool along with the
-    eigenvectors to avoid calculating them again in functions like
-    ``diagonalize``."""
-
-    if not M.is_square:
-        return False, []
-
-    eigenvecs = M.eigenvects(simplify=True, dotprodsimp=dotprodsimp)
-
-    for val, mult, basis in eigenvecs:
-        if reals_only and not val.is_real: # if we have a complex eigenvalue
-            return False, eigenvecs
-
-        if mult != len(basis): # if the geometric multiplicity doesn't equal the algebraic
-            return False, eigenvecs
-
-    return True, eigenvecs
-
-def _is_diagonalizable(M, reals_only=False, dotprodsimp=None, **kwargs):
-    """Returns ``True`` if a matrix is diagonalizable.
-
-    Parameters
-    ==========
-
-    reals_only : bool, optional
-        If ``True``, it tests whether the matrix can be diagonalized
-        to contain only real numbers on the diagonal.
-
-
-        If ``False``, it tests whether the matrix can be diagonalized
-        at all, even with numbers that may not be real.
-
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification
-        is used during matrix multiplications to control expression
-        blowup and thus speed up calculation.
-
-    Examples
-    ========
-
-    Example of a diagonalizable matrix:
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2, 0], [0, 3, 0], [2, -4, 2]])
-    >>> M.is_diagonalizable()
-    True
-
-    Example of a non-diagonalizable matrix:
-
-    >>> M = Matrix([[0, 1], [0, 0]])
-    >>> M.is_diagonalizable()
-    False
-
-    Example of a matrix that is diagonalized in terms of non-real entries:
-
-    >>> M = Matrix([[0, 1], [-1, 0]])
-    >>> M.is_diagonalizable(reals_only=False)
-    True
-    >>> M.is_diagonalizable(reals_only=True)
-    False
-
-    See Also
-    ========
-
-    is_diagonal
-    diagonalize
-    """
-
-    if 'clear_cache' in kwargs:
-        SymPyDeprecationWarning(
-            feature='clear_cache',
-            deprecated_since_version=1.4,
-            issue=15887
-        ).warn()
-
-    if 'clear_subproducts' in kwargs:
-        SymPyDeprecationWarning(
-            feature='clear_subproducts',
-            deprecated_since_version=1.4,
-            issue=15887
-        ).warn()
-
-    if not M.is_square:
-        return False
-
-    if all(e.is_real for e in M) and M.is_symmetric():
-        return True
-
-    if all(e.is_complex for e in M) and M.is_hermitian:
-        return True
-
-    return _is_diagonalizable_with_eigen(M, reals_only=reals_only,
-            dotprodsimp=dotprodsimp)[0]
-
-
-def _diagonalize(M, reals_only=False, sort=False, normalize=False,
-        dotprodsimp=None):
-    """
-    Return (P, D), where D is diagonal and
-
-        D = P^-1 * M * P
-
-    where M is current matrix.
-
-    Parameters
-    ==========
-
-    reals_only : bool. Whether to throw an error if complex numbers are need
-                    to diagonalize. (Default: False)
-
-    sort : bool. Sort the eigenvalues along the diagonal. (Default: False)
-
-    normalize : bool. If True, normalize the columns of P. (Default: False)
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
-    >>> M
-    Matrix([
-    [1,  2, 0],
-    [0,  3, 0],
-    [2, -4, 2]])
-    >>> (P, D) = M.diagonalize()
-    >>> D
-    Matrix([
-    [1, 0, 0],
-    [0, 2, 0],
-    [0, 0, 3]])
-    >>> P
-    Matrix([
-    [-1, 0, -1],
-    [ 0, 0, -1],
-    [ 2, 1,  2]])
-    >>> P.inv() * M * P
-    Matrix([
-    [1, 0, 0],
-    [0, 2, 0],
-    [0, 0, 3]])
-
-    See Also
-    ========
-
-    is_diagonal
-    is_diagonalizable
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    is_diagonalizable, eigenvecs = _is_diagonalizable_with_eigen(M,
-                reals_only=reals_only, dotprodsimp=dotprodsimp)
-
-    if not is_diagonalizable:
-        raise MatrixError("Matrix is not diagonalizable")
-
-    if sort:
-        eigenvecs = sorted(eigenvecs, key=default_sort_key)
-
-    p_cols, diag = [], []
-
-    for val, mult, basis in eigenvecs:
-        diag   += [val] * mult
-        p_cols += basis
-
-    if normalize:
-        p_cols = [v / v.norm() for v in p_cols]
-
-    return M.hstack(*p_cols), M.diag(*diag)
-
-
-def _eval_is_positive_definite(M, method="eigen", dotprodsimp=None):
-    """Algorithm dump for computing positive-definiteness of a
-    matrix.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    method : str, optional
-        Specifies the method for computing positive-definiteness of
-        a matrix.
-
-        If ``'eigen'``, it computes the full eigenvalues and decides
-        if the matrix is positive-definite.
-
-        If ``'CH'``, it attempts computing the Cholesky
-        decomposition to detect the definitiveness.
-
-        If ``'LDL'``, it attempts computing the LDL
-        decomposition to detect the definitiveness.
-    """
-
-    if M.is_hermitian:
-        if method == 'eigen':
-            eigen = M.eigenvals(dotprodsimp=dotprodsimp)
-            args  = [x.is_positive for x in eigen.keys()]
-
-            return fuzzy_and(args)
-
-        elif method == 'CH':
-            try:
-                M.cholesky(hermitian=True)
-            except NonPositiveDefiniteMatrixError:
-                return False
-
-            return True
-
-        elif method == 'LDL':
-            try:
-                M.LDLdecomposition(hermitian=True)
-            except NonPositiveDefiniteMatrixError:
-                return False
-
-            return True
-
-        else:
-            raise NotImplementedError()
-
-    elif M.is_square:
-        M_H = (M + M.H) / 2
-
-        return M_H._eval_is_positive_definite(method=method,
-                dotprodsimp=dotprodsimp)
-
-def _is_positive_definite(M):
-    return M._eval_is_positive_definite()
-
-def _is_positive_semidefinite(M):
-    if M.is_hermitian:
-        eigen = M.eigenvals()
-        args  = [x.is_nonnegative for x in eigen.keys()]
-
-        return fuzzy_and(args)
-
-    elif M.is_square:
-        return ((M + M.H) / 2).is_positive_semidefinite
-
-    return None
-
-def _is_negative_definite(M):
-    if M.is_hermitian:
-        eigen = M.eigenvals()
-        args  = [x.is_negative for x in eigen.keys()]
-
-        return fuzzy_and(args)
-
-    elif M.is_square:
-        return ((M + M.H) / 2).is_negative_definite
-
-    return None
-
-def _is_negative_semidefinite(M):
-    if M.is_hermitian:
-        eigen = M.eigenvals()
-        args  = [x.is_nonpositive for x in eigen.keys()]
-
-        return fuzzy_and(args)
-
-    elif M.is_square:
-        return ((M + M.H) / 2).is_negative_semidefinite
-
-    return None
-
-def _is_indefinite(M):
-    if M.is_hermitian:
-        eigen = M.eigenvals()
-        args1        = [x.is_positive for x in eigen.keys()]
-        any_positive = fuzzy_or(args1)
-        args2        = [x.is_negative for x in eigen.keys()]
-        any_negative = fuzzy_or(args2)
-
-        return fuzzy_and([any_positive, any_negative])
-
-    elif M.is_square:
-        return ((M + M.H) / 2).is_indefinite
-
-    return None
-
-_doc_positive_definite = \
-    r"""Finds out the definiteness of a matrix.
-
-    Examples
-    ========
-
-    An example of numeric positive definite matrix:
-
-    >>> from sympy import Matrix
-    >>> A = Matrix([[1, -2], [-2, 6]])
-    >>> A.is_positive_definite
-    True
-    >>> A.is_positive_semidefinite
-    True
-    >>> A.is_negative_definite
-    False
-    >>> A.is_negative_semidefinite
-    False
-    >>> A.is_indefinite
-    False
-
-    An example of numeric negative definite matrix:
-
-    >>> A = Matrix([[-1, 2], [2, -6]])
-    >>> A.is_positive_definite
-    False
-    >>> A.is_positive_semidefinite
-    False
-    >>> A.is_negative_definite
-    True
-    >>> A.is_negative_semidefinite
-    True
-    >>> A.is_indefinite
-    False
-
-    An example of numeric indefinite matrix:
-
-    >>> A = Matrix([[1, 2], [2, 1]])
-    >>> A.is_positive_definite
-    False
-    >>> A.is_positive_semidefinite
-    False
-    >>> A.is_negative_definite
-    False
-    >>> A.is_negative_semidefinite
-    False
-    >>> A.is_indefinite
-    True
-
-    Notes
-    =====
-
-    Definitiveness is not very commonly discussed for non-hermitian
-    matrices.
-
-    However, computing the definitiveness of a matrix can be
-    generalized over any real matrix by taking the symmetric part:
-
-    `A_S = 1/2 (A + A^{T})`
-
-    Or over any complex matrix by taking the hermitian part:
-
-    `A_H = 1/2 (A + A^{H})`
-
-    And computing the eigenvalues.
-
-    References
-    ==========
-
-    .. [1] https://en.wikipedia.org/wiki/Definiteness_of_a_matrix#Eigenvalues
-
-    .. [2] http://mathworld.wolfram.com/PositiveDefiniteMatrix.html
-
-    .. [3] Johnson, C. R. "Positive Definite Matrices." Amer.
-        Math. Monthly 77, 259-264 1970.
-    """
-
-_is_positive_definite.__doc__     = _doc_positive_definite
-_is_positive_semidefinite.__doc__ = _doc_positive_definite
-_is_negative_definite.__doc__     = _doc_positive_definite
-_is_negative_semidefinite.__doc__ = _doc_positive_definite
-_is_indefinite.__doc__            = _doc_positive_definite
-
-
-def _jordan_form(M, calc_transform=True, dotprodsimp=None, **kwargs):
-    """Return ``(P, J)`` where `J` is a Jordan block
-    matrix and `P` is a matrix such that
-
-        ``M == P*J*P**-1``
-
-    Parameters
-    ==========
-
-    calc_transform : bool
-        If ``False``, then only `J` is returned.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    chop : bool
-        All matrices are converted to exact types when computing
-        eigenvalues and eigenvectors.  As a result, there may be
-        approximation errors.  If ``chop==True``, these errors
-        will be truncated.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix([[ 6,  5, -2, -3], [-3, -1,  3,  3], [ 2,  1, -2, -3], [-1,  1,  5,  5]])
-    >>> P, J = M.jordan_form()
-    >>> J
-    Matrix([
-    [2, 1, 0, 0],
-    [0, 2, 0, 0],
-    [0, 0, 2, 1],
-    [0, 0, 0, 2]])
-
-    See Also
-    ========
-
-    jordan_block
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError("Only square matrices have Jordan forms")
-
-    chop       = kwargs.pop('chop', False)
-    mat        = M
-    has_floats = M.has(Float)
-
-    if has_floats:
-        try:
-            max_prec = max(term._prec for term in M._mat if isinstance(term, Float))
-        except ValueError:
-            # if no term in the matrix is explicitly a Float calling max()
-            # will throw a error so setting max_prec to default value of 53
-            max_prec = 53
-
-        # setting minimum max_dps to 15 to prevent loss of precision in
-        # matrix containing non evaluated expressions
-        max_dps = max(prec_to_dps(max_prec), 15)
-
-    def restore_floats(*args):
-        """If ``has_floats`` is `True`, cast all ``args`` as
-        matrices of floats."""
-
-        if has_floats:
-            args = [m.evalf(n=max_dps, chop=chop) for m in args]
-        if len(args) == 1:
-            return args[0]
-
-        return args
-
-    # cache calculations for some speedup
-    mat_cache = {}
-
-    def eig_mat(val, pow):
-        """Cache computations of ``(M - val*I)**pow`` for quick
-        retrieval"""
-
-        if (val, pow) in mat_cache:
-            return mat_cache[(val, pow)]
-
-        if (val, pow - 1) in mat_cache:
-            mat_cache[(val, pow)] = mat_cache[(val, pow - 1)].multiply(
-                    mat_cache[(val, 1)], dotprodsimp=dotprodsimp)
-        else:
-            mat_cache[(val, pow)] = (mat - val*M.eye(M.rows)).pow(pow,
-                    dotprodsimp=dotprodsimp)
-
-        return mat_cache[(val, pow)]
-
-    # helper functions
-    def nullity_chain(val, algebraic_multiplicity):
-        """Calculate the sequence  [0, nullity(E), nullity(E**2), ...]
-        until it is constant where ``E = M - val*I``"""
-
-        # mat.rank() is faster than computing the null space,
-        # so use the rank-nullity theorem
-        cols    = M.cols
-        ret     = [0]
-        nullity = cols - eig_mat(val, 1).rank(dotprodsimp=dotprodsimp)
-        i       = 2
-
-        while nullity != ret[-1]:
-            ret.append(nullity)
-
-            if nullity == algebraic_multiplicity:
-                break
-
-            nullity  = cols - eig_mat(val, i).rank(dotprodsimp=dotprodsimp)
-            i       += 1
-
-            # Due to issues like #7146 and #15872, SymPy sometimes
-            # gives the wrong rank. In this case, raise an error
-            # instead of returning an incorrect matrix
-            if nullity < ret[-1] or nullity > algebraic_multiplicity:
-                raise MatrixError(
-                    "SymPy had encountered an inconsistent "
-                    "result while computing Jordan block: "
-                    "{}".format(M))
-
-        return ret
-
-    def blocks_from_nullity_chain(d):
-        """Return a list of the size of each Jordan block.
-        If d_n is the nullity of E**n, then the number
-        of Jordan blocks of size n is
-
-            2*d_n - d_(n-1) - d_(n+1)"""
-
-        # d[0] is always the number of columns, so skip past it
-        mid = [2*d[n] - d[n - 1] - d[n + 1] for n in range(1, len(d) - 1)]
-        # d is assumed to plateau with "d[ len(d) ] == d[-1]", so
-        # 2*d_n - d_(n-1) - d_(n+1) == d_n - d_(n-1)
-        end = [d[-1] - d[-2]] if len(d) > 1 else [d[0]]
-
-        return mid + end
-
-    def pick_vec(small_basis, big_basis):
-        """Picks a vector from big_basis that isn't in
-        the subspace spanned by small_basis"""
-
-        if len(small_basis) == 0:
-            return big_basis[0]
-
-        for v in big_basis:
-            _, pivots = M.hstack(*(small_basis + [v])).echelon_form(
-                    with_pivots=True, dotprodsimp=dotprodsimp)
-
-            if pivots[-1] == len(small_basis):
-                return v
-
-    # roots doesn't like Floats, so replace them with Rationals
-    if has_floats:
-        mat = mat.applyfunc(lambda x: nsimplify(x, rational=True))
-
-    # first calculate the jordan block structure
-    eigs = mat.eigenvals(dotprodsimp=dotprodsimp)
-
-    # make sure that we found all the roots by counting
-    # the algebraic multiplicity
-    if sum(m for m in eigs.values()) != mat.cols:
-        raise MatrixError("Could not compute eigenvalues for {}".format(mat))
-
-    # most matrices have distinct eigenvalues
-    # and so are diagonalizable.  In this case, don't
-    # do extra work!
-    if len(eigs.keys()) == mat.cols:
-        blocks     = list(sorted(eigs.keys(), key=default_sort_key))
-        jordan_mat = mat.diag(*blocks)
-
-        if not calc_transform:
-            return restore_floats(jordan_mat)
-
-        jordan_basis = [eig_mat(eig, 1).nullspace(dotprodsimp=dotprodsimp)[0]
-                for eig in blocks]
-        basis_mat    = mat.hstack(*jordan_basis)
-
-        return restore_floats(basis_mat, jordan_mat)
-
-    block_structure = []
-
-    for eig in sorted(eigs.keys(), key=default_sort_key):
-        algebraic_multiplicity = eigs[eig]
-        chain = nullity_chain(eig, algebraic_multiplicity)
-        block_sizes = blocks_from_nullity_chain(chain)
-
-        # if block_sizes =       = [a, b, c, ...], then the number of
-        # Jordan blocks of size 1 is a, of size 2 is b, etc.
-        # create an array that has (eig, block_size) with one
-        # entry for each block
-        size_nums = [(i+1, num) for i, num in enumerate(block_sizes)]
-
-        # we expect larger Jordan blocks to come earlier
-        size_nums.reverse()
-
-        block_structure.extend(
-            (eig, size) for size, num in size_nums for _ in range(num))
-
-    jordan_form_size = sum(size for eig, size in block_structure)
-
-    if jordan_form_size != M.rows:
-        raise MatrixError(
-            "SymPy had encountered an inconsistent result while "
-            "computing Jordan block. : {}".format(M))
-
-    blocks     = (mat.jordan_block(size=size, eigenvalue=eig) for eig, size in block_structure)
-    jordan_mat = mat.diag(*blocks)
-
-    if not calc_transform:
-        return restore_floats(jordan_mat)
-
-    # For each generalized eigenspace, calculate a basis.
-    # We start by looking for a vector in null( (A - eig*I)**n )
-    # which isn't in null( (A - eig*I)**(n-1) ) where n is
-    # the size of the Jordan block
-    #
-    # Ideally we'd just loop through block_structure and
-    # compute each generalized eigenspace.  However, this
-    # causes a lot of unneeded computation.  Instead, we
-    # go through the eigenvalues separately, since we know
-    # their generalized eigenspaces must have bases that
-    # are linearly independent.
-    jordan_basis = []
-
-    for eig in sorted(eigs.keys(), key=default_sort_key):
-        eig_basis = []
-
-        for block_eig, size in block_structure:
-            if block_eig != eig:
-                continue
-
-            null_big   = (eig_mat(eig, size)).nullspace(dotprodsimp=dotprodsimp)
-            null_small = (eig_mat(eig, size - 1)).nullspace(dotprodsimp=dotprodsimp)
-
-            # we want to pick something that is in the big basis
-            # and not the small, but also something that is independent
-            # of any other generalized eigenvectors from a different
-            # generalized eigenspace sharing the same eigenvalue.
-            vec      = pick_vec(null_small + eig_basis, null_big)
-            new_vecs = [eig_mat(eig, i).multiply(vec, dotprodsimp=dotprodsimp)
-                    for i in range(size)]
-
-            eig_basis.extend(new_vecs)
-            jordan_basis.extend(reversed(new_vecs))
-
-    basis_mat = mat.hstack(*jordan_basis)
-
-    return restore_floats(basis_mat, jordan_mat)
-
-
-def _left_eigenvects(M, **flags):
-    """Returns left eigenvectors and eigenvalues.
-
-    This function returns the list of triples (eigenval, multiplicity,
-    basis) for the left eigenvectors. Options are the same as for
-    eigenvects(), i.e. the ``**flags`` arguments gets passed directly to
-    eigenvects().
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix([[0, 1, 1], [1, 0, 0], [1, 1, 1]])
-    >>> M.eigenvects()
-    [(-1, 1, [Matrix([
-    [-1],
-    [ 1],
-    [ 0]])]), (0, 1, [Matrix([
-    [ 0],
-    [-1],
-    [ 1]])]), (2, 1, [Matrix([
-    [2/3],
-    [1/3],
-    [  1]])])]
-    >>> M.left_eigenvects()
-    [(-1, 1, [Matrix([[-2, 1, 1]])]), (0, 1, [Matrix([[-1, -1, 1]])]), (2,
-    1, [Matrix([[1, 1, 1]])])]
-
-    """
-
-    eigs = M.transpose().eigenvects(**flags)
-
-    return [(val, mult, [l.transpose() for l in basis]) for val, mult, basis in eigs]
-
-
-def _singular_values(M, dotprodsimp=None):
-    """Compute the singular values of a Matrix
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix, Symbol
-    >>> x = Symbol('x', real=True)
-    >>> M = Matrix([[0, 1, 0], [0, x, 0], [-1, 0, 0]])
-    >>> M.singular_values()
-    [sqrt(x**2 + 1), 1, 0]
-
-    See Also
-    ========
-
-    condition_number
-    """
-
-    if M.rows >= M.cols:
-        valmultpairs = M.H.multiply(M, dotprodsimp=dotprodsimp) \
-                .eigenvals(dotprodsimp=dotprodsimp)
-    else:
-        valmultpairs = M.multiply(M.H, dotprodsimp=dotprodsimp) \
-                .eigenvals(dotprodsimp=dotprodsimp)
-
-    # Expands result from eigenvals into a simple list
-    vals = []
-
-    for k, v in valmultpairs.items():
-        vals += [sqrt(k)] * v  # dangerous! same k in several spots!
-
-    # Pad with zeros if singular values are computed in reverse way,
-    # to give consistent format.
-    if len(vals) < M.cols:
-        vals += [M.zero] * (M.cols - len(vals))
-
-    # sort them in descending order
-    vals.sort(reverse=True, key=default_sort_key)
-
-    return vals
+from __future__ import division, print_function
+
+from types import FunctionType
+
+from mpmath.libmp.libmpf import prec_to_dps
+
+from sympy.core.compatibility import default_sort_key
+from sympy.core.logic import fuzzy_and, fuzzy_or
+from sympy.core.numbers import Float
+from sympy.core.symbol import Dummy
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.polys import roots
+from sympy.simplify import nsimplify, simplify as _simplify
+from sympy.utilities.exceptions import SymPyDeprecationWarning
+
+from .common import (MatrixError, NonSquareMatrixError,
+    NonPositiveDefiniteMatrixError)
+
+from .utilities import _iszero
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _eigenvals(M, error_when_incomplete=True, dotprodsimp=None, **flags):
+    r"""Return eigenvalues using the Berkowitz agorithm to compute
+    the characteristic polynomial.
+
+    Parameters
+    ==========
+
+    error_when_incomplete : bool, optional
+        If it is set to ``True``, it will raise an error if not all
+        eigenvalues are computed. This is caused by ``roots`` not returning
+        a full list of eigenvalues.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    simplify : bool or function, optional
+        If it is set to ``True``, it attempts to return the most
+        simplified form of expressions returned by applying default
+        simplification method in every routine.
+
+        If it is set to ``False``, it will skip simplification in this
+        particular routine to save computation resources.
+
+        If a function is passed to, it will attempt to apply
+        the particular function as simplification method.
+
+    rational : bool, optional
+        If it is set to ``True``, every floating point numbers would be
+        replaced with rationals before computation. It can solve some
+        issues of ``roots`` routine not working well with floats.
+
+    multiple : bool, optional
+        If it is set to ``True``, the result will be in the form of a
+        list.
+
+        If it is set to ``False``, the result will be in the form of a
+        dictionary.
+
+    Returns
+    =======
+
+    eigs : list or dict
+        Eigenvalues of a matrix. The return format would be specified by
+        the key ``multiple``.
+
+    Raises
+    ======
+
+    MatrixError
+        If not enough roots had got computed.
+
+    NonSquareMatrixError
+        If attempted to compute eigenvalues from a non-square matrix.
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1])
+    >>> M.eigenvals()
+    {-1: 1, 0: 1, 2: 1}
+
+    See Also
+    ========
+
+    MatrixDeterminant.charpoly
+    eigenvects
+
+    Notes
+    =====
+
+    Eigenvalues of a matrix `A` can be computed by solving a matrix
+    equation `\det(A - \lambda I) = 0`
+    """
+
+    simplify = flags.get('simplify', False) # Collect simplify flag before popped up, to reuse later in the routine.
+    multiple = flags.get('multiple', False) # Collect multiple flag to decide whether return as a dict or list.
+    rational = flags.pop('rational', True)
+
+    if not M:
+        return {}
+
+    if rational:
+        M = M.applyfunc(
+            lambda x: nsimplify(x, rational=True) if x.has(Float) else x)
+
+    if M.is_upper or M.is_lower:
+        if not M.is_square:
+            raise NonSquareMatrixError()
+
+        diagonal_entries = [M[i, i] for i in range(M.rows)]
+
+        if multiple:
+            eigs = diagonal_entries
+
+        else:
+            eigs = {}
+
+            for diagonal_entry in diagonal_entries:
+                if diagonal_entry not in eigs:
+                    eigs[diagonal_entry] = 0
+
+                eigs[diagonal_entry] += 1
+
+    else:
+        flags.pop('simplify', None)  # pop unsupported flag
+
+        if isinstance(simplify, FunctionType):
+            eigs = roots(M.charpoly(x=Dummy('x'), simplify=simplify,
+                    dotprodsimp=dotprodsimp), **flags)
+        else:
+            eigs = roots(M.charpoly(x=Dummy('x'), dotprodsimp=dotprodsimp), **flags)
+
+    # make sure the algebraic multiplicity sums to the
+    # size of the matrix
+    if error_when_incomplete and (sum(eigs.values()) if
+            isinstance(eigs, dict) else len(eigs)) != M.cols:
+        raise MatrixError("Could not compute eigenvalues for {}".format(M))
+
+    # Since 'simplify' flag is unsupported in roots()
+    # simplify() function will be applied once at the end of the routine.
+    if not simplify:
+        return eigs
+    if not isinstance(simplify, FunctionType):
+        simplify = _simplify
+
+    # With 'multiple' flag set true, simplify() will be mapped for the list
+    # Otherwise, simplify() will be mapped for the keys of the dictionary
+    if not multiple:
+        return {simplify(key): value for key, value in eigs.items()}
+    else:
+        return [simplify(value) for value in eigs]
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero,
+        dotprodsimp=None, **flags):
+    """Return list of triples (eigenval, multiplicity, eigenspace).
+
+    Parameters
+    ==========
+
+    error_when_incomplete : bool, optional
+        Raise an error when not all eigenvalues are computed. This is
+        caused by ``roots`` not returning a full list of eigenvalues.
+
+    iszerofunc : function, optional
+        Specifies a zero testing function to be used in ``rref``.
+
+        Default value is ``_iszero``, which uses SymPy's naive and fast
+        default assumption handler.
+
+        It can also accept any user-specified zero testing function, if it
+        is formatted as a function which accepts a single symbolic argument
+        and returns ``True`` if it is tested as zero and ``False`` if it
+        is tested as non-zero, and ``None`` if it is undecidable.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    simplify : bool or function, optional
+        If ``True``, ``as_content_primitive()`` will be used to tidy up
+        normalization artifacts.
+
+        It will also be used by the ``nullspace`` routine.
+
+    chop : bool or positive number, optional
+        If the matrix contains any Floats, they will be changed to Rationals
+        for computation purposes, but the answers will be returned after
+        being evaluated with evalf. The ``chop`` flag is passed to ``evalf``.
+        When ``chop=True`` a default precision will be used; a number will
+        be interpreted as the desired level of precision.
+
+    Returns
+    =======
+    ret : [(eigenval, multiplicity, eigenspace), ...]
+        A ragged list containing tuples of data obtained by ``eigenvals``
+        and ``nullspace``.
+
+        ``eigenspace`` is a list containing the ``eigenvector`` for each
+        eigenvalue.
+
+        ``eigenvector`` is a vector in the form of a ``Matrix``. e.g.
+        a vector of length 3 is returned as ``Matrix([a_1, a_2, a_3])``.
+
+    Raises
+    ======
+
+    NotImplementedError
+        If failed to compute nullspace.
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1])
+    >>> M.eigenvects()
+    [(-1, 1, [Matrix([
+    [-1],
+    [ 1],
+    [ 0]])]), (0, 1, [Matrix([
+    [ 0],
+    [-1],
+    [ 1]])]), (2, 1, [Matrix([
+    [2/3],
+    [1/3],
+    [  1]])])]
+
+    See Also
+    ========
+
+    eigenvals
+    MatrixSubspaces.nullspace
+    """
+
+    def eigenspace(eigenval):
+        """Get a basis for the eigenspace for a particular eigenvalue"""
+
+        m   = M - M.eye(M.rows) * eigenval
+        ret = m.nullspace(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
+
+        # the nullspace for a real eigenvalue should be
+        # non-trivial.  If we didn't find an eigenvector, try once
+        # more a little harder
+        if len(ret) == 0 and simplify:
+            ret = m.nullspace(iszerofunc=iszerofunc, simplify=True, dotprodsimp=dotprodsimp)
+        if len(ret) == 0:
+            raise NotImplementedError(
+                    "Can't evaluate eigenvector for eigenvalue %s" % eigenval)
+
+        return ret
+
+    simplify = flags.get('simplify', True)
+
+    if not isinstance(simplify, FunctionType):
+        simpfunc = _simplify if simplify else lambda x: x
+
+    primitive = flags.get('simplify', False)
+    chop      = flags.pop('chop', False)
+
+    flags.pop('multiple', None)  # remove this if it's there
+
+    has_floats = M.has(Float) # roots doesn't like Floats, so replace them with Rationals
+
+    if has_floats:
+        M = M.applyfunc(lambda x: nsimplify(x, rational=True))
+
+    eigenvals = M.eigenvals(rational=False,
+            error_when_incomplete=error_when_incomplete,
+            dotprodsimp=dotprodsimp, **flags)
+
+    ret = [(val, mult, eigenspace(val)) for val, mult in
+                sorted(eigenvals.items(), key=default_sort_key)]
+
+    if primitive:
+        # if the primitive flag is set, get rid of any common
+        # integer denominators
+        def denom_clean(l):
+            from sympy import gcd
+            return [(v / gcd(list(v))).applyfunc(simpfunc) for v in l]
+
+        ret = [(val, mult, denom_clean(es)) for val, mult, es in ret]
+
+    if has_floats:
+        # if we had floats to start with, turn the eigenvectors to floats
+        ret = [(val.evalf(chop=chop), mult, [v.evalf(chop=chop) for v in es])
+                for val, mult, es in ret]
+
+    return ret
+
+
+def _is_diagonalizable_with_eigen(M, reals_only=False, dotprodsimp=None):
+    """See _is_diagonalizable. This function returns the bool along with the
+    eigenvectors to avoid calculating them again in functions like
+    ``diagonalize``."""
+
+    if not M.is_square:
+        return False, []
+
+    eigenvecs = M.eigenvects(simplify=True, dotprodsimp=dotprodsimp)
+
+    for val, mult, basis in eigenvecs:
+        if reals_only and not val.is_real: # if we have a complex eigenvalue
+            return False, eigenvecs
+
+        if mult != len(basis): # if the geometric multiplicity doesn't equal the algebraic
+            return False, eigenvecs
+
+    return True, eigenvecs
+
+def _is_diagonalizable(M, reals_only=False, dotprodsimp=None, **kwargs):
+    """Returns ``True`` if a matrix is diagonalizable.
+
+    Parameters
+    ==========
+
+    reals_only : bool, optional
+        If ``True``, it tests whether the matrix can be diagonalized
+        to contain only real numbers on the diagonal.
+
+
+        If ``False``, it tests whether the matrix can be diagonalized
+        at all, even with numbers that may not be real.
+
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification
+        is used during matrix multiplications to control expression
+        blowup and thus speed up calculation.
+
+    Examples
+    ========
+
+    Example of a diagonalizable matrix:
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2, 0], [0, 3, 0], [2, -4, 2]])
+    >>> M.is_diagonalizable()
+    True
+
+    Example of a non-diagonalizable matrix:
+
+    >>> M = Matrix([[0, 1], [0, 0]])
+    >>> M.is_diagonalizable()
+    False
+
+    Example of a matrix that is diagonalized in terms of non-real entries:
+
+    >>> M = Matrix([[0, 1], [-1, 0]])
+    >>> M.is_diagonalizable(reals_only=False)
+    True
+    >>> M.is_diagonalizable(reals_only=True)
+    False
+
+    See Also
+    ========
+
+    is_diagonal
+    diagonalize
+    """
+
+    if 'clear_cache' in kwargs:
+        SymPyDeprecationWarning(
+            feature='clear_cache',
+            deprecated_since_version=1.4,
+            issue=15887
+        ).warn()
+
+    if 'clear_subproducts' in kwargs:
+        SymPyDeprecationWarning(
+            feature='clear_subproducts',
+            deprecated_since_version=1.4,
+            issue=15887
+        ).warn()
+
+    if not M.is_square:
+        return False
+
+    if all(e.is_real for e in M) and M.is_symmetric():
+        return True
+
+    if all(e.is_complex for e in M) and M.is_hermitian:
+        return True
+
+    return _is_diagonalizable_with_eigen(M, reals_only=reals_only,
+            dotprodsimp=dotprodsimp)[0]
+
+
+def _diagonalize(M, reals_only=False, sort=False, normalize=False,
+        dotprodsimp=None):
+    """
+    Return (P, D), where D is diagonal and
+
+        D = P^-1 * M * P
+
+    where M is current matrix.
+
+    Parameters
+    ==========
+
+    reals_only : bool. Whether to throw an error if complex numbers are need
+                    to diagonalize. (Default: False)
+
+    sort : bool. Sort the eigenvalues along the diagonal. (Default: False)
+
+    normalize : bool. If True, normalize the columns of P. (Default: False)
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
+    >>> M
+    Matrix([
+    [1,  2, 0],
+    [0,  3, 0],
+    [2, -4, 2]])
+    >>> (P, D) = M.diagonalize()
+    >>> D
+    Matrix([
+    [1, 0, 0],
+    [0, 2, 0],
+    [0, 0, 3]])
+    >>> P
+    Matrix([
+    [-1, 0, -1],
+    [ 0, 0, -1],
+    [ 2, 1,  2]])
+    >>> P.inv() * M * P
+    Matrix([
+    [1, 0, 0],
+    [0, 2, 0],
+    [0, 0, 3]])
+
+    See Also
+    ========
+
+    is_diagonal
+    is_diagonalizable
+    """
+
+    if not M.is_square:
+        raise NonSquareMatrixError()
+
+    is_diagonalizable, eigenvecs = _is_diagonalizable_with_eigen(M,
+                reals_only=reals_only, dotprodsimp=dotprodsimp)
+
+    if not is_diagonalizable:
+        raise MatrixError("Matrix is not diagonalizable")
+
+    if sort:
+        eigenvecs = sorted(eigenvecs, key=default_sort_key)
+
+    p_cols, diag = [], []
+
+    for val, mult, basis in eigenvecs:
+        diag   += [val] * mult
+        p_cols += basis
+
+    if normalize:
+        p_cols = [v / v.norm() for v in p_cols]
+
+    return M.hstack(*p_cols), M.diag(*diag)
+
+
+def _eval_is_positive_definite(M, method="eigen", dotprodsimp=None):
+    """Algorithm dump for computing positive-definiteness of a
+    matrix.
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    method : str, optional
+        Specifies the method for computing positive-definiteness of
+        a matrix.
+
+        If ``'eigen'``, it computes the full eigenvalues and decides
+        if the matrix is positive-definite.
+
+        If ``'CH'``, it attempts computing the Cholesky
+        decomposition to detect the definitiveness.
+
+        If ``'LDL'``, it attempts computing the LDL
+        decomposition to detect the definitiveness.
+    """
+
+    if M.is_hermitian:
+        if method == 'eigen':
+            eigen = M.eigenvals(dotprodsimp=dotprodsimp)
+            args  = [x.is_positive for x in eigen.keys()]
+
+            return fuzzy_and(args)
+
+        elif method == 'CH':
+            try:
+                M.cholesky(hermitian=True)
+            except NonPositiveDefiniteMatrixError:
+                return False
+
+            return True
+
+        elif method == 'LDL':
+            try:
+                M.LDLdecomposition(hermitian=True)
+            except NonPositiveDefiniteMatrixError:
+                return False
+
+            return True
+
+        else:
+            raise NotImplementedError()
+
+    elif M.is_square:
+        M_H = (M + M.H) / 2
+
+        return M_H._eval_is_positive_definite(method=method,
+                dotprodsimp=dotprodsimp)
+
+def _is_positive_definite(M):
+    return M._eval_is_positive_definite()
+
+def _is_positive_semidefinite(M):
+    if M.is_hermitian:
+        eigen = M.eigenvals()
+        args  = [x.is_nonnegative for x in eigen.keys()]
+
+        return fuzzy_and(args)
+
+    elif M.is_square:
+        return ((M + M.H) / 2).is_positive_semidefinite
+
+    return None
+
+def _is_negative_definite(M):
+    if M.is_hermitian:
+        eigen = M.eigenvals()
+        args  = [x.is_negative for x in eigen.keys()]
+
+        return fuzzy_and(args)
+
+    elif M.is_square:
+        return ((M + M.H) / 2).is_negative_definite
+
+    return None
+
+def _is_negative_semidefinite(M):
+    if M.is_hermitian:
+        eigen = M.eigenvals()
+        args  = [x.is_nonpositive for x in eigen.keys()]
+
+        return fuzzy_and(args)
+
+    elif M.is_square:
+        return ((M + M.H) / 2).is_negative_semidefinite
+
+    return None
+
+def _is_indefinite(M):
+    if M.is_hermitian:
+        eigen = M.eigenvals()
+        args1        = [x.is_positive for x in eigen.keys()]
+        any_positive = fuzzy_or(args1)
+        args2        = [x.is_negative for x in eigen.keys()]
+        any_negative = fuzzy_or(args2)
+
+        return fuzzy_and([any_positive, any_negative])
+
+    elif M.is_square:
+        return ((M + M.H) / 2).is_indefinite
+
+    return None
+
+_doc_positive_definite = \
+    r"""Finds out the definiteness of a matrix.
+
+    Examples
+    ========
+
+    An example of numeric positive definite matrix:
+
+    >>> from sympy import Matrix
+    >>> A = Matrix([[1, -2], [-2, 6]])
+    >>> A.is_positive_definite
+    True
+    >>> A.is_positive_semidefinite
+    True
+    >>> A.is_negative_definite
+    False
+    >>> A.is_negative_semidefinite
+    False
+    >>> A.is_indefinite
+    False
+
+    An example of numeric negative definite matrix:
+
+    >>> A = Matrix([[-1, 2], [2, -6]])
+    >>> A.is_positive_definite
+    False
+    >>> A.is_positive_semidefinite
+    False
+    >>> A.is_negative_definite
+    True
+    >>> A.is_negative_semidefinite
+    True
+    >>> A.is_indefinite
+    False
+
+    An example of numeric indefinite matrix:
+
+    >>> A = Matrix([[1, 2], [2, 1]])
+    >>> A.is_positive_definite
+    False
+    >>> A.is_positive_semidefinite
+    False
+    >>> A.is_negative_definite
+    False
+    >>> A.is_negative_semidefinite
+    False
+    >>> A.is_indefinite
+    True
+
+    Notes
+    =====
+
+    Definitiveness is not very commonly discussed for non-hermitian
+    matrices.
+
+    However, computing the definitiveness of a matrix can be
+    generalized over any real matrix by taking the symmetric part:
+
+    `A_S = 1/2 (A + A^{T})`
+
+    Or over any complex matrix by taking the hermitian part:
+
+    `A_H = 1/2 (A + A^{H})`
+
+    And computing the eigenvalues.
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Definiteness_of_a_matrix#Eigenvalues
+
+    .. [2] http://mathworld.wolfram.com/PositiveDefiniteMatrix.html
+
+    .. [3] Johnson, C. R. "Positive Definite Matrices." Amer.
+        Math. Monthly 77, 259-264 1970.
+    """
+
+_is_positive_definite.__doc__     = _doc_positive_definite
+_is_positive_semidefinite.__doc__ = _doc_positive_definite
+_is_negative_definite.__doc__     = _doc_positive_definite
+_is_negative_semidefinite.__doc__ = _doc_positive_definite
+_is_indefinite.__doc__            = _doc_positive_definite
+
+
+def _jordan_form(M, calc_transform=True, dotprodsimp=None, **kwargs):
+    """Return ``(P, J)`` where `J` is a Jordan block
+    matrix and `P` is a matrix such that
+
+        ``M == P*J*P**-1``
+
+    Parameters
+    ==========
+
+    calc_transform : bool
+        If ``False``, then only `J` is returned.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    chop : bool
+        All matrices are converted to exact types when computing
+        eigenvalues and eigenvectors.  As a result, there may be
+        approximation errors.  If ``chop==True``, these errors
+        will be truncated.
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix([[ 6,  5, -2, -3], [-3, -1,  3,  3], [ 2,  1, -2, -3], [-1,  1,  5,  5]])
+    >>> P, J = M.jordan_form()
+    >>> J
+    Matrix([
+    [2, 1, 0, 0],
+    [0, 2, 0, 0],
+    [0, 0, 2, 1],
+    [0, 0, 0, 2]])
+
+    See Also
+    ========
+
+    jordan_block
+    """
+
+    if not M.is_square:
+        raise NonSquareMatrixError("Only square matrices have Jordan forms")
+
+    chop       = kwargs.pop('chop', False)
+    mat        = M
+    has_floats = M.has(Float)
+
+    if has_floats:
+        try:
+            max_prec = max(term._prec for term in M._mat if isinstance(term, Float))
+        except ValueError:
+            # if no term in the matrix is explicitly a Float calling max()
+            # will throw a error so setting max_prec to default value of 53
+            max_prec = 53
+
+        # setting minimum max_dps to 15 to prevent loss of precision in
+        # matrix containing non evaluated expressions
+        max_dps = max(prec_to_dps(max_prec), 15)
+
+    def restore_floats(*args):
+        """If ``has_floats`` is `True`, cast all ``args`` as
+        matrices of floats."""
+
+        if has_floats:
+            args = [m.evalf(n=max_dps, chop=chop) for m in args]
+        if len(args) == 1:
+            return args[0]
+
+        return args
+
+    # cache calculations for some speedup
+    mat_cache = {}
+
+    def eig_mat(val, pow):
+        """Cache computations of ``(M - val*I)**pow`` for quick
+        retrieval"""
+
+        if (val, pow) in mat_cache:
+            return mat_cache[(val, pow)]
+
+        if (val, pow - 1) in mat_cache:
+            mat_cache[(val, pow)] = mat_cache[(val, pow - 1)].multiply(
+                    mat_cache[(val, 1)], dotprodsimp=dotprodsimp)
+        else:
+            mat_cache[(val, pow)] = (mat - val*M.eye(M.rows)).pow(pow,
+                    dotprodsimp=dotprodsimp)
+
+        return mat_cache[(val, pow)]
+
+    # helper functions
+    def nullity_chain(val, algebraic_multiplicity):
+        """Calculate the sequence  [0, nullity(E), nullity(E**2), ...]
+        until it is constant where ``E = M - val*I``"""
+
+        # mat.rank() is faster than computing the null space,
+        # so use the rank-nullity theorem
+        cols    = M.cols
+        ret     = [0]
+        nullity = cols - eig_mat(val, 1).rank(dotprodsimp=dotprodsimp)
+        i       = 2
+
+        while nullity != ret[-1]:
+            ret.append(nullity)
+
+            if nullity == algebraic_multiplicity:
+                break
+
+            nullity  = cols - eig_mat(val, i).rank(dotprodsimp=dotprodsimp)
+            i       += 1
+
+            # Due to issues like #7146 and #15872, SymPy sometimes
+            # gives the wrong rank. In this case, raise an error
+            # instead of returning an incorrect matrix
+            if nullity < ret[-1] or nullity > algebraic_multiplicity:
+                raise MatrixError(
+                    "SymPy had encountered an inconsistent "
+                    "result while computing Jordan block: "
+                    "{}".format(M))
+
+        return ret
+
+    def blocks_from_nullity_chain(d):
+        """Return a list of the size of each Jordan block.
+        If d_n is the nullity of E**n, then the number
+        of Jordan blocks of size n is
+
+            2*d_n - d_(n-1) - d_(n+1)"""
+
+        # d[0] is always the number of columns, so skip past it
+        mid = [2*d[n] - d[n - 1] - d[n + 1] for n in range(1, len(d) - 1)]
+        # d is assumed to plateau with "d[ len(d) ] == d[-1]", so
+        # 2*d_n - d_(n-1) - d_(n+1) == d_n - d_(n-1)
+        end = [d[-1] - d[-2]] if len(d) > 1 else [d[0]]
+
+        return mid + end
+
+    def pick_vec(small_basis, big_basis):
+        """Picks a vector from big_basis that isn't in
+        the subspace spanned by small_basis"""
+
+        if len(small_basis) == 0:
+            return big_basis[0]
+
+        for v in big_basis:
+            _, pivots = M.hstack(*(small_basis + [v])).echelon_form(
+                    with_pivots=True, dotprodsimp=dotprodsimp)
+
+            if pivots[-1] == len(small_basis):
+                return v
+
+    # roots doesn't like Floats, so replace them with Rationals
+    if has_floats:
+        mat = mat.applyfunc(lambda x: nsimplify(x, rational=True))
+
+    # first calculate the jordan block structure
+    eigs = mat.eigenvals(dotprodsimp=dotprodsimp)
+
+    # make sure that we found all the roots by counting
+    # the algebraic multiplicity
+    if sum(m for m in eigs.values()) != mat.cols:
+        raise MatrixError("Could not compute eigenvalues for {}".format(mat))
+
+    # most matrices have distinct eigenvalues
+    # and so are diagonalizable.  In this case, don't
+    # do extra work!
+    if len(eigs.keys()) == mat.cols:
+        blocks     = list(sorted(eigs.keys(), key=default_sort_key))
+        jordan_mat = mat.diag(*blocks)
+
+        if not calc_transform:
+            return restore_floats(jordan_mat)
+
+        jordan_basis = [eig_mat(eig, 1).nullspace(dotprodsimp=dotprodsimp)[0]
+                for eig in blocks]
+        basis_mat    = mat.hstack(*jordan_basis)
+
+        return restore_floats(basis_mat, jordan_mat)
+
+    block_structure = []
+
+    for eig in sorted(eigs.keys(), key=default_sort_key):
+        algebraic_multiplicity = eigs[eig]
+        chain = nullity_chain(eig, algebraic_multiplicity)
+        block_sizes = blocks_from_nullity_chain(chain)
+
+        # if block_sizes =       = [a, b, c, ...], then the number of
+        # Jordan blocks of size 1 is a, of size 2 is b, etc.
+        # create an array that has (eig, block_size) with one
+        # entry for each block
+        size_nums = [(i+1, num) for i, num in enumerate(block_sizes)]
+
+        # we expect larger Jordan blocks to come earlier
+        size_nums.reverse()
+
+        block_structure.extend(
+            (eig, size) for size, num in size_nums for _ in range(num))
+
+    jordan_form_size = sum(size for eig, size in block_structure)
+
+    if jordan_form_size != M.rows:
+        raise MatrixError(
+            "SymPy had encountered an inconsistent result while "
+            "computing Jordan block. : {}".format(M))
+
+    blocks     = (mat.jordan_block(size=size, eigenvalue=eig) for eig, size in block_structure)
+    jordan_mat = mat.diag(*blocks)
+
+    if not calc_transform:
+        return restore_floats(jordan_mat)
+
+    # For each generalized eigenspace, calculate a basis.
+    # We start by looking for a vector in null( (A - eig*I)**n )
+    # which isn't in null( (A - eig*I)**(n-1) ) where n is
+    # the size of the Jordan block
+    #
+    # Ideally we'd just loop through block_structure and
+    # compute each generalized eigenspace.  However, this
+    # causes a lot of unneeded computation.  Instead, we
+    # go through the eigenvalues separately, since we know
+    # their generalized eigenspaces must have bases that
+    # are linearly independent.
+    jordan_basis = []
+
+    for eig in sorted(eigs.keys(), key=default_sort_key):
+        eig_basis = []
+
+        for block_eig, size in block_structure:
+            if block_eig != eig:
+                continue
+
+            null_big   = (eig_mat(eig, size)).nullspace(dotprodsimp=dotprodsimp)
+            null_small = (eig_mat(eig, size - 1)).nullspace(dotprodsimp=dotprodsimp)
+
+            # we want to pick something that is in the big basis
+            # and not the small, but also something that is independent
+            # of any other generalized eigenvectors from a different
+            # generalized eigenspace sharing the same eigenvalue.
+            vec      = pick_vec(null_small + eig_basis, null_big)
+            new_vecs = [eig_mat(eig, i).multiply(vec, dotprodsimp=dotprodsimp)
+                    for i in range(size)]
+
+            eig_basis.extend(new_vecs)
+            jordan_basis.extend(reversed(new_vecs))
+
+    basis_mat = mat.hstack(*jordan_basis)
+
+    return restore_floats(basis_mat, jordan_mat)
+
+
+def _left_eigenvects(M, **flags):
+    """Returns left eigenvectors and eigenvalues.
+
+    This function returns the list of triples (eigenval, multiplicity,
+    basis) for the left eigenvectors. Options are the same as for
+    eigenvects(), i.e. the ``**flags`` arguments gets passed directly to
+    eigenvects().
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix([[0, 1, 1], [1, 0, 0], [1, 1, 1]])
+    >>> M.eigenvects()
+    [(-1, 1, [Matrix([
+    [-1],
+    [ 1],
+    [ 0]])]), (0, 1, [Matrix([
+    [ 0],
+    [-1],
+    [ 1]])]), (2, 1, [Matrix([
+    [2/3],
+    [1/3],
+    [  1]])])]
+    >>> M.left_eigenvects()
+    [(-1, 1, [Matrix([[-2, 1, 1]])]), (0, 1, [Matrix([[-1, -1, 1]])]), (2,
+    1, [Matrix([[1, 1, 1]])])]
+
+    """
+
+    eigs = M.transpose().eigenvects(**flags)
+
+    return [(val, mult, [l.transpose() for l in basis]) for val, mult, basis in eigs]
+
+
+def _singular_values(M, dotprodsimp=None):
+    """Compute the singular values of a Matrix
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix, Symbol
+    >>> x = Symbol('x', real=True)
+    >>> M = Matrix([[0, 1, 0], [0, x, 0], [-1, 0, 0]])
+    >>> M.singular_values()
+    [sqrt(x**2 + 1), 1, 0]
+
+    See Also
+    ========
+
+    condition_number
+    """
+
+    if M.rows >= M.cols:
+        valmultpairs = M.H.multiply(M, dotprodsimp=dotprodsimp) \
+                .eigenvals(dotprodsimp=dotprodsimp)
+    else:
+        valmultpairs = M.multiply(M.H, dotprodsimp=dotprodsimp) \
+                .eigenvals(dotprodsimp=dotprodsimp)
+
+    # Expands result from eigenvals into a simple list
+    vals = []
+
+    for k, v in valmultpairs.items():
+        vals += [sqrt(k)] * v  # dangerous! same k in several spots!
+
+    # Pad with zeros if singular values are computed in reverse way,
+    # to give consistent format.
+    if len(vals) < M.cols:
+        vals += [M.zero] * (M.cols - len(vals))
+
+    # sort them in descending order
+    vals.sort(reverse=True, key=default_sort_key)
+
+    return vals
diff --git a/sympy/matrices/inverse.py b/sympy/matrices/inverse.py
deleted file mode 100644
index e4c369c60c..0000000000
--- a/sympy/matrices/inverse.py
+++ /dev/null
@@ -1,498 +0,0 @@
-from __future__ import division, print_function
-
-from sympy.core.numbers import mod_inverse
-
-from .common import MatrixError, NonSquareMatrixError, NonInvertibleMatrixError
-from .utilities import _iszero
-
-
-def _pinv_full_rank(M, dotprodsimp=None):
-    """Subroutine for full row or column rank matrices.
-
-    For full row rank matrices, inverse of ``A * A.H`` Exists.
-    For full column rank matrices, inverse of ``A.H * A`` Exists.
-
-    This routine can apply for both cases by checking the shape
-    and have small decision.
-    """
-
-    if M.is_zero_matrix:
-        return M.H
-
-    if M.rows >= M.cols:
-        return M.H.multiply(M, dotprodsimp=dotprodsimp) \
-                .inv(dotprodsimp=dotprodsimp) \
-                .multiply(M.H, dotprodsimp=dotprodsimp)
-    else:
-        return M.H.multiply(M.multiply(M.H, dotprodsimp=dotprodsimp)
-                .inv(dotprodsimp=dotprodsimp), dotprodsimp=dotprodsimp)
-
-def _pinv_rank_decomposition(M, dotprodsimp=None):
-    """Subroutine for rank decomposition
-
-    With rank decompositions, `A` can be decomposed into two full-
-    rank matrices, and each matrix can take pseudoinverse
-    individually.
-    """
-
-    if M.is_zero_matrix:
-        return M.H
-
-    B, C = M.rank_decomposition(dotprodsimp=dotprodsimp)
-
-    Bp = _pinv_full_rank(B, dotprodsimp=dotprodsimp)
-    Cp = _pinv_full_rank(C, dotprodsimp=dotprodsimp)
-
-    return Cp.multiply(Bp, dotprodsimp=dotprodsimp)
-
-def _pinv_diagonalization(M, dotprodsimp=None):
-    """Subroutine using diagonalization
-
-    This routine can sometimes fail if SymPy's eigenvalue
-    computation is not reliable.
-    """
-
-    if M.is_zero_matrix:
-        return M.H
-
-    A  = M
-    AH = M.H
-
-    try:
-        if M.rows >= M.cols:
-            P, D   = AH.multiply(A, dotprodsimp=dotprodsimp).diagonalize(
-                        normalize=True, dotprodsimp=dotprodsimp)
-            D_pinv = D.applyfunc(lambda x: 0 if _iszero(x) else 1 / x)
-
-            return P.multiply(D_pinv, dotprodsimp=dotprodsimp) \
-                    .multiply(P.H, dotprodsimp=dotprodsimp) \
-                    .multiply(AH, dotprodsimp=dotprodsimp)
-
-        else:
-            P, D   = A.multiply(AH, dotprodsimp=dotprodsimp).diagonalize(
-                        normalize=True, dotprodsimp=dotprodsimp)
-            D_pinv = D.applyfunc(lambda x: 0 if _iszero(x) else 1 / x)
-
-            return AH.multiply(P, dotprodsimp=dotprodsimp) \
-                    .multiply(D_pinv, dotprodsimp=dotprodsimp) \
-                    .multiply(P.H, dotprodsimp=dotprodsimp)
-
-    except MatrixError:
-        raise NotImplementedError(
-            'pinv for rank-deficient matrices where '
-            'diagonalization of A.H*A fails is not supported yet.')
-
-def _pinv(M, method='RD', dotprodsimp=None):
-    """Calculate the Moore-Penrose pseudoinverse of the matrix.
-
-    The Moore-Penrose pseudoinverse exists and is unique for any matrix.
-    If the matrix is invertible, the pseudoinverse is the same as the
-    inverse.
-
-    Parameters
-    ==========
-
-    method : String, optional
-        Specifies the method for computing the pseudoinverse.
-
-        If ``'RD'``, Rank-Decomposition will be used.
-
-        If ``'ED'``, Diagonalization will be used.
-
-    Examples
-    ========
-
-    Computing pseudoinverse by rank decomposition :
-
-    >>> from sympy import Matrix
-    >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
-    >>> A.pinv()
-    Matrix([
-    [-17/18,  4/9],
-    [  -1/9,  1/9],
-    [ 13/18, -2/9]])
-
-    Computing pseudoinverse by diagonalization :
-
-    >>> B = A.pinv(method='ED')
-    >>> B.simplify()
-    >>> B
-    Matrix([
-    [-17/18,  4/9],
-    [  -1/9,  1/9],
-    [ 13/18, -2/9]])
-
-    See Also
-    ========
-
-    inv
-    pinv_solve
-
-    References
-    ==========
-
-    .. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse
-
-    """
-
-    # Trivial case: pseudoinverse of all-zero matrix is its transpose.
-    if M.is_zero_matrix:
-        return M.H
-
-    if method == 'RD':
-        return _pinv_rank_decomposition(M, dotprodsimp=dotprodsimp)
-    elif method == 'ED':
-        return _pinv_diagonalization(M, dotprodsimp=dotprodsimp)
-    else:
-        raise ValueError('invalid pinv method %s' % repr(method))
-
-
-def _inv_mod(M, m):
-    r"""
-    Returns the inverse of the matrix `K` (mod `m`), if it exists.
-
-    Method to find the matrix inverse of `K` (mod `m`) implemented in this function:
-
-    * Compute `\mathrm{adj}(K) = \mathrm{cof}(K)^t`, the adjoint matrix of `K`.
-
-    * Compute `r = 1/\mathrm{det}(K) \pmod m`.
-
-    * `K^{-1} = r\cdot \mathrm{adj}(K) \pmod m`.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> A = Matrix(2, 2, [1, 2, 3, 4])
-    >>> A.inv_mod(5)
-    Matrix([
-    [3, 1],
-    [4, 2]])
-    >>> A.inv_mod(3)
-    Matrix([
-    [1, 1],
-    [0, 1]])
-
-    """
-
-    if not M.is_square:
-        raise NonSquareMatrixError()
-
-    N       = M.cols
-    det_K   = M.det()
-    det_inv = None
-
-    try:
-        det_inv = mod_inverse(det_K, m)
-    except ValueError:
-        raise NonInvertibleMatrixError('Matrix is not invertible (mod %d)' % m)
-
-    K_adj = M.adjugate()
-    K_inv = M.__class__(N, N,
-            [det_inv * K_adj[i, j] % m for i in range(N) for j in range(N)])
-
-    return K_inv
-
-
-def _verify_invertible(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Initial check to see if a matrix is invertible. Raises or returns
-    determinant for use in _inv_ADJ."""
-
-    if not M.is_square:
-        raise NonSquareMatrixError("A Matrix must be square to invert.")
-
-    d    = M.det(method='berkowitz', dotprodsimp=dotprodsimp)
-    zero = d.equals(0)
-
-    if zero is None: # if equals() can't decide, will rref be able to?
-        ok   = M.rref(simplify=True, dotprodsimp=dotprodsimp)[0]
-        zero = any(iszerofunc(ok[j, j]) for j in range(ok.rows))
-
-    if zero:
-        raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
-
-    return d
-
-def _inv_ADJ(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Calculates the inverse using the adjugate matrix and a determinant.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    See Also
-    ========
-
-    inv
-    inverse_GE
-    inverse_LU
-    inverse_CH
-    inverse_LDL
-    """
-
-    d = _verify_invertible(M, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    return M.adjugate(dotprodsimp=dotprodsimp) / d
-
-def _inv_GE(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Calculates the inverse using Gaussian elimination.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    See Also
-    ========
-
-    inv
-    inverse_ADJ
-    inverse_LU
-    inverse_CH
-    inverse_LDL
-    """
-
-    from .dense import Matrix
-
-    if not M.is_square:
-        raise NonSquareMatrixError("A Matrix must be square to invert.")
-
-    big = Matrix.hstack(M.as_mutable(), Matrix.eye(M.rows))
-    red = big.rref(iszerofunc=iszerofunc, simplify=True, dotprodsimp=dotprodsimp)[0]
-
-    if any(iszerofunc(red[j, j]) for j in range(red.rows)):
-        raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
-
-    return M._new(red[:, big.rows:])
-
-def _inv_LU(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Calculates the inverse using LU decomposition.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    See Also
-    ========
-
-    inv
-    inverse_ADJ
-    inverse_GE
-    inverse_CH
-    inverse_LDL
-    """
-
-    _verify_invertible(M, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    return M.LUsolve(M.eye(M.rows), iszerofunc=_iszero,
-            dotprodsimp=dotprodsimp)
-
-def _inv_CH(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Calculates the inverse using cholesky decomposition.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    See Also
-    ========
-
-    inv
-    inverse_ADJ
-    inverse_GE
-    inverse_LU
-    inverse_LDL
-    """
-
-    _verify_invertible(M, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    return M.cholesky_solve(M.eye(M.rows), dotprodsimp=dotprodsimp)
-
-def _inv_LDL(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Calculates the inverse using LDL decomposition.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    See Also
-    ========
-
-    inv
-    inverse_ADJ
-    inverse_GE
-    inverse_LU
-    inverse_CH
-    """
-
-    _verify_invertible(M, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    return M.LDLsolve(M.eye(M.rows), dotprodsimp=dotprodsimp)
-
-def _inv_QR(M, iszerofunc=_iszero, dotprodsimp=None):
-    """Calculates the inverse using QR decomposition.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    See Also
-    ========
-
-    inv
-    inverse_ADJ
-    inverse_GE
-    inverse_CH
-    inverse_LDL
-    """
-
-    _verify_invertible(M, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    return M.QRsolve(M.eye(M.rows), dotprodsimp=dotprodsimp)
-
-def _inv(M, method=None, iszerofunc=_iszero, try_block_diag=False,
-        dotprodsimp=None):
-    """
-    Return the inverse of a matrix using the method indicated. Default for
-    dense matrices is is Gauss elimination, default for sparse matrices is LDL.
-
-    Parameters
-    ==========
-
-    method : ('GE', 'LU', 'ADJ', 'CH', 'LDL')
-
-    iszerofunc : function, optional
-        Zero-testing function to use.
-
-    try_block_diag : bool, optional
-        If True then will try to form block diagonal matrices using the
-        method get_diag_blocks(), invert these individually, and then
-        reconstruct the full inverse matrix.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy import SparseMatrix, Matrix
-    >>> A = SparseMatrix([
-    ... [ 2, -1,  0],
-    ... [-1,  2, -1],
-    ... [ 0,  0,  2]])
-    >>> A.inv('CH')
-    Matrix([
-    [2/3, 1/3, 1/6],
-    [1/3, 2/3, 1/3],
-    [  0,   0, 1/2]])
-    >>> A.inv(method='LDL') # use of 'method=' is optional
-    Matrix([
-    [2/3, 1/3, 1/6],
-    [1/3, 2/3, 1/3],
-    [  0,   0, 1/2]])
-    >>> A * _
-    Matrix([
-    [1, 0, 0],
-    [0, 1, 0],
-    [0, 0, 1]])
-    >>> A = Matrix(A)
-    >>> A.inv('CH')
-    Matrix([
-    [2/3, 1/3, 1/6],
-    [1/3, 2/3, 1/3],
-    [  0,   0, 1/2]])
-    >>> A.inv('ADJ') == A.inv('GE') == A.inv('LU') == A.inv('CH') == A.inv('LDL') == A.inv('QR')
-    True
-
-    Notes
-    =====
-
-    According to the ``method`` keyword, it calls the appropriate method:
-
-        GE .... inverse_GE(); default for dense matrices
-        LU .... inverse_LU()
-        ADJ ... inverse_ADJ()
-        CH ... inverse_CH()
-        LDL ... inverse_LDL(); default for sparse matrices
-        QR ... inverse_QR()
-
-    Note, the GE and LU methods may require the matrix to be simplified
-    before it is inverted in order to properly detect zeros during
-    pivoting. In difficult cases a custom zero detection function can
-    be provided by setting the ``iszerofunc`` argument to a function that
-    should return True if its argument is zero. The ADJ routine computes
-    the determinant and uses that to detect singular matrices in addition
-    to testing for zeros on the diagonal.
-
-    See Also
-    ========
-
-    inverse_ADJ
-    inverse_GE
-    inverse_LU
-    inverse_CH
-    inverse_LDL
-
-    Raises
-    ======
-
-    ValueError
-        If the determinant of the matrix is zero.
-    """
-
-    from sympy.matrices import diag, SparseMatrix
-
-    if method is None:
-        method = 'LDL' if isinstance(M, SparseMatrix) else 'GE'
-
-    if try_block_diag:
-        blocks = M.get_diag_blocks()
-        r      = []
-
-        for block in blocks:
-            r.append(block.inv(method=method, iszerofunc=iszerofunc,
-                    dotprodsimp=dotprodsimp))
-
-        return diag(*r)
-
-    if method == "GE":
-        rv = M.inverse_GE(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    elif method == "LU":
-        rv = M.inverse_LU(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    elif method == "ADJ":
-        rv = M.inverse_ADJ(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    elif method == "CH":
-        rv = M.inverse_CH(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    elif method == "LDL":
-        rv = M.inverse_LDL(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    elif method == "QR":
-        rv = M.inverse_QR(iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-    else:
-        raise ValueError("Inversion method unrecognized")
-
-    return M._new(rv)
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 373f7563b1..5241bcdfce 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -8,9 +8,10 @@
     Callable, NotIterable, as_int, is_sequence)
 from sympy.core.decorators import deprecated
 from sympy.core.expr import Expr
+from sympy.core.numbers import mod_inverse
 from sympy.core.power import Pow
 from sympy.core.singleton import S
-from sympy.core.symbol import Dummy, Symbol, _uniquely_named_symbol
+from sympy.core.symbol import Dummy, Symbol, _uniquely_named_symbol, symbols
 from sympy.core.sympify import sympify
 from sympy.functions import exp, factorial, log
 from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
@@ -50,12 +51,7 @@
 
 from .solvers import (
     _diagonal_solve, _lower_triangular_solve, _upper_triangular_solve,
-    _cholesky_solve, _LDLsolve, _LUsolve, _QRsolve, _gauss_jordan_solve,
-    _pinv_solve, _solve, _solve_least_squares)
-
-from .inverse import (
-    _pinv, _inv_mod, _inv_ADJ, _inv_GE, _inv_LU, _inv_CH, _inv_LDL, _inv_QR,
-    _inv)
+    _cholesky_solve, _LDLsolve, _LUsolve, _QRsolve, _gauss_jordan_solve)
 
 
 class DeferredVector(Symbol, NotIterable):
@@ -1677,6 +1673,211 @@ def log(self, simplify=cancel):
 
         return ret
 
+    def inv_mod(self, m):
+        r"""
+        Returns the inverse of the matrix `K` (mod `m`), if it exists.
+
+        Method to find the matrix inverse of `K` (mod `m`) implemented in this function:
+
+        * Compute `\mathrm{adj}(K) = \mathrm{cof}(K)^t`, the adjoint matrix of `K`.
+
+        * Compute `r = 1/\mathrm{det}(K) \pmod m`.
+
+        * `K^{-1} = r\cdot \mathrm{adj}(K) \pmod m`.
+
+        Examples
+        ========
+
+        >>> from sympy import Matrix
+        >>> A = Matrix(2, 2, [1, 2, 3, 4])
+        >>> A.inv_mod(5)
+        Matrix([
+        [3, 1],
+        [4, 2]])
+        >>> A.inv_mod(3)
+        Matrix([
+        [1, 1],
+        [0, 1]])
+
+        """
+        if not self.is_square:
+            raise NonSquareMatrixError()
+        N = self.cols
+        det_K = self.det()
+        det_inv = None
+
+        try:
+            det_inv = mod_inverse(det_K, m)
+        except ValueError:
+            raise NonInvertibleMatrixError('Matrix is not invertible (mod %d)' % m)
+
+        K_adj = self.adjugate()
+        K_inv = self.__class__(N, N,
+                               [det_inv * K_adj[i, j] % m for i in range(N) for
+                                j in range(N)])
+        return K_inv
+
+    def inverse_ADJ(self, iszerofunc=_iszero, dotprodsimp=None):
+        """Calculates the inverse using the adjugate matrix and a determinant.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        See Also
+        ========
+
+        inv
+        inverse_LU
+        inverse_GE
+        """
+        if not self.is_square:
+            raise NonSquareMatrixError("A Matrix must be square to invert.")
+
+        d = self.det(method='berkowitz', dotprodsimp=dotprodsimp)
+        zero = d.equals(0)
+        if zero is None:
+            # if equals() can't decide, will rref be able to?
+            ok = self.rref(simplify=True, dotprodsimp=dotprodsimp)[0]
+            zero = any(iszerofunc(ok[j, j]) for j in range(ok.rows))
+        if zero:
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
+
+        return self.adjugate(dotprodsimp=dotprodsimp) / d
+
+    def inverse_GE(self, iszerofunc=_iszero, dotprodsimp=None):
+        """Calculates the inverse using Gaussian elimination.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        See Also
+        ========
+
+        inv
+        inverse_LU
+        inverse_ADJ
+        """
+        from .dense import Matrix
+        if not self.is_square:
+            raise NonSquareMatrixError("A Matrix must be square to invert.")
+
+        big = Matrix.hstack(self.as_mutable(), Matrix.eye(self.rows))
+        red = big.rref(iszerofunc=iszerofunc, simplify=True, dotprodsimp=dotprodsimp)[0]
+        if any(iszerofunc(red[j, j]) for j in range(red.rows)):
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
+
+        return self._new(red[:, big.rows:])
+
+    def inverse_LU(self, iszerofunc=_iszero, dotprodsimp=None):
+        """Calculates the inverse using LU decomposition.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        See Also
+        ========
+
+        inv
+        inverse_GE
+        inverse_ADJ
+        """
+        if not self.is_square:
+            raise NonSquareMatrixError()
+
+        ok = self.rref(simplify=True, dotprodsimp=dotprodsimp)[0]
+        if any(iszerofunc(ok[j, j]) for j in range(ok.rows)):
+            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
+
+        return self.LUsolve(self.eye(self.rows), iszerofunc=_iszero,
+                dotprodsimp=dotprodsimp)
+
+    def inv(self, method=None, dotprodsimp=None, **kwargs):
+        """
+        Return the inverse of a matrix.
+
+        CASE 1: If the matrix is a dense matrix.
+
+        Return the matrix inverse using the method indicated (default
+        is Gauss elimination).
+
+        Parameters
+        ==========
+
+        method : ('GE', 'LU', or 'ADJ') for dense matrices,
+                 ('LDL', 'CH') for sparse
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        Notes
+        =====
+
+        According to the ``method`` keyword, it calls the appropriate method:
+
+          GE .... inverse_GE(); default
+          LU .... inverse_LU()
+          ADJ ... inverse_ADJ()
+
+        See Also
+        ========
+
+        inverse_LU
+        inverse_GE
+        inverse_ADJ
+
+        Raises
+        ------
+        ValueError
+            If the determinant of the matrix is zero.
+
+        CASE 2: If the matrix is a sparse matrix.
+
+        Return the matrix inverse using Cholesky or LDL (default).
+
+        kwargs
+        ======
+
+        method : ('CH', 'LDL')
+
+        Notes
+        =====
+
+        According to the ``method`` keyword, it calls the appropriate method:
+
+          LDL ... inverse_LDL(); default
+          CH .... inverse_CH()
+
+        Raises
+        ------
+        ValueError
+            If the determinant of the matrix is zero.
+
+        """
+        if not self.is_square:
+            raise NonSquareMatrixError()
+        if method is not None:
+            kwargs['method'] = method
+        if dotprodsimp is not None:
+            kwargs['dotprodsimp'] = dotprodsimp
+        return self._eval_inverse(**kwargs)
+
     def is_nilpotent(self, dotprodsimp=None):
         """Checks if a matrix is nilpotent.
 
@@ -1906,6 +2107,216 @@ def norm(self, ord=None):
             else:
                 raise NotImplementedError("Matrix Norms under development")
 
+    def pinv_solve(self, B, arbitrary_matrix=None):
+        """Solve ``Ax = B`` using the Moore-Penrose pseudoinverse.
+
+        There may be zero, one, or infinite solutions.  If one solution
+        exists, it will be returned.  If infinite solutions exist, one will
+        be returned based on the value of arbitrary_matrix.  If no solutions
+        exist, the least-squares solution is returned.
+
+        Parameters
+        ==========
+
+        B : Matrix
+            The right hand side of the equation to be solved for.  Must have
+            the same number of rows as matrix A.
+        arbitrary_matrix : Matrix
+            If the system is underdetermined (e.g. A has more columns than
+            rows), infinite solutions are possible, in terms of an arbitrary
+            matrix.  This parameter may be set to a specific matrix to use
+            for that purpose; if so, it must be the same shape as x, with as
+            many rows as matrix A has columns, and as many columns as matrix
+            B.  If left as None, an appropriate matrix containing dummy
+            symbols in the form of ``wn_m`` will be used, with n and m being
+            row and column position of each symbol.
+
+        Returns
+        =======
+
+        x : Matrix
+            The matrix that will satisfy ``Ax = B``.  Will have as many rows as
+            matrix A has columns, and as many columns as matrix B.
+
+        Examples
+        ========
+
+        >>> from sympy import Matrix
+        >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
+        >>> B = Matrix([7, 8])
+        >>> A.pinv_solve(B)
+        Matrix([
+        [ _w0_0/6 - _w1_0/3 + _w2_0/6 - 55/18],
+        [-_w0_0/3 + 2*_w1_0/3 - _w2_0/3 + 1/9],
+        [ _w0_0/6 - _w1_0/3 + _w2_0/6 + 59/18]])
+        >>> A.pinv_solve(B, arbitrary_matrix=Matrix([0, 0, 0]))
+        Matrix([
+        [-55/18],
+        [   1/9],
+        [ 59/18]])
+
+        See Also
+        ========
+
+        sympy.matrices.dense.DenseMatrix.lower_triangular_solve
+        sympy.matrices.dense.DenseMatrix.upper_triangular_solve
+        gauss_jordan_solve
+        cholesky_solve
+        diagonal_solve
+        LDLsolve
+        LUsolve
+        QRsolve
+        pinv
+
+        Notes
+        =====
+
+        This may return either exact solutions or least squares solutions.
+        To determine which, check ``A * A.pinv() * B == B``.  It will be
+        True if exact solutions exist, and False if only a least-squares
+        solution exists.  Be aware that the left hand side of that equation
+        may need to be simplified to correctly compare to the right hand
+        side.
+
+        References
+        ==========
+
+        .. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#Obtaining_all_solutions_of_a_linear_system
+
+        """
+        from sympy.matrices import eye
+        A = self
+        A_pinv = self.pinv()
+        if arbitrary_matrix is None:
+            rows, cols = A.cols, B.cols
+            w = symbols('w:{0}_:{1}'.format(rows, cols), cls=Dummy)
+            arbitrary_matrix = self.__class__(cols, rows, w).T
+        return A_pinv * B + (eye(A.cols) - A_pinv * A) * arbitrary_matrix
+
+    def _eval_pinv_full_rank(self):
+        """Subroutine for full row or column rank matrices.
+
+        For full row rank matrices, inverse of ``A * A.H`` Exists.
+        For full column rank matrices, inverse of ``A.H * A`` Exists.
+
+        This routine can apply for both cases by checking the shape
+        and have small decision.
+        """
+        if self.is_zero_matrix:
+            return self.H
+
+        if self.rows >= self.cols:
+            return (self.H * self).inv() * self.H
+        else:
+            return self.H * (self * self.H).inv()
+
+    def _eval_pinv_rank_decomposition(self):
+        """Subroutine for rank decomposition
+
+        With rank decompositions, `A` can be decomposed into two full-
+        rank matrices, and each matrix can take pseudoinverse
+        individually.
+        """
+        if self.is_zero_matrix:
+            return self.H
+
+        B, C = self.rank_decomposition()
+
+        Bp = B._eval_pinv_full_rank()
+        Cp = C._eval_pinv_full_rank()
+
+        return Cp * Bp
+
+    def _eval_pinv_diagonalization(self):
+        """Subroutine using diagonalization
+
+        This routine can sometimes fail if SymPy's eigenvalue
+        computation is not reliable.
+        """
+        if self.is_zero_matrix:
+            return self.H
+
+        A = self
+        AH = self.H
+
+        try:
+            if self.rows >= self.cols:
+                P, D = (AH * A).diagonalize(normalize=True)
+                D_pinv = D.applyfunc(lambda x: 0 if _iszero(x) else 1 / x)
+                return P * D_pinv * P.H * AH
+            else:
+                P, D = (A * AH).diagonalize(normalize=True)
+                D_pinv = D.applyfunc(lambda x: 0 if _iszero(x) else 1 / x)
+                return AH * P * D_pinv * P.H
+        except MatrixError:
+            raise NotImplementedError(
+                'pinv for rank-deficient matrices where '
+                'diagonalization of A.H*A fails is not supported yet.')
+
+
+    def pinv(self, method='RD'):
+        """Calculate the Moore-Penrose pseudoinverse of the matrix.
+
+        The Moore-Penrose pseudoinverse exists and is unique for any matrix.
+        If the matrix is invertible, the pseudoinverse is the same as the
+        inverse.
+
+        Parameters
+        ==========
+
+        method : String, optional
+            Specifies the method for computing the pseudoinverse.
+
+            If ``'RD'``, Rank-Decomposition will be used.
+
+            If ``'ED'``, Diagonalization will be used.
+
+        Examples
+        ========
+
+        Computing pseudoinverse by rank decomposition :
+
+        >>> from sympy import Matrix
+        >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
+        >>> A.pinv()
+        Matrix([
+        [-17/18,  4/9],
+        [  -1/9,  1/9],
+        [ 13/18, -2/9]])
+
+        Computing pseudoinverse by diagonalization :
+
+        >>> B = A.pinv(method='ED')
+        >>> B.simplify()
+        >>> B
+        Matrix([
+        [-17/18,  4/9],
+        [  -1/9,  1/9],
+        [ 13/18, -2/9]])
+
+        See Also
+        ========
+
+        inv
+        pinv_solve
+
+        References
+        ==========
+
+        .. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse
+
+        """
+        # Trivial case: pseudoinverse of all-zero matrix is its transpose.
+        if self.is_zero_matrix:
+            return self.H
+
+        if method == 'RD':
+            return self._eval_pinv_rank_decomposition()
+        elif method == 'ED':
+            return self._eval_pinv_diagonalization()
+        else:
+            raise ValueError()
+
     def print_nonzero(self, symb="X"):
         """Shows location of non-zero entries for fast shape lookup.
 
@@ -1956,6 +2367,168 @@ def project(self, v):
         """
         return v * (self.dot(v) / v.dot(v))
 
+    def solve_least_squares(self, rhs, method='CH'):
+        """Return the least-square fit to the data.
+
+        Parameters
+        ==========
+
+        rhs : Matrix
+            Vector representing the right hand side of the linear equation.
+
+        method : string or boolean, optional
+            If set to ``'CH'``, ``cholesky_solve`` routine will be used.
+
+            If set to ``'LDL'``, ``LDLsolve`` routine will be used.
+
+            If set to ``'QR'``, ``QRsolve`` routine will be used.
+
+            If set to ``'PINV'``, ``pinv_solve`` routine will be used.
+
+            Otherwise, the conjugate of ``self`` will be used to create a system
+            of equations that is passed to ``solve`` along with the hint
+            defined by ``method``.
+
+        Returns
+        =======
+
+        solutions : Matrix
+            Vector representing the solution.
+
+        Examples
+        ========
+
+        >>> from sympy.matrices import Matrix, ones
+        >>> A = Matrix([1, 2, 3])
+        >>> B = Matrix([2, 3, 4])
+        >>> S = Matrix(A.row_join(B))
+        >>> S
+        Matrix([
+        [1, 2],
+        [2, 3],
+        [3, 4]])
+
+        If each line of S represent coefficients of Ax + By
+        and x and y are [2, 3] then S*xy is:
+
+        >>> r = S*Matrix([2, 3]); r
+        Matrix([
+        [ 8],
+        [13],
+        [18]])
+
+        But let's add 1 to the middle value and then solve for the
+        least-squares value of xy:
+
+        >>> xy = S.solve_least_squares(Matrix([8, 14, 18])); xy
+        Matrix([
+        [ 5/3],
+        [10/3]])
+
+        The error is given by S*xy - r:
+
+        >>> S*xy - r
+        Matrix([
+        [1/3],
+        [1/3],
+        [1/3]])
+        >>> _.norm().n(2)
+        0.58
+
+        If a different xy is used, the norm will be higher:
+
+        >>> xy += ones(2, 1)/10
+        >>> (S*xy - r).norm().n(2)
+        1.5
+
+        """
+        if method == 'CH':
+            return self.cholesky_solve(rhs)
+        elif method == 'QR':
+            return self.QRsolve(rhs)
+        elif method == 'LDL':
+            return self.LDLsolve(rhs)
+        elif method == 'PINV':
+            return self.pinv_solve(rhs)
+        else:
+            t = self.H
+            return (t * self).solve(t * rhs, method=method)
+
+    def solve(self, rhs, method='GJ', dotprodsimp=None):
+        """Solves linear equation where the unique solution exists.
+
+        Parameters
+        ==========
+
+        rhs : Matrix
+            Vector representing the right hand side of the linear equation.
+
+        method : string, optional
+           If set to ``'GJ'``, the Gauss-Jordan elimination will be used, which
+           is implemented in the routine ``gauss_jordan_solve``.
+
+           If set to ``'LU'``, ``LUsolve`` routine will be used.
+
+           If set to ``'QR'``, ``QRsolve`` routine will be used.
+
+           If set to ``'PINV'``, ``pinv_solve`` routine will be used.
+
+           It also supports the methods available for special linear systems
+
+           For positive definite systems:
+
+           If set to ``'CH'``, ``cholesky_solve`` routine will be used.
+
+           If set to ``'LDL'``, ``LDLsolve`` routine will be used.
+
+           To use a different method and to compute the solution via the
+           inverse, use a method defined in the .inv() docstring.
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        Returns
+        =======
+
+        solutions : Matrix
+            Vector representing the solution.
+
+        Raises
+        ======
+
+        ValueError
+            If there is not a unique solution then a ``ValueError`` will be
+            raised.
+
+            If ``self`` is not square, a ``ValueError`` and a different routine
+            for solving the system will be suggested.
+        """
+
+        if method == 'GJ':
+            try:
+                soln, param = self.gauss_jordan_solve(rhs, dotprodsimp=dotprodsimp)
+                if param:
+                    raise NonInvertibleMatrixError("Matrix det == 0; not invertible. "
+                    "Try ``self.gauss_jordan_solve(rhs)`` to obtain a parametric solution.")
+            except ValueError:
+                raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
+            return soln
+        elif method == 'LU':
+            return self.LUsolve(rhs, dotprodsimp=dotprodsimp)
+        elif method == 'CH':
+            return self.cholesky_solve(rhs, dotprodsimp=dotprodsimp)
+        elif method == 'QR':
+            return self.QRsolve(rhs)
+        elif method == 'LDL':
+            return self.LDLsolve(rhs, dotprodsimp=dotprodsimp)
+        elif method == 'PINV':
+            return self.pinv_solve(rhs)
+        else:
+            return self.inv(method=method, dotprodsimp=dotprodsimp) \
+                    .multiply(rhs, dotprodsimp=dotprodsimp)
+
     def table(self, printer, rowstart='[', rowend=']', rowsep='\n',
               colsep=', ', align='right'):
         r"""
@@ -2089,10 +2662,9 @@ def vech(self, diagonal=True, check_symmetry=True):
                     count += 1
         return v
 
-    def rank_decomposition(self, iszerofunc=_iszero, simplify=False,
-            dotprodsimp=None):
+    def rank_decomposition(self, iszerofunc=_iszero, simplify=False):
         return _rank_decomposition(self, iszerofunc=iszerofunc,
-                simplify=simplify, dotprodsimp=dotprodsimp)
+                simplify=simplify)
 
     def cholesky(self, hermitian=True, dotprodsimp=None):
         raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix')
@@ -2142,45 +2714,6 @@ def gauss_jordan_solve(self, B, freevar=False, dotprodsimp=None):
         return _gauss_jordan_solve(self, B, freevar=freevar,
                 dotprodsimp=dotprodsimp)
 
-    def pinv_solve(self, B, arbitrary_matrix=None, dotprodsimp=None):
-        return _pinv_solve(self, B, arbitrary_matrix=arbitrary_matrix,
-                dotprodsimp=dotprodsimp)
-
-    def solve(self, rhs, method='GJ', dotprodsimp=None):
-        return _solve(self, rhs, method=method, dotprodsimp=dotprodsimp)
-
-    def solve_least_squares(self, rhs, method='CH'):
-        return _solve_least_squares(self, rhs, method=method)
-
-    def pinv(self, method='RD', dotprodsimp=None):
-        return _pinv(self, method=method, dotprodsimp=dotprodsimp)
-
-    def inv_mod(self, m):
-        return _inv_mod(self, m)
-
-    def inverse_ADJ(self, iszerofunc=_iszero, dotprodsimp=None):
-        return _inv_ADJ(self, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    def inverse_GE(self, iszerofunc=_iszero, dotprodsimp=None):
-        return _inv_GE(self, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    def inverse_LU(self, iszerofunc=_iszero, dotprodsimp=None):
-        return _inv_LU(self, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    def inverse_CH(self, iszerofunc=_iszero, dotprodsimp=None):
-        return _inv_CH(self, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    def inverse_LDL(self, iszerofunc=_iszero, dotprodsimp=None):
-        return _inv_LDL(self, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    def inverse_QR(self, iszerofunc=_iszero, dotprodsimp=None):
-        return _inv_QR(self, iszerofunc=iszerofunc, dotprodsimp=dotprodsimp)
-
-    def inv(self, method=None, iszerofunc=_iszero, try_block_diag=False,
-            dotprodsimp=None):
-        return _inv(self, method=method, iszerofunc=iszerofunc,
-                try_block_diag=try_block_diag, dotprodsimp=dotprodsimp)
-
     rank_decomposition.__doc__     = _rank_decomposition.__doc__
     cholesky.__doc__               = _cholesky.__doc__
     LDLdecomposition.__doc__       = _LDLdecomposition.__doc__
@@ -2197,19 +2730,6 @@ def inv(self, method=None, iszerofunc=_iszero, try_block_diag=False,
     LUsolve.__doc__                = _LUsolve.__doc__
     QRsolve.__doc__                = _QRsolve.__doc__
     gauss_jordan_solve.__doc__     = _gauss_jordan_solve.__doc__
-    pinv_solve.__doc__             = _pinv_solve.__doc__
-    solve.__doc__                  = _solve.__doc__
-    solve_least_squares.__doc__    = _solve_least_squares.__doc__
-
-    pinv.__doc__                   = _pinv.__doc__
-    inv_mod.__doc__                = _inv_mod.__doc__
-    inverse_ADJ.__doc__            = _inv_ADJ.__doc__
-    inverse_GE.__doc__             = _inv_GE.__doc__
-    inverse_LU.__doc__             = _inv_LU.__doc__
-    inverse_CH.__doc__             = _inv_CH.__doc__
-    inverse_LDL.__doc__            = _inv_LDL.__doc__
-    inverse_QR.__doc__             = _inv_QR.__doc__
-    inv.__doc__                    = _inv.__doc__
 
 
 @deprecated(
diff --git a/sympy/matrices/reductions.py b/sympy/matrices/reductions.py
index e021b397c4..53436c8163 100644
--- a/sympy/matrices/reductions.py
+++ b/sympy/matrices/reductions.py
@@ -1,335 +1,335 @@
-from __future__ import division, print_function
-
-from types import FunctionType
-
-from sympy.simplify.simplify import (
-    simplify as _simplify, dotprodsimp as _dotprodsimp)
-
-from .utilities import _iszero
-from .determinant import _find_reasonable_pivot
-
-
-def _row_reduce_list(mat, rows, cols, one, iszerofunc, simpfunc,
-                normalize_last=True, normalize=True, zero_above=True,
-                dotprodsimp=None):
-    """Row reduce a flat list representation of a matrix and return a tuple
-    (rref_matrix, pivot_cols, swaps) where ``rref_matrix`` is a flat list,
-    ``pivot_cols`` are the pivot columns and ``swaps`` are any row swaps that
-    were used in the process of row reduction.
-
-    Parameters
-    ==========
-
-    mat : list
-        list of matrix elements, must be ``rows`` * ``cols`` in length
-
-    rows, cols : integer
-        number of rows and columns in flat list representation
-
-    one : SymPy object
-        represents the value one, from ``Matrix.one``
-
-    iszerofunc : determines if an entry can be used as a pivot
-
-    simpfunc : used to simplify elements and test if they are
-        zero if ``iszerofunc`` returns `None`
-
-    normalize_last : indicates where all row reduction should
-        happen in a fraction-free manner and then the rows are
-        normalized (so that the pivots are 1), or whether
-        rows should be normalized along the way (like the naive
-        row reduction algorithm)
-
-    normalize : whether pivot rows should be normalized so that
-        the pivot value is 1
-
-    zero_above : whether entries above the pivot should be zeroed.
-        If ``zero_above=False``, an echelon matrix will be returned.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-    """
-
-    def get_col(i):
-        return mat[i::cols]
-
-    def row_swap(i, j):
-        mat[i*cols:(i + 1)*cols], mat[j*cols:(j + 1)*cols] = \
-            mat[j*cols:(j + 1)*cols], mat[i*cols:(i + 1)*cols]
-
-    def cross_cancel(a, i, b, j):
-        """Does the row op row[i] = a*row[i] - b*row[j]"""
-        q = (j - i)*cols
-        for p in range(i*cols, (i + 1)*cols):
-            mat[p] = dps(a*mat[p] - b*mat[p + q])
-
-    dps = _dotprodsimp if dotprodsimp else lambda e: e
-    piv_row, piv_col = 0, 0
-    pivot_cols = []
-    swaps = []
-
-    # use a fraction free method to zero above and below each pivot
-    while piv_col < cols and piv_row < rows:
-        pivot_offset, pivot_val, \
-        assumed_nonzero, newly_determined = _find_reasonable_pivot(
-                get_col(piv_col)[piv_row:], iszerofunc, simpfunc)
-
-        # _find_reasonable_pivot may have simplified some things
-        # in the process.  Let's not let them go to waste
-        for (offset, val) in newly_determined:
-            offset += piv_row
-            mat[offset*cols + piv_col] = val
-
-        if pivot_offset is None:
-            piv_col += 1
-            continue
-
-        pivot_cols.append(piv_col)
-        if pivot_offset != 0:
-            row_swap(piv_row, pivot_offset + piv_row)
-            swaps.append((piv_row, pivot_offset + piv_row))
-
-        # if we aren't normalizing last, we normalize
-        # before we zero the other rows
-        if normalize_last is False:
-            i, j = piv_row, piv_col
-            mat[i*cols + j] = one
-            for p in range(i*cols + j + 1, (i + 1)*cols):
-                mat[p] = dps(mat[p] / pivot_val)
-            # after normalizing, the pivot value is 1
-            pivot_val = one
-
-        # zero above and below the pivot
-        for row in range(rows):
-            # don't zero our current row
-            if row == piv_row:
-                continue
-            # don't zero above the pivot unless we're told.
-            if zero_above is False and row < piv_row:
-                continue
-            # if we're already a zero, don't do anything
-            val = mat[row*cols + piv_col]
-            if iszerofunc(val):
-                continue
-
-            cross_cancel(pivot_val, row, val, piv_row)
-        piv_row += 1
-
-    # normalize each row
-    if normalize_last is True and normalize is True:
-        for piv_i, piv_j in enumerate(pivot_cols):
-            pivot_val = mat[piv_i*cols + piv_j]
-            mat[piv_i*cols + piv_j] = one
-            for p in range(piv_i*cols + piv_j + 1, (piv_i + 1)*cols):
-                mat[p] = dps(mat[p] / pivot_val)
-
-    return mat, tuple(pivot_cols), tuple(swaps)
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _row_reduce(M, iszerofunc, simpfunc, normalize_last=True,
-                normalize=True, zero_above=True, dotprodsimp=None):
-
-    mat, pivot_cols, swaps = _row_reduce_list(list(M), M.rows, M.cols, M.one,
-            iszerofunc, simpfunc, normalize_last=normalize_last,
-            normalize=normalize, zero_above=zero_above, dotprodsimp=dotprodsimp)
-
-    return M._new(M.rows, M.cols, mat), pivot_cols, swaps
-
-
-def _is_echelon(M, iszerofunc=_iszero):
-    """Returns `True` if the matrix is in echelon form. That is, all rows of
-    zeros are at the bottom, and below each leading non-zero in a row are
-    exclusively zeros."""
-
-    if M.rows <= 0 or M.cols <= 0:
-        return True
-
-    zeros_below = all(iszerofunc(t) for t in M[1:, 0])
-
-    if iszerofunc(M[0, 0]):
-        return zeros_below and _is_echelon(M[:, 1:], iszerofunc)
-
-    return zeros_below and _is_echelon(M[1:, 1:], iszerofunc)
-
-
-def _echelon_form(M, iszerofunc=_iszero, simplify=False, with_pivots=False,
-        dotprodsimp=None):
-    """Returns a matrix row-equivalent to ``M`` that is in echelon form. Note
-    that echelon form of a matrix is *not* unique, however, properties like the
-    row space and the null space are preserved.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix([[1, 2], [3, 4]])
-    >>> M.echelon_form()
-    Matrix([
-    [1,  2],
-    [0, -2]])
-    """
-
-    simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
-
-    mat, pivots, _ = _row_reduce(M, iszerofunc, simpfunc,
-            normalize_last=True, normalize=False, zero_above=False,
-            dotprodsimp=dotprodsimp)
-
-    if with_pivots:
-        return mat, pivots
-
-    return mat
-
-
-# This functions is a candidate for caching if it gets implemented for matrices.
-def _rank(M, iszerofunc=_iszero, simplify=False, dotprodsimp=None):
-    """Returns the rank of a matrix.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> from sympy.abc import x
-    >>> m = Matrix([[1, 2], [x, 1 - 1/x]])
-    >>> m.rank()
-    2
-    >>> n = Matrix(3, 3, range(1, 10))
-    >>> n.rank()
-    2
-    """
-
-    def _permute_complexity_right(M, iszerofunc):
-        """Permute columns with complicated elements as
-        far right as they can go.  Since the ``sympy`` row reduction
-        algorithms start on the left, having complexity right-shifted
-        speeds things up.
-
-        Returns a tuple (mat, perm) where perm is a permutation
-        of the columns to perform to shift the complex columns right, and mat
-        is the permuted matrix."""
-
-        def complexity(i):
-            # the complexity of a column will be judged by how many
-            # element's zero-ness cannot be determined
-            return sum(1 if iszerofunc(e) is None else 0 for e in M[:, i])
-
-        complex = [(complexity(i), i) for i in range(M.cols)]
-        perm    = [j for (i, j) in sorted(complex)]
-
-        return (M.permute(perm, orientation='cols'), perm)
-
-    simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
-
-    # for small matrices, we compute the rank explicitly
-    # if is_zero on elements doesn't answer the question
-    # for small matrices, we fall back to the full routine.
-    if M.rows <= 0 or M.cols <= 0:
-        return 0
-
-    if M.rows <= 1 or M.cols <= 1:
-        zeros = [iszerofunc(x) for x in M]
-
-        if False in zeros:
-            return 1
-
-    if M.rows == 2 and M.cols == 2:
-        zeros = [iszerofunc(x) for x in M]
-
-        if not False in zeros and not None in zeros:
-            return 0
-
-        d = M.det(dotprodsimp=dotprodsimp)
-
-        if iszerofunc(d) and False in zeros:
-            return 1
-        if iszerofunc(d) is False:
-            return 2
-
-    mat, _       = _permute_complexity_right(M, iszerofunc=iszerofunc)
-    _, pivots, _ = _row_reduce(mat, iszerofunc, simpfunc, normalize_last=True,
-            normalize=False, zero_above=False, dotprodsimp=dotprodsimp)
-
-    return len(pivots)
-
-
-def _rref(M, iszerofunc=_iszero, simplify=False, pivots=True,
-        normalize_last=True, dotprodsimp=None):
-    """Return reduced row-echelon form of matrix and indices of pivot vars.
-
-    Parameters
-    ==========
-
-    iszerofunc : Function
-        A function used for detecting whether an element can
-        act as a pivot.  ``lambda x: x.is_zero`` is used by default.
-
-    simplify : Function
-        A function used to simplify elements when looking for a pivot.
-        By default SymPy's ``simplify`` is used.
-
-    pivots : True or False
-        If ``True``, a tuple containing the row-reduced matrix and a tuple
-        of pivot columns is returned.  If ``False`` just the row-reduced
-        matrix is returned.
-
-    normalize_last : True or False
-        If ``True``, no pivots are normalized to `1` until after all
-        entries above and below each pivot are zeroed.  This means the row
-        reduction algorithm is fraction free until the very last step.
-        If ``False``, the naive row reduction procedure is used where
-        each pivot is normalized to be `1` before row operations are
-        used to zero above and below the pivot.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> from sympy.abc import x
-    >>> m = Matrix([[1, 2], [x, 1 - 1/x]])
-    >>> m.rref()
-    (Matrix([
-    [1, 0],
-    [0, 1]]), (0, 1))
-    >>> rref_matrix, rref_pivots = m.rref()
-    >>> rref_matrix
-    Matrix([
-    [1, 0],
-    [0, 1]])
-    >>> rref_pivots
-    (0, 1)
-
-    Notes
-    =====
-
-    The default value of ``normalize_last=True`` can provide significant
-    speedup to row reduction, especially on matrices with symbols.  However,
-    if you depend on the form row reduction algorithm leaves entries
-    of the matrix, set ``noramlize_last=False``
-    """
-
-    simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
-
-    mat, pivot_cols, _ = _row_reduce(M, iszerofunc, simpfunc,
-            normalize_last, normalize=True, zero_above=True,
-            dotprodsimp=dotprodsimp)
-
-    if pivots:
-        mat = (mat, pivot_cols)
-
-    return mat
+from __future__ import division, print_function
+
+from types import FunctionType
+
+from sympy.simplify.simplify import (
+    simplify as _simplify, dotprodsimp as _dotprodsimp)
+
+from .utilities import _iszero
+from .determinant import _find_reasonable_pivot
+
+
+def _row_reduce_list(mat, rows, cols, one, iszerofunc, simpfunc,
+                normalize_last=True, normalize=True, zero_above=True,
+                dotprodsimp=None):
+    """Row reduce a flat list representation of a matrix and return a tuple
+    (rref_matrix, pivot_cols, swaps) where ``rref_matrix`` is a flat list,
+    ``pivot_cols`` are the pivot columns and ``swaps`` are any row swaps that
+    were used in the process of row reduction.
+
+    Parameters
+    ==========
+
+    mat : list
+        list of matrix elements, must be ``rows`` * ``cols`` in length
+
+    rows, cols : integer
+        number of rows and columns in flat list representation
+
+    one : SymPy object
+        represents the value one, from ``Matrix.one``
+
+    iszerofunc : determines if an entry can be used as a pivot
+
+    simpfunc : used to simplify elements and test if they are
+        zero if ``iszerofunc`` returns `None`
+
+    normalize_last : indicates where all row reduction should
+        happen in a fraction-free manner and then the rows are
+        normalized (so that the pivots are 1), or whether
+        rows should be normalized along the way (like the naive
+        row reduction algorithm)
+
+    normalize : whether pivot rows should be normalized so that
+        the pivot value is 1
+
+    zero_above : whether entries above the pivot should be zeroed.
+        If ``zero_above=False``, an echelon matrix will be returned.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+    """
+
+    def get_col(i):
+        return mat[i::cols]
+
+    def row_swap(i, j):
+        mat[i*cols:(i + 1)*cols], mat[j*cols:(j + 1)*cols] = \
+            mat[j*cols:(j + 1)*cols], mat[i*cols:(i + 1)*cols]
+
+    def cross_cancel(a, i, b, j):
+        """Does the row op row[i] = a*row[i] - b*row[j]"""
+        q = (j - i)*cols
+        for p in range(i*cols, (i + 1)*cols):
+            mat[p] = dps(a*mat[p] - b*mat[p + q])
+
+    dps = _dotprodsimp if dotprodsimp else lambda e: e
+    piv_row, piv_col = 0, 0
+    pivot_cols = []
+    swaps = []
+
+    # use a fraction free method to zero above and below each pivot
+    while piv_col < cols and piv_row < rows:
+        pivot_offset, pivot_val, \
+        assumed_nonzero, newly_determined = _find_reasonable_pivot(
+                get_col(piv_col)[piv_row:], iszerofunc, simpfunc)
+
+        # _find_reasonable_pivot may have simplified some things
+        # in the process.  Let's not let them go to waste
+        for (offset, val) in newly_determined:
+            offset += piv_row
+            mat[offset*cols + piv_col] = val
+
+        if pivot_offset is None:
+            piv_col += 1
+            continue
+
+        pivot_cols.append(piv_col)
+        if pivot_offset != 0:
+            row_swap(piv_row, pivot_offset + piv_row)
+            swaps.append((piv_row, pivot_offset + piv_row))
+
+        # if we aren't normalizing last, we normalize
+        # before we zero the other rows
+        if normalize_last is False:
+            i, j = piv_row, piv_col
+            mat[i*cols + j] = one
+            for p in range(i*cols + j + 1, (i + 1)*cols):
+                mat[p] = dps(mat[p] / pivot_val)
+            # after normalizing, the pivot value is 1
+            pivot_val = one
+
+        # zero above and below the pivot
+        for row in range(rows):
+            # don't zero our current row
+            if row == piv_row:
+                continue
+            # don't zero above the pivot unless we're told.
+            if zero_above is False and row < piv_row:
+                continue
+            # if we're already a zero, don't do anything
+            val = mat[row*cols + piv_col]
+            if iszerofunc(val):
+                continue
+
+            cross_cancel(pivot_val, row, val, piv_row)
+        piv_row += 1
+
+    # normalize each row
+    if normalize_last is True and normalize is True:
+        for piv_i, piv_j in enumerate(pivot_cols):
+            pivot_val = mat[piv_i*cols + piv_j]
+            mat[piv_i*cols + piv_j] = one
+            for p in range(piv_i*cols + piv_j + 1, (piv_i + 1)*cols):
+                mat[p] = dps(mat[p] / pivot_val)
+
+    return mat, tuple(pivot_cols), tuple(swaps)
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _row_reduce(M, iszerofunc, simpfunc, normalize_last=True,
+                normalize=True, zero_above=True, dotprodsimp=None):
+
+    mat, pivot_cols, swaps = _row_reduce_list(list(M), M.rows, M.cols, M.one,
+            iszerofunc, simpfunc, normalize_last=normalize_last,
+            normalize=normalize, zero_above=zero_above, dotprodsimp=dotprodsimp)
+
+    return M._new(M.rows, M.cols, mat), pivot_cols, swaps
+
+
+def _is_echelon(M, iszerofunc=_iszero):
+    """Returns `True` if the matrix is in echelon form. That is, all rows of
+    zeros are at the bottom, and below each leading non-zero in a row are
+    exclusively zeros."""
+
+    if M.rows <= 0 or M.cols <= 0:
+        return True
+
+    zeros_below = all(iszerofunc(t) for t in M[1:, 0])
+
+    if iszerofunc(M[0, 0]):
+        return zeros_below and _is_echelon(M[:, 1:], iszerofunc)
+
+    return zeros_below and _is_echelon(M[1:, 1:], iszerofunc)
+
+
+def _echelon_form(M, iszerofunc=_iszero, simplify=False, with_pivots=False,
+        dotprodsimp=None):
+    """Returns a matrix row-equivalent to ``M`` that is in echelon form. Note
+    that echelon form of a matrix is *not* unique, however, properties like the
+    row space and the null space are preserved.
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix([[1, 2], [3, 4]])
+    >>> M.echelon_form()
+    Matrix([
+    [1,  2],
+    [0, -2]])
+    """
+
+    simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
+
+    mat, pivots, _ = _row_reduce(M, iszerofunc, simpfunc,
+            normalize_last=True, normalize=False, zero_above=False,
+            dotprodsimp=dotprodsimp)
+
+    if with_pivots:
+        return mat, pivots
+
+    return mat
+
+
+# This functions is a candidate for caching if it gets implemented for matrices.
+def _rank(M, iszerofunc=_iszero, simplify=False, dotprodsimp=None):
+    """Returns the rank of a matrix.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> from sympy.abc import x
+    >>> m = Matrix([[1, 2], [x, 1 - 1/x]])
+    >>> m.rank()
+    2
+    >>> n = Matrix(3, 3, range(1, 10))
+    >>> n.rank()
+    2
+    """
+
+    def _permute_complexity_right(M, iszerofunc):
+        """Permute columns with complicated elements as
+        far right as they can go.  Since the ``sympy`` row reduction
+        algorithms start on the left, having complexity right-shifted
+        speeds things up.
+
+        Returns a tuple (mat, perm) where perm is a permutation
+        of the columns to perform to shift the complex columns right, and mat
+        is the permuted matrix."""
+
+        def complexity(i):
+            # the complexity of a column will be judged by how many
+            # element's zero-ness cannot be determined
+            return sum(1 if iszerofunc(e) is None else 0 for e in M[:, i])
+
+        complex = [(complexity(i), i) for i in range(M.cols)]
+        perm    = [j for (i, j) in sorted(complex)]
+
+        return (M.permute(perm, orientation='cols'), perm)
+
+    simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
+
+    # for small matrices, we compute the rank explicitly
+    # if is_zero on elements doesn't answer the question
+    # for small matrices, we fall back to the full routine.
+    if M.rows <= 0 or M.cols <= 0:
+        return 0
+
+    if M.rows <= 1 or M.cols <= 1:
+        zeros = [iszerofunc(x) for x in M]
+
+        if False in zeros:
+            return 1
+
+    if M.rows == 2 and M.cols == 2:
+        zeros = [iszerofunc(x) for x in M]
+
+        if not False in zeros and not None in zeros:
+            return 0
+
+        d = M.det(dotprodsimp=dotprodsimp)
+
+        if iszerofunc(d) and False in zeros:
+            return 1
+        if iszerofunc(d) is False:
+            return 2
+
+    mat, _       = _permute_complexity_right(M, iszerofunc=iszerofunc)
+    _, pivots, _ = _row_reduce(mat, iszerofunc, simpfunc, normalize_last=True,
+            normalize=False, zero_above=False, dotprodsimp=dotprodsimp)
+
+    return len(pivots)
+
+
+def _rref(M, iszerofunc=_iszero, simplify=False, pivots=True,
+        normalize_last=True, dotprodsimp=None):
+    """Return reduced row-echelon form of matrix and indices of pivot vars.
+
+    Parameters
+    ==========
+
+    iszerofunc : Function
+        A function used for detecting whether an element can
+        act as a pivot.  ``lambda x: x.is_zero`` is used by default.
+
+    simplify : Function
+        A function used to simplify elements when looking for a pivot.
+        By default SymPy's ``simplify`` is used.
+
+    pivots : True or False
+        If ``True``, a tuple containing the row-reduced matrix and a tuple
+        of pivot columns is returned.  If ``False`` just the row-reduced
+        matrix is returned.
+
+    normalize_last : True or False
+        If ``True``, no pivots are normalized to `1` until after all
+        entries above and below each pivot are zeroed.  This means the row
+        reduction algorithm is fraction free until the very last step.
+        If ``False``, the naive row reduction procedure is used where
+        each pivot is normalized to be `1` before row operations are
+        used to zero above and below the pivot.
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> from sympy.abc import x
+    >>> m = Matrix([[1, 2], [x, 1 - 1/x]])
+    >>> m.rref()
+    (Matrix([
+    [1, 0],
+    [0, 1]]), (0, 1))
+    >>> rref_matrix, rref_pivots = m.rref()
+    >>> rref_matrix
+    Matrix([
+    [1, 0],
+    [0, 1]])
+    >>> rref_pivots
+    (0, 1)
+
+    Notes
+    =====
+
+    The default value of ``normalize_last=True`` can provide significant
+    speedup to row reduction, especially on matrices with symbols.  However,
+    if you depend on the form row reduction algorithm leaves entries
+    of the matrix, set ``noramlize_last=False``
+    """
+
+    simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
+
+    mat, pivot_cols, _ = _row_reduce(M, iszerofunc, simpfunc,
+            normalize_last, normalize=True, zero_above=True,
+            dotprodsimp=dotprodsimp)
+
+    if pivots:
+        mat = (mat, pivot_cols)
+
+    return mat
diff --git a/sympy/matrices/solvers.py b/sympy/matrices/solvers.py
index 2299def01f..fa45573c4f 100644
--- a/sympy/matrices/solvers.py
+++ b/sympy/matrices/solvers.py
@@ -1,11 +1,11 @@
 from __future__ import division, print_function
 
 from sympy.core.function import expand_mul
-from sympy.core.symbol import Dummy, _uniquely_named_symbol, symbols
+from sympy.core.symbol import _uniquely_named_symbol
 from sympy.simplify.simplify import dotprodsimp as _dotprodsimp
 from sympy.utilities.iterables import numbered_symbols
 
-from .common import ShapeError, NonSquareMatrixError, NonInvertibleMatrixError
+from .common import ShapeError, NonSquareMatrixError
 from .utilities import _iszero
 
 
@@ -263,25 +263,21 @@ def _cholesky_solve(M, rhs, dotprodsimp=None):
     pinv_solve
     """
 
-    if M.rows < M.cols:
-        raise NotImplementedError(
-            'Under-determined System. Try M.gauss_jordan_solve(rhs)')
-
     hermitian = True
-    reform    = False
 
     if M.is_symmetric():
         hermitian = False
-    elif not M.is_hermitian:
-        reform = True
-
-    if reform or M.is_positive_definite is False:
-        H         = M.H
-        M         = H.multiply(M, dotprodsimp=dotprodsimp)
-        rhs       = H.multiply(rhs, dotprodsimp=dotprodsimp)
-        hermitian = not M.is_symmetric()
+        L = M.cholesky(hermitian=hermitian, dotprodsimp=dotprodsimp)
+    elif M.is_hermitian:
+        L = M.cholesky(hermitian=hermitian, dotprodsimp=dotprodsimp)
+    elif M.rows >= M.cols:
+        L = M.H.multiply(M, dotprodsimp=dotprodsimp).cholesky(
+                hermitian=hermitian, dotprodsimp=dotprodsimp)
+        rhs = M.H.multiply(rhs, dotprodsimp=dotprodsimp)
+    else:
+        raise NotImplementedError('Under-determined System. '
+                                    'Try M.gauss_jordan_solve(rhs)')
 
-    L = M.cholesky(hermitian=hermitian, dotprodsimp=dotprodsimp)
     Y = L.lower_triangular_solve(rhs, dotprodsimp=dotprodsimp)
 
     if hermitian:
@@ -328,27 +324,23 @@ def _LDLsolve(M, rhs, dotprodsimp=None):
     pinv_solve
     """
 
-    if M.rows < M.cols:
-        raise NotImplementedError(
-            'Under-determined System. Try M.gauss_jordan_solve(rhs)')
-
     hermitian = True
-    reform    = False
 
     if M.is_symmetric():
         hermitian = False
-    elif not M.is_hermitian:
-        reform = True
-
-    if reform or M.is_positive_definite is False:
-        H         = M.H
-        M         = H.multiply(M, dotprodsimp=dotprodsimp)
-        rhs       = H.multiply(rhs, dotprodsimp=dotprodsimp)
-        hermitian = not M.is_symmetric()
+        L, D = M.LDLdecomposition(hermitian=hermitian, dotprodsimp=dotprodsimp)
+    elif M.is_hermitian:
+        L, D = M.LDLdecomposition(hermitian=hermitian, dotprodsimp=dotprodsimp)
+    elif M.rows >= M.cols:
+        L, D = M.H.multiply(M, dotprodsimp=dotprodsimp) \
+                .LDLdecomposition(hermitian=hermitian, dotprodsimp=dotprodsimp)
+        rhs = M.H.multiply(rhs, dotprodsimp=dotprodsimp)
+    else:
+        raise NotImplementedError('Under-determined System. '
+                                    'Try M.gauss_jordan_solve(rhs)')
 
-    L, D = M.LDLdecomposition(hermitian=hermitian, dotprodsimp=dotprodsimp)
-    Y    = L.lower_triangular_solve(rhs, dotprodsimp=dotprodsimp)
-    Z    = D.diagonal_solve(Y)
+    Y = L.lower_triangular_solve(rhs, dotprodsimp=dotprodsimp)
+    Z = D.diagonal_solve(Y)
 
     if hermitian:
         return (L.H).upper_triangular_solve(Z, dotprodsimp=dotprodsimp)
@@ -652,265 +644,3 @@ def _gauss_jordan_solve(M, B, freevar=False, dotprodsimp=None):
         return sol, tau, free_var_index
     else:
         return sol, tau
-
-
-def _pinv_solve(M, B, arbitrary_matrix=None, dotprodsimp=None):
-    """Solve ``Ax = B`` using the Moore-Penrose pseudoinverse.
-
-    There may be zero, one, or infinite solutions.  If one solution
-    exists, it will be returned.  If infinite solutions exist, one will
-    be returned based on the value of arbitrary_matrix.  If no solutions
-    exist, the least-squares solution is returned.
-
-    Parameters
-    ==========
-
-    B : Matrix
-        The right hand side of the equation to be solved for.  Must have
-        the same number of rows as matrix A.
-    arbitrary_matrix : Matrix
-        If the system is underdetermined (e.g. A has more columns than
-        rows), infinite solutions are possible, in terms of an arbitrary
-        matrix.  This parameter may be set to a specific matrix to use
-        for that purpose; if so, it must be the same shape as x, with as
-        many rows as matrix A has columns, and as many columns as matrix
-        B.  If left as None, an appropriate matrix containing dummy
-        symbols in the form of ``wn_m`` will be used, with n and m being
-        row and column position of each symbol.
-
-    Returns
-    =======
-
-    x : Matrix
-        The matrix that will satisfy ``Ax = B``.  Will have as many rows as
-        matrix A has columns, and as many columns as matrix B.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
-    >>> B = Matrix([7, 8])
-    >>> A.pinv_solve(B)
-    Matrix([
-    [ _w0_0/6 - _w1_0/3 + _w2_0/6 - 55/18],
-    [-_w0_0/3 + 2*_w1_0/3 - _w2_0/3 + 1/9],
-    [ _w0_0/6 - _w1_0/3 + _w2_0/6 + 59/18]])
-    >>> A.pinv_solve(B, arbitrary_matrix=Matrix([0, 0, 0]))
-    Matrix([
-    [-55/18],
-    [   1/9],
-    [ 59/18]])
-
-    See Also
-    ========
-
-    sympy.matrices.dense.DenseMatrix.lower_triangular_solve
-    sympy.matrices.dense.DenseMatrix.upper_triangular_solve
-    gauss_jordan_solve
-    cholesky_solve
-    diagonal_solve
-    LDLsolve
-    LUsolve
-    QRsolve
-    pinv
-
-    Notes
-    =====
-
-    This may return either exact solutions or least squares solutions.
-    To determine which, check ``A * A.pinv() * B == B``.  It will be
-    True if exact solutions exist, and False if only a least-squares
-    solution exists.  Be aware that the left hand side of that equation
-    may need to be simplified to correctly compare to the right hand
-    side.
-
-    References
-    ==========
-
-    .. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#Obtaining_all_solutions_of_a_linear_system
-
-    """
-
-    from sympy.matrices import eye
-
-    A      = M
-    A_pinv = M.pinv(dotprodsimp=dotprodsimp)
-
-    if arbitrary_matrix is None:
-        rows, cols       = A.cols, B.cols
-        w                = symbols('w:{0}_:{1}'.format(rows, cols), cls=Dummy)
-        arbitrary_matrix = M.__class__(cols, rows, w).T
-
-    return A_pinv.multiply(B, dotprodsimp=dotprodsimp) + (eye(A.cols) -
-            A_pinv.multiply(A, dotprodsimp=dotprodsimp)).multiply(
-            arbitrary_matrix, dotprodsimp=dotprodsimp)
-
-
-def _solve(M, rhs, method='GJ', dotprodsimp=None):
-    """Solves linear equation where the unique solution exists.
-
-    Parameters
-    ==========
-
-    rhs : Matrix
-        Vector representing the right hand side of the linear equation.
-
-    method : string, optional
-        If set to ``'GJ'`` or ``'GE'``, the Gauss-Jordan elimination will be
-        used, which is implemented in the routine ``gauss_jordan_solve``.
-
-        If set to ``'LU'``, ``LUsolve`` routine will be used.
-
-        If set to ``'QR'``, ``QRsolve`` routine will be used.
-
-        If set to ``'PINV'``, ``pinv_solve`` routine will be used.
-
-        It also supports the methods available for special linear systems
-
-        For positive definite systems:
-
-        If set to ``'CH'``, ``cholesky_solve`` routine will be used.
-
-        If set to ``'LDL'``, ``LDLsolve`` routine will be used.
-
-        To use a different method and to compute the solution via the
-        inverse, use a method defined in the .inv() docstring.
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Returns
-    =======
-
-    solutions : Matrix
-        Vector representing the solution.
-
-    Raises
-    ======
-
-    ValueError
-        If there is not a unique solution then a ``ValueError`` will be
-        raised.
-
-        If ``M`` is not square, a ``ValueError`` and a different routine
-        for solving the system will be suggested.
-    """
-
-    if method == 'GJ' or method == 'GE':
-        try:
-            soln, param = M.gauss_jordan_solve(rhs, dotprodsimp=dotprodsimp)
-
-            if param:
-                raise NonInvertibleMatrixError("Matrix det == 0; not invertible. "
-                "Try ``M.gauss_jordan_solve(rhs)`` to obtain a parametric solution.")
-
-        except ValueError:
-            raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
-
-        return soln
-
-    elif method == 'LU':
-        return M.LUsolve(rhs, dotprodsimp=dotprodsimp)
-    elif method == 'CH':
-        return M.cholesky_solve(rhs, dotprodsimp=dotprodsimp)
-    elif method == 'QR':
-        return M.QRsolve(rhs, dotprodsimp=dotprodsimp)
-    elif method == 'LDL':
-        return M.LDLsolve(rhs, dotprodsimp=dotprodsimp)
-    elif method == 'PINV':
-        return M.pinv_solve(rhs, dotprodsimp=dotprodsimp)
-    else:
-        return M.inv(method=method, dotprodsimp=dotprodsimp) \
-                .multiply(rhs, dotprodsimp=dotprodsimp)
-
-
-def _solve_least_squares(M, rhs, method='CH'):
-    """Return the least-square fit to the data.
-
-    Parameters
-    ==========
-
-    rhs : Matrix
-        Vector representing the right hand side of the linear equation.
-
-    method : string or boolean, optional
-        If set to ``'CH'``, ``cholesky_solve`` routine will be used.
-
-        If set to ``'LDL'``, ``LDLsolve`` routine will be used.
-
-        If set to ``'QR'``, ``QRsolve`` routine will be used.
-
-        If set to ``'PINV'``, ``pinv_solve`` routine will be used.
-
-        Otherwise, the conjugate of ``M`` will be used to create a system
-        of equations that is passed to ``solve`` along with the hint
-        defined by ``method``.
-
-    Returns
-    =======
-
-    solutions : Matrix
-        Vector representing the solution.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix, ones
-    >>> A = Matrix([1, 2, 3])
-    >>> B = Matrix([2, 3, 4])
-    >>> S = Matrix(A.row_join(B))
-    >>> S
-    Matrix([
-    [1, 2],
-    [2, 3],
-    [3, 4]])
-
-    If each line of S represent coefficients of Ax + By
-    and x and y are [2, 3] then S*xy is:
-
-    >>> r = S*Matrix([2, 3]); r
-    Matrix([
-    [ 8],
-    [13],
-    [18]])
-
-    But let's add 1 to the middle value and then solve for the
-    least-squares value of xy:
-
-    >>> xy = S.solve_least_squares(Matrix([8, 14, 18])); xy
-    Matrix([
-    [ 5/3],
-    [10/3]])
-
-    The error is given by S*xy - r:
-
-    >>> S*xy - r
-    Matrix([
-    [1/3],
-    [1/3],
-    [1/3]])
-    >>> _.norm().n(2)
-    0.58
-
-    If a different xy is used, the norm will be higher:
-
-    >>> xy += ones(2, 1)/10
-    >>> (S*xy - r).norm().n(2)
-    1.5
-
-    """
-
-    if method == 'CH':
-        return M.cholesky_solve(rhs)
-    elif method == 'QR':
-        return M.QRsolve(rhs)
-    elif method == 'LDL':
-        return M.LDLsolve(rhs)
-    elif method == 'PINV':
-        return M.pinv_solve(rhs)
-    else:
-        t = M.H
-        return (t * M).solve(t * rhs, method=method)
diff --git a/sympy/matrices/sparse.py b/sympy/matrices/sparse.py
index f42e8783be..b72a3bfcc5 100644
--- a/sympy/matrices/sparse.py
+++ b/sympy/matrices/sparse.py
@@ -15,8 +15,6 @@
 from .dense import Matrix
 from .matrices import MatrixBase, ShapeError
 
-from .utilities import _iszero
-
 from .decompositions import (
     _liupc, _row_structure_symbolic_cholesky, _cholesky_sparse,
     _LDLdecomposition_sparse)
@@ -276,11 +274,57 @@ def __getitem__(self, key):
     def __setitem__(self, key, value):
         raise NotImplementedError()
 
-    def _eval_inverse(self, **kwargs):
-        return self.inv(method=kwargs.get('method', 'LDL'),
-                        iszerofunc=kwargs.get('iszerofunc', _iszero),
-                        dotprodsimp=kwargs.get('dotprodsimp', None),
-                        try_block_diag=kwargs.get('try_block_diag', False))
+    def _eval_inverse(self, method='LDL', dotprodsimp=None):
+        """Return the matrix inverse using Cholesky or LDL (default)
+        decomposition as selected with the ``method`` keyword: 'CH' or 'LDL',
+        respectively.
+
+        Parameters
+        ==========
+
+        dotprodsimp : bool, optional
+            Specifies whether intermediate term algebraic simplification is used
+            during matrix multiplications to control expression blowup and thus
+            speed up calculation.
+
+        Examples
+        ========
+
+        >>> from sympy import SparseMatrix, Matrix
+        >>> A = SparseMatrix([
+        ... [ 2, -1,  0],
+        ... [-1,  2, -1],
+        ... [ 0,  0,  2]])
+        >>> A.inv('CH')
+        Matrix([
+        [2/3, 1/3, 1/6],
+        [1/3, 2/3, 1/3],
+        [  0,   0, 1/2]])
+        >>> A.inv(method='LDL') # use of 'method=' is optional
+        Matrix([
+        [2/3, 1/3, 1/6],
+        [1/3, 2/3, 1/3],
+        [  0,   0, 1/2]])
+        >>> A * _
+        Matrix([
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1]])
+        """
+        M = self.as_mutable()
+        I = M.eye(M.rows)
+        if not self.is_symmetric():
+            I = M.T
+            M = I.multiply(M, dotprodsimp=dotprodsimp)
+
+        if method == "LDL":
+            rv = M.LDLsolve(I, dotprodsimp=dotprodsimp)
+        elif method == "CH":
+            rv = M.cholesky_solve(I, dotprodsimp=dotprodsimp)
+        else:
+            raise NotImplementedError(
+                'Method may be "CH" or "LDL", not %s.' % method)
+        return self._new(rv)
 
     def _eval_Abs(self):
         return self.applyfunc(lambda x: Abs(x))
diff --git a/sympy/matrices/subspaces.py b/sympy/matrices/subspaces.py
index fbe9ee6dc6..e78cac9d94 100644
--- a/sympy/matrices/subspaces.py
+++ b/sympy/matrices/subspaces.py
@@ -1,222 +1,222 @@
-from __future__ import division, print_function
-
-from sympy.core.compatibility import reduce
-
-from .utilities import _iszero
-
-
-def _columnspace(M, simplify=False, dotprodsimp=None):
-    """Returns a list of vectors (Matrix objects) that span columnspace of ``M``
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
-    >>> M
-    Matrix([
-    [ 1,  3, 0],
-    [-2, -6, 0],
-    [ 3,  9, 6]])
-    >>> M.columnspace()
-    [Matrix([
-    [ 1],
-    [-2],
-    [ 3]]), Matrix([
-    [0],
-    [0],
-    [6]])]
-
-    See Also
-    ========
-
-    nullspace
-    rowspace
-    """
-
-    reduced, pivots = M.echelon_form(simplify=simplify, with_pivots=True,
-            dotprodsimp=dotprodsimp)
-
-    return [M.col(i) for i in pivots]
-
-
-def _nullspace(M, simplify=False, iszerofunc=_iszero, dotprodsimp=None):
-    """Returns list of vectors (Matrix objects) that span nullspace of ``M``
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy.matrices import Matrix
-    >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
-    >>> M
-    Matrix([
-    [ 1,  3, 0],
-    [-2, -6, 0],
-    [ 3,  9, 6]])
-    >>> M.nullspace()
-    [Matrix([
-    [-3],
-    [ 1],
-    [ 0]])]
-
-    See Also
-    ========
-
-    columnspace
-    rowspace
-    """
-
-    reduced, pivots = M.rref(iszerofunc=iszerofunc, simplify=simplify,
-            dotprodsimp=dotprodsimp)
-
-    free_vars = [i for i in range(M.cols) if i not in pivots]
-    basis     = []
-
-    for free_var in free_vars:
-        # for each free variable, we will set it to 1 and all others
-        # to 0.  Then, we will use back substitution to solve the system
-        vec           = [M.zero] * M.cols
-        vec[free_var] = M.one
-
-        for piv_row, piv_col in enumerate(pivots):
-            vec[piv_col] -= reduced[piv_row, free_var]
-
-        basis.append(vec)
-
-    return [M._new(M.cols, 1, b) for b in basis]
-
-
-def _rowspace(M, simplify=False, dotprodsimp=None):
-    """Returns a list of vectors that span the row space of ``M``.
-
-    Parameters
-    ==========
-
-    dotprodsimp : bool, optional
-        Specifies whether intermediate term algebraic simplification is used
-        during matrix multiplications to control expression blowup and thus
-        speed up calculation.
-
-    Examples
-    ========
-
-    >>> from sympy import Matrix
-    >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
-    >>> M
-    Matrix([
-    [ 1,  3, 0],
-    [-2, -6, 0],
-    [ 3,  9, 6]])
-    >>> M.rowspace()
-    [Matrix([[1, 3, 0]]), Matrix([[0, 0, 6]])]
-    """
-
-    reduced, pivots = M.echelon_form(simplify=simplify, with_pivots=True,
-            dotprodsimp=dotprodsimp)
-
-    return [reduced.row(i) for i in range(len(pivots))]
-
-
-def _orthogonalize(cls, *vecs, **kwargs):
-    """Apply the Gram-Schmidt orthogonalization procedure
-    to vectors supplied in ``vecs``.
-
-    Parameters
-    ==========
-
-    vecs
-        vectors to be made orthogonal
-
-    normalize : bool
-        If ``True``, return an orthonormal basis.
-
-    rankcheck : bool
-        If ``True``, the computation does not stop when encountering
-        linearly dependent vectors.
-
-        If ``False``, it will raise ``ValueError`` when any zero
-        or linearly dependent vectors are found.
-
-    Returns
-    =======
-
-    list
-        List of orthogonal (or orthonormal) basis vectors.
-
-    Examples
-    ========
-
-    >>> from sympy import I, Matrix
-    >>> v = [Matrix([1, I]), Matrix([1, -I])]
-    >>> Matrix.orthogonalize(*v)
-    [Matrix([
-    [1],
-    [I]]), Matrix([
-    [ 1],
-    [-I]])]
-
-    See Also
-    ========
-
-    MatrixBase.QRdecomposition
-
-    References
-    ==========
-
-    .. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
-    """
-
-    normalize = kwargs.get('normalize', False)
-    rankcheck = kwargs.get('rankcheck', False)
-
-    def project(a, b):
-        return b * (a.dot(b, hermitian=True) / b.dot(b, hermitian=True))
-
-    def perp_to_subspace(vec, basis):
-        """projects vec onto the subspace given
-        by the orthogonal basis ``basis``"""
-
-        components = [project(vec, b) for b in basis]
-
-        if len(basis) == 0:
-            return vec
-
-        return vec - reduce(lambda a, b: a + b, components)
-
-    ret  = []
-    vecs = list(vecs) # make sure we start with a non-zero vector
-
-    while len(vecs) > 0 and vecs[0].is_zero_matrix:
-        if rankcheck is False:
-            del vecs[0]
-        else:
-            raise ValueError("GramSchmidt: vector set not linearly independent")
-
-    for vec in vecs:
-        perp = perp_to_subspace(vec, ret)
-
-        if not perp.is_zero_matrix:
-            ret.append(cls(perp))
-        elif rankcheck is True:
-            raise ValueError("GramSchmidt: vector set not linearly independent")
-
-    if normalize:
-        ret = [vec / vec.norm() for vec in ret]
-
-    return ret
+from __future__ import division, print_function
+
+from sympy.core.compatibility import reduce
+
+from .utilities import _iszero
+
+
+def _columnspace(M, simplify=False, dotprodsimp=None):
+    """Returns a list of vectors (Matrix objects) that span columnspace of ``M``
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
+    >>> M
+    Matrix([
+    [ 1,  3, 0],
+    [-2, -6, 0],
+    [ 3,  9, 6]])
+    >>> M.columnspace()
+    [Matrix([
+    [ 1],
+    [-2],
+    [ 3]]), Matrix([
+    [0],
+    [0],
+    [6]])]
+
+    See Also
+    ========
+
+    nullspace
+    rowspace
+    """
+
+    reduced, pivots = M.echelon_form(simplify=simplify, with_pivots=True,
+            dotprodsimp=dotprodsimp)
+
+    return [M.col(i) for i in pivots]
+
+
+def _nullspace(M, simplify=False, iszerofunc=_iszero, dotprodsimp=None):
+    """Returns list of vectors (Matrix objects) that span nullspace of ``M``
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy.matrices import Matrix
+    >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
+    >>> M
+    Matrix([
+    [ 1,  3, 0],
+    [-2, -6, 0],
+    [ 3,  9, 6]])
+    >>> M.nullspace()
+    [Matrix([
+    [-3],
+    [ 1],
+    [ 0]])]
+
+    See Also
+    ========
+
+    columnspace
+    rowspace
+    """
+
+    reduced, pivots = M.rref(iszerofunc=iszerofunc, simplify=simplify,
+            dotprodsimp=dotprodsimp)
+
+    free_vars = [i for i in range(M.cols) if i not in pivots]
+    basis     = []
+
+    for free_var in free_vars:
+        # for each free variable, we will set it to 1 and all others
+        # to 0.  Then, we will use back substitution to solve the system
+        vec           = [M.zero] * M.cols
+        vec[free_var] = M.one
+
+        for piv_row, piv_col in enumerate(pivots):
+            vec[piv_col] -= reduced[piv_row, free_var]
+
+        basis.append(vec)
+
+    return [M._new(M.cols, 1, b) for b in basis]
+
+
+def _rowspace(M, simplify=False, dotprodsimp=None):
+    """Returns a list of vectors that span the row space of ``M``.
+
+    Parameters
+    ==========
+
+    dotprodsimp : bool, optional
+        Specifies whether intermediate term algebraic simplification is used
+        during matrix multiplications to control expression blowup and thus
+        speed up calculation.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
+    >>> M
+    Matrix([
+    [ 1,  3, 0],
+    [-2, -6, 0],
+    [ 3,  9, 6]])
+    >>> M.rowspace()
+    [Matrix([[1, 3, 0]]), Matrix([[0, 0, 6]])]
+    """
+
+    reduced, pivots = M.echelon_form(simplify=simplify, with_pivots=True,
+            dotprodsimp=dotprodsimp)
+
+    return [reduced.row(i) for i in range(len(pivots))]
+
+
+def _orthogonalize(cls, *vecs, **kwargs):
+    """Apply the Gram-Schmidt orthogonalization procedure
+    to vectors supplied in ``vecs``.
+
+    Parameters
+    ==========
+
+    vecs
+        vectors to be made orthogonal
+
+    normalize : bool
+        If ``True``, return an orthonormal basis.
+
+    rankcheck : bool
+        If ``True``, the computation does not stop when encountering
+        linearly dependent vectors.
+
+        If ``False``, it will raise ``ValueError`` when any zero
+        or linearly dependent vectors are found.
+
+    Returns
+    =======
+
+    list
+        List of orthogonal (or orthonormal) basis vectors.
+
+    Examples
+    ========
+
+    >>> from sympy import I, Matrix
+    >>> v = [Matrix([1, I]), Matrix([1, -I])]
+    >>> Matrix.orthogonalize(*v)
+    [Matrix([
+    [1],
+    [I]]), Matrix([
+    [ 1],
+    [-I]])]
+
+    See Also
+    ========
+
+    MatrixBase.QRdecomposition
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
+    """
+
+    normalize = kwargs.get('normalize', False)
+    rankcheck = kwargs.get('rankcheck', False)
+
+    def project(a, b):
+        return b * (a.dot(b, hermitian=True) / b.dot(b, hermitian=True))
+
+    def perp_to_subspace(vec, basis):
+        """projects vec onto the subspace given
+        by the orthogonal basis ``basis``"""
+
+        components = [project(vec, b) for b in basis]
+
+        if len(basis) == 0:
+            return vec
+
+        return vec - reduce(lambda a, b: a + b, components)
+
+    ret  = []
+    vecs = list(vecs) # make sure we start with a non-zero vector
+
+    while len(vecs) > 0 and vecs[0].is_zero_matrix:
+        if rankcheck is False:
+            del vecs[0]
+        else:
+            raise ValueError("GramSchmidt: vector set not linearly independent")
+
+    for vec in vecs:
+        perp = perp_to_subspace(vec, ret)
+
+        if not perp.is_zero_matrix:
+            ret.append(cls(perp))
+        elif rankcheck is True:
+            raise ValueError("GramSchmidt: vector set not linearly independent")
+
+    if normalize:
+        ret = [vec / vec.norm() for vec in ret]
+
+    return ret
diff --git a/sympy/matrices/utilities.py b/sympy/matrices/utilities.py
index 01e7925b95..6024028318 100644
--- a/sympy/matrices/utilities.py
+++ b/sympy/matrices/utilities.py
@@ -1,14 +1,14 @@
-from __future__ import division, print_function
-
-from sympy.core.function import expand_mul
-
-
-def _iszero(x):
-    """Returns True if x is zero."""
-    return getattr(x, 'is_zero', None)
-
-
-def _is_zero_after_expand_mul(x):
-    """Tests by expand_mul only, suitable for polynomials and rational
-    functions."""
-    return expand_mul(x) == 0
+from __future__ import division, print_function
+
+from sympy.core.function import expand_mul
+
+
+def _iszero(x):
+    """Returns True if x is zero."""
+    return getattr(x, 'is_zero', None)
+
+
+def _is_zero_after_expand_mul(x):
+    """Tests by expand_mul only, suitable for polynomials and rational
+    functions."""
+    return expand_mul(x) == 0
diff --git a/sympy/solvers/recurr.py b/sympy/solvers/recurr.py
index 37c852d661..4af9ba4a64 100644
--- a/sympy/solvers/recurr.py
+++ b/sympy/solvers/recurr.py
@@ -813,10 +813,10 @@ def rsolve(f, y, init=None):
                     i = int(k.args[0])
                 else:
                     raise ValueError("Integer or term expected, got '%s'" % k)
-
-            eq = solution.subs(n, i) - v
-            if eq.has(S.NaN):
+            try:
                 eq = solution.limit(n, i) - v
+            except NotImplementedError:
+                eq = solution.subs(n, i) - v
             equations.append(eq)
 
         result = solve(equations, *symbols)


Test patch:
diff --git a/sympy/core/tests/test_function.py b/sympy/core/tests/test_function.py
index 97bb7d3e11..6ed938d4b2 100644
--- a/sympy/core/tests/test_function.py
+++ b/sympy/core/tests/test_function.py
@@ -1350,3 +1350,8 @@ def test_Derivative_free_symbols():
     f = Function('f')
     n = Symbol('n', integer=True, positive=True)
     assert diff(f(x), (x, n)).free_symbols == {n, x}
+
+
+def test_issue_10503():
+    f = exp(x**3)*cos(x**6)
+    assert f.series(x, 0, 14) == 1 + x**3 + x**6/2 + x**9/6 - 11*x**12/24 + O(x**14)
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index ffc5e896de..ec579006f1 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -1116,9 +1116,6 @@ def test_inverse():
     assert A.inv() == eye(4)
     assert A.inv(method="LU") == eye(4)
     assert A.inv(method="ADJ") == eye(4)
-    assert A.inv(method="CH") == eye(4)
-    assert A.inv(method="LDL") == eye(4)
-    assert A.inv(method="QR") == eye(4)
     A = Matrix([[2, 3, 5],
                 [3, 6, 2],
                 [8, 3, 6]])
@@ -1126,21 +1123,18 @@ def test_inverse():
     assert A*Ainv == eye(3)
     assert A.inv(method="LU") == Ainv
     assert A.inv(method="ADJ") == Ainv
-    assert A.inv(method="CH") == Ainv
-    assert A.inv(method="LDL") == Ainv
-    assert A.inv(method="QR") == Ainv
 
     # test that immutability is not a problem
     cls = ImmutableMatrix
     m = cls([[48, 49, 31],
              [ 9, 71, 94],
              [59, 28, 65]])
-    assert all(type(m.inv(s)) is cls for s in 'GE ADJ LU CH LDL QR'.split())
+    assert all(type(m.inv(s)) is cls for s in 'GE ADJ LU'.split())
     cls = ImmutableSparseMatrix
     m = cls([[48, 49, 31],
              [ 9, 71, 94],
              [59, 28, 65]])
-    assert all(type(m.inv(s)) is cls for s in 'GE ADJ LU CH LDL QR'.split())
+    assert all(type(m.inv(s)) is cls for s in 'CH LDL'.split())
 
 
 def test_matrix_inverse_mod():
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index ec6d61f56f..48887a3db9 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -533,22 +533,6 @@ def test_issue_14411():
     assert limit(3*sec(4*pi*x - x/3), x, 3*pi/(24*pi - 2)) is -oo
 
 
-def test_issue_13382():
-    assert limit(x*(((x + 1)**2 + 1)/(x**2 + 1) - 1), x, oo) == 2
-
-
-def test_issue_13403():
-    assert limit(x*(-1 + (x + log(x + 1) + 1)/(x + log(x))), x ,oo) == 1
-
-
-def test_issue_13416():
-    assert limit((-x**3*log(x)**3 + (x - 1)*(x + 1)**2*log(x + 1)**3)/(x**2*log(x)**3), x ,oo) == 1
-
-
-def test_issue_13462():
-    assert limit(n**2*(2*n*(-(1 - 1/(2*n))**x + 1) - x - (-x**2/4 + x/4)/n), n, oo) == x*(x**2 - 3*x + 2)/24
-
-
 def test_issue_14574():
     assert limit(sqrt(x)*cos(x - x**2) / (x + 1), x, oo) == 0