text: import mpmath as mp
  from collections.abc import Callable
  
  
  from sympy.core.add import Add
  from sympy.core.basic import Basic
  from sympy.core.function import diff
  from sympy.core.expr import Expr
  from sympy.core.kind import _NumberKind, UndefinedKind
  from sympy.core.mul import Mul
  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.sympify import sympify, _sympify
  from sympy.functions.combinatorial.factorials import binomial, factorial
  from sympy.functions.elementary.complexes import re
  from sympy.functions.elementary.exponential import exp, log
  from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
  from sympy.functions.special.tensor_functions import KroneckerDelta, LeviCivita
  from sympy.polys import cancel
  from sympy.printing import sstr
  from sympy.printing.defaults import Printable
  from sympy.printing.str import StrPrinter
  from sympy.utilities.decorator import deprecated
  from sympy.utilities.exceptions import sympy_deprecation_warning
  from sympy.utilities.iterables import flatten, NotIterable, is_sequence, reshape
  from sympy.utilities.misc import as_int, filldedent
  
  from .common import (
      MatrixCommon, MatrixError, NonSquareMatrixError, NonInvertibleMatrixError,
      ShapeError, MatrixKind)
  
  from .utilities import _iszero, _is_zero_after_expand_mul, _simplify
  
  from .determinant import (
      _find_reasonable_pivot, _find_reasonable_pivot_naive,
      _adjugate, _charpoly, _cofactor, _cofactor_matrix, _per,
      _det, _det_bareiss, _det_berkowitz, _det_LU, _minor, _minor_submatrix)
  
  from .reductions import _is_echelon, _echelon_form, _rank, _rref
  from .subspaces import _columnspace, _nullspace, _rowspace, _orthogonalize
  
  from .eigen import (
      _eigenvals, _eigenvects,
      _bidiagonalize, _bidiagonal_decomposition,
      _is_diagonalizable, _diagonalize,
      _is_positive_definite, _is_positive_semidefinite,
      _is_negative_definite, _is_negative_semidefinite, _is_indefinite,
      _jordan_form, _left_eigenvects, _singular_values)
  
  from .decompositions import (
      _rank_decomposition, _cholesky, _LDLdecomposition,
      _LUdecomposition, _LUdecomposition_Simple, _LUdecompositionFF,
      _singular_value_decomposition, _QRdecomposition, _upper_hessenberg_decomposition)
  
  from .graph import (
      _connected_components, _connected_components_decomposition,
      _strongly_connected_components, _strongly_connected_components_decomposition)
  
  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, _inv_block)
  
  
  class DeferredVector(Symbol, NotIterable):
      """A vector whose components are deferred (e.g. for use with lambdify)
  
      Examples
      ========
  
      >>> from sympy import DeferredVector, lambdify
      >>> X = DeferredVector( 'X' )
      >>> X
      X
      >>> expr = (X[0] + 2, X[2] + 3)
      >>> func = lambdify( X, expr)
      >>> func( [1, 2, 3] )
      (3, 6)
      """
  
      def __getitem__(self, i):
          if i == -0:
              i = 0
          if i < 0:
              raise IndexError('DeferredVector index out of range')
          component_name = '%s[%d]' % (self.name, i)
          return Symbol(component_name)
  
      def __str__(self):
          return sstr(self)
  
      def __repr__(self):
          return "DeferredVector('%s')" % self.name
  
  
  class MatrixDeterminant(MatrixCommon):
      """Provides basic matrix determinant operations. Should not be instantiated
      directly. See ``determinant.py`` for their implementations."""
  
      def _eval_det_bareiss(self, iszerofunc=_is_zero_after_expand_mul):
          return _det_bareiss(self, iszerofunc=iszerofunc)
  
      def _eval_det_berkowitz(self):
          return _det_berkowitz(self)
  
      def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None):
          return _det_LU(self, iszerofunc=iszerofunc, simpfunc=simpfunc)
  
      def _eval_determinant(self): # for expressions.determinant.Determinant
          return _det(self)
  
      def adjugate(self, method="berkowitz"):
          return _adjugate(self, method=method)
  
      def charpoly(self, x='lambda', simplify=_simplify):
          return _charpoly(self, x=x, simplify=simplify)
  
      def cofactor(self, i, j, method="berkowitz"):
          return _cofactor(self, i, j, method=method)
  
      def cofactor_matrix(self, method="berkowitz"):
          return _cofactor_matrix(self, method=method)
  
      def det(self, method="bareiss", iszerofunc=None):
          return _det(self, method=method, iszerofunc=iszerofunc)
  
      def per(self):
          return _per(self)
  
      def minor(self, i, j, method="berkowitz"):
          return _minor(self, i, j, method=method)
  
      def minor_submatrix(self, i, j):
          return _minor_submatrix(self, i, j)
  
      _find_reasonable_pivot.__doc__       = _find_reasonable_pivot.__doc__
      _find_reasonable_pivot_naive.__doc__ = _find_reasonable_pivot_naive.__doc__
      _eval_det_bareiss.__doc__            = _det_bareiss.__doc__
      _eval_det_berkowitz.__doc__          = _det_berkowitz.__doc__
      _eval_det_lu.__doc__                 = _det_LU.__doc__
      _eval_determinant.__doc__            = _det.__doc__
      adjugate.__doc__                     = _adjugate.__doc__
      charpoly.__doc__                     = _charpoly.__doc__
      cofactor.__doc__                     = _cofactor.__doc__
      cofactor_matrix.__doc__              = _cofactor_matrix.__doc__
      det.__doc__                          = _det.__doc__
      per.__doc__                          = _per.__doc__
      minor.__doc__                        = _minor.__doc__
      minor_submatrix.__doc__              = _minor_submatrix.__doc__
  
  
  class MatrixReductions(MatrixDeterminant):
      """Provides basic matrix row/column operations. Should not be instantiated
      directly. See ``reductions.py`` for some of their implementations."""
  
      def echelon_form(self, iszerofunc=_iszero, simplify=False, with_pivots=False):
          return _echelon_form(self, iszerofunc=iszerofunc, simplify=simplify,
                  with_pivots=with_pivots)
  
      @property
      def is_echelon(self):
          return _is_echelon(self)
  
      def rank(self, iszerofunc=_iszero, simplify=False):
          return _rank(self, iszerofunc=iszerofunc, simplify=simplify)
  
      def rref(self, iszerofunc=_iszero, simplify=False, pivots=True,
              normalize_last=True):
          return _rref(self, iszerofunc=iszerofunc, simplify=simplify,
              pivots=pivots, normalize_last=normalize_last)
  
      echelon_form.__doc__ = _echelon_form.__doc__
      is_echelon.__doc__   = _is_echelon.__doc__
      rank.__doc__         = _rank.__doc__
      rref.__doc__         = _rref.__doc__
  
      def _normalize_op_args(self, op, col, k, col1, col2, error_str="col"):
          """Validate the arguments for a row/column operation.  ``error_str``
          can be one of "row" or "col" depending on the arguments being parsed."""
          if op not in ["n->kn", "n<->m", "n->n+km"]:
              raise ValueError("Unknown {} operation '{}'. Valid col operations "
                               "are 'n->kn', 'n<->m', 'n->n+km'".format(error_str, op))
  
          # define self_col according to error_str
          self_cols = self.cols if error_str == 'col' else self.rows
  
          # normalize and validate the arguments
          if op == "n->kn":
              col = col if col is not None else col1
              if col is None or k is None:
                  raise ValueError("For a {0} operation 'n->kn' you must provide the "
                                   "kwargs `{0}` and `k`".format(error_str))
              if not 0 <= col < self_cols:
                  raise ValueError("This matrix does not have a {} '{}'".format(error_str, col))
  
          elif op == "n<->m":
              # we need two cols to swap. It does not matter
              # how they were specified, so gather them together and
              # remove `None`
              cols = {col, k, col1, col2}.difference([None])
              if len(cols) > 2:
                  # maybe the user left `k` by mistake?
                  cols = {col, col1, col2}.difference([None])
              if len(cols) != 2:
                  raise ValueError("For a {0} operation 'n<->m' you must provide the "
                                   "kwargs `{0}1` and `{0}2`".format(error_str))
              col1, col2 = cols
              if not 0 <= col1 < self_cols:
                  raise ValueError("This matrix does not have a {} '{}'".format(error_str, col1))
              if not 0 <= col2 < self_cols:
                  raise ValueError("This matrix does not have a {} '{}'".format(error_str, col2))
  
          elif op == "n->n+km":
              col = col1 if col is None else col
              col2 = col1 if col2 is None else col2
              if col is None or col2 is None or k is None:
                  raise ValueError("For a {0} operation 'n->n+km' you must provide the "
                                   "kwargs `{0}`, `k`, and `{0}2`".format(error_str))
              if col == col2:
                  raise ValueError("For a {0} operation 'n->n+km' `{0}` and `{0}2` must "
                                   "be different.".format(error_str))
              if not 0 <= col < self_cols:
                  raise ValueError("This matrix does not have a {} '{}'".format(error_str, col))
              if not 0 <= col2 < self_cols:
                  raise ValueError("This matrix does not have a {} '{}'".format(error_str, col2))
  
          else:
              raise ValueError('invalid operation %s' % repr(op))
  
          return op, col, k, col1, col2
  
      def _eval_col_op_multiply_col_by_const(self, col, k):
          def entry(i, j):
              if j == col:
                  return k * self[i, j]
              return self[i, j]
          return self._new(self.rows, self.cols, entry)
  
      def _eval_col_op_swap(self, col1, col2):
          def entry(i, j):
              if j == col1:
                  return self[i, col2]
              elif j == col2:
                  return self[i, col1]
              return self[i, j]
          return self._new(self.rows, self.cols, entry)
  
      def _eval_col_op_add_multiple_to_other_col(self, col, k, col2):
          def entry(i, j):
              if j == col:
                  return self[i, j] + k * self[i, col2]
              return self[i, j]
          return self._new(self.rows, self.cols, entry)
  
      def _eval_row_op_swap(self, row1, row2):
          def entry(i, j):
              if i == row1:
                  return self[row2, j]
              elif i == row2:
                  return self[row1, j]
              return self[i, j]
          return self._new(self.rows, self.cols, entry)
  
      def _eval_row_op_multiply_row_by_const(self, row, k):
          def entry(i, j):
              if i == row:
                  return k * self[i, j]
              return self[i, j]
          return self._new(self.rows, self.cols, entry)
  
      def _eval_row_op_add_multiple_to_other_row(self, row, k, row2):
          def entry(i, j):
              if i == row:
                  return self[i, j] + k * self[row2, j]
              return self[i, j]
          return self._new(self.rows, self.cols, entry)
  
      def elementary_col_op(self, op="n->kn", col=None, k=None, col1=None, col2=None):
          """Performs the elementary column operation `op`.
  
          `op` may be one of
  
              * ``"n->kn"`` (column n goes to k*n)
              * ``"n<->m"`` (swap column n and column m)
              * ``"n->n+km"`` (column n goes to column n + k*column m)
  
          Parameters
          ==========
  
          op : string; the elementary row operation
          col : the column to apply the column operation
          k : the multiple to apply in the column operation
          col1 : one column of a column swap
          col2 : second column of a column swap or column "m" in the column operation
                 "n->n+km"
          """
  
          op, col, k, col1, col2 = self._normalize_op_args(op, col, k, col1, col2, "col")
  
          # now that we've validated, we're all good to dispatch
          if op == "n->kn":
              return self._eval_col_op_multiply_col_by_const(col, k)
          if op == "n<->m":
              return self._eval_col_op_swap(col1, col2)
          if op == "n->n+km":
              return self._eval_col_op_add_multiple_to_other_col(col, k, col2)
  
      def elementary_row_op(self, op="n->kn", row=None, k=None, row1=None, row2=None):
          """Performs the elementary row operation `op`.
  
          `op` may be one of
  
              * ``"n->kn"`` (row n goes to k*n)
              * ``"n<->m"`` (swap row n and row m)
              * ``"n->n+km"`` (row n goes to row n + k*row m)
  
          Parameters
          ==========
  
          op : string; the elementary row operation
          row : the row to apply the row operation
          k : the multiple to apply in the row operation
          row1 : one row of a row swap
          row2 : second row of a row swap or row "m" in the row operation
                 "n->n+km"
          """
  
          op, row, k, row1, row2 = self._normalize_op_args(op, row, k, row1, row2, "row")
  
          # now that we've validated, we're all good to dispatch
          if op == "n->kn":
              return self._eval_row_op_multiply_row_by_const(row, k)
          if op == "n<->m":
              return self._eval_row_op_swap(row1, row2)
          if op == "n->n+km":
              return self._eval_row_op_add_multiple_to_other_row(row, k, row2)
  
  
  class MatrixSubspaces(MatrixReductions):
      """Provides methods relating to the fundamental subspaces of a matrix.
      Should not be instantiated directly. See ``subspaces.py`` for their
      implementations."""
  
      def columnspace(self, simplify=False):
          return _columnspace(self, simplify=simplify)
  
      def nullspace(self, simplify=False, iszerofunc=_iszero):
          return _nullspace(self, simplify=simplify, iszerofunc=iszerofunc)
  
      def rowspace(self, simplify=False):
          return _rowspace(self, simplify=simplify)
  
      # This is a classmethod but is converted to such later in order to allow
      # assignment of __doc__ since that does not work for already wrapped
      # classmethods in Python 3.6.
      def orthogonalize(cls, *vecs, **kwargs):
          return _orthogonalize(cls, *vecs, **kwargs)
  
      columnspace.__doc__   = _columnspace.__doc__
      nullspace.__doc__     = _nullspace.__doc__
      rowspace.__doc__      = _rowspace.__doc__
      orthogonalize.__doc__ = _orthogonalize.__doc__
  
      orthogonalize         = classmethod(orthogonalize)  # type:ignore
  
  
  class MatrixEigen(MatrixSubspaces):
      """Provides basic matrix eigenvalue/vector operations.
      Should not be instantiated directly. See ``eigen.py`` for their
      implementations."""
  
      def eigenvals(self, error_when_incomplete=True, **flags):
          return _eigenvals(self, error_when_incomplete=error_when_incomplete, **flags)
  
      def eigenvects(self, error_when_incomplete=True, iszerofunc=_iszero, **flags):
          return _eigenvects(self, error_when_incomplete=error_when_incomplete,
                  iszerofunc=iszerofunc, **flags)
  
      def is_diagonalizable(self, reals_only=False, **kwargs):
          return _is_diagonalizable(self, reals_only=reals_only, **kwargs)
  
      def diagonalize(self, reals_only=False, sort=False, normalize=False):
          return _diagonalize(self, reals_only=reals_only, sort=sort,
                  normalize=normalize)
  
      def bidiagonalize(self, upper=True):
          return _bidiagonalize(self, upper=upper)
  
      def bidiagonal_decomposition(self, upper=True):
          return _bidiagonal_decomposition(self, upper=upper)
  
      @property
      def is_positive_definite(self):
          return _is_positive_definite(self)
  
      @property
      def is_positive_semidefinite(self):
          return _is_positive_semidefinite(self)
  
      @property
      def is_negative_definite(self):
          return _is_negative_definite(self)
  
      @property
      def is_negative_semidefinite(self):
          return _is_negative_semidefinite(self)
  
      @property
      def is_indefinite(self):
          return _is_indefinite(self)
  
      def jordan_form(self, calc_transform=True, **kwargs):
          return _jordan_form(self, calc_transform=calc_transform, **kwargs)
  
      def left_eigenvects(self, **flags):
          return _left_eigenvects(self, **flags)
  
      def singular_values(self):
          return _singular_values(self)
  
      eigenvals.__doc__                  = _eigenvals.__doc__
      eigenvects.__doc__                 = _eigenvects.__doc__
      is_diagonalizable.__doc__          = _is_diagonalizable.__doc__
      diagonalize.__doc__                = _diagonalize.__doc__
      is_positive_definite.__doc__       = _is_positive_definite.__doc__
      is_positive_semidefinite.__doc__   = _is_positive_semidefinite.__doc__
      is_negative_definite.__doc__       = _is_negative_definite.__doc__
      is_negative_semidefinite.__doc__   = _is_negative_semidefinite.__doc__
      is_indefinite.__doc__              = _is_indefinite.__doc__
      jordan_form.__doc__                = _jordan_form.__doc__
      left_eigenvects.__doc__            = _left_eigenvects.__doc__
      singular_values.__doc__            = _singular_values.__doc__
      bidiagonalize.__doc__              = _bidiagonalize.__doc__
      bidiagonal_decomposition.__doc__   = _bidiagonal_decomposition.__doc__
  
  
  class MatrixCalculus(MatrixCommon):
      """Provides calculus-related matrix operations."""
  
      def diff(self, *args, **kwargs):
          """Calculate the derivative of each element in the matrix.
          ``args`` will be passed to the ``integrate`` function.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> from sympy.abc import x, y
          >>> M = Matrix([[x, y], [1, 0]])
          >>> M.diff(x)
          Matrix([
          [1, 0],
          [0, 0]])
  
          See Also
          ========
  
          integrate
          limit
          """
          # XXX this should be handled here rather than in Derivative
          from sympy.tensor.array.array_derivatives import ArrayDerivative
          kwargs.setdefault('evaluate', True)
          deriv = ArrayDerivative(self, *args, evaluate=True)
          if not isinstance(self, Basic):
              return deriv.as_mutable()
          else:
              return deriv
  
      def _eval_derivative(self, arg):
          return self.applyfunc(lambda x: x.diff(arg))
  
      def integrate(self, *args, **kwargs):
          """Integrate each element of the matrix.  ``args`` will
          be passed to the ``integrate`` function.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> from sympy.abc import x, y
          >>> M = Matrix([[x, y], [1, 0]])
          >>> M.integrate((x, ))
          Matrix([
          [x**2/2, x*y],
          [     x,   0]])
          >>> M.integrate((x, 0, 2))
          Matrix([
          [2, 2*y],
          [2,   0]])
  
          See Also
          ========
  
          limit
          diff
          """
          return self.applyfunc(lambda x: x.integrate(*args, **kwargs))
  
      def jacobian(self, X):
          """Calculates the Jacobian matrix (derivative of a vector-valued function).
  
          Parameters
          ==========
  
          ``self`` : vector of expressions representing functions f_i(x_1, ..., x_n).
          X : set of x_i's in order, it can be a list or a Matrix
  
          Both ``self`` and X can be a row or a column matrix in any order
          (i.e., jacobian() should always work).
  
          Examples
          ========
  
          >>> from sympy import sin, cos, Matrix
          >>> from sympy.abc import rho, phi
          >>> X = Matrix([rho*cos(phi), rho*sin(phi), rho**2])
          >>> Y = Matrix([rho, phi])
          >>> X.jacobian(Y)
          Matrix([
          [cos(phi), -rho*sin(phi)],
          [sin(phi),  rho*cos(phi)],
          [   2*rho,             0]])
          >>> X = Matrix([rho*cos(phi), rho*sin(phi)])
          >>> X.jacobian(Y)
          Matrix([
          [cos(phi), -rho*sin(phi)],
          [sin(phi),  rho*cos(phi)]])
  
          See Also
          ========
  
          hessian
          wronskian
          """
          if not isinstance(X, MatrixBase):
              X = self._new(X)
          # Both X and ``self`` can be a row or a column matrix, so we need to make
          # sure all valid combinations work, but everything else fails:
          if self.shape[0] == 1:
              m = self.shape[1]
          elif self.shape[1] == 1:
              m = self.shape[0]
          else:
              raise TypeError("``self`` must be a row or a column matrix")
          if X.shape[0] == 1:
              n = X.shape[1]
          elif X.shape[1] == 1:
              n = X.shape[0]
          else:
              raise TypeError("X must be a row or a column matrix")
  
          # m is the number of functions and n is the number of variables
          # computing the Jacobian is now easy:
          return self._new(m, n, lambda j, i: self[j].diff(X[i]))
  
      def limit(self, *args):
          """Calculate the limit of each element in the matrix.
          ``args`` will be passed to the ``limit`` function.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> from sympy.abc import x, y
          >>> M = Matrix([[x, y], [1, 0]])
          >>> M.limit(x, 2)
          Matrix([
          [2, y],
          [1, 0]])
  
          See Also
          ========
  
          integrate
          diff
          """
          return self.applyfunc(lambda x: x.limit(*args))
  
  
  # https://github.com/sympy/sympy/pull/12854
  class MatrixDeprecated(MatrixCommon):
      """A class to house deprecated matrix methods."""
      def _legacy_array_dot(self, b):
          """Compatibility function for deprecated behavior of ``matrix.dot(vector)``
          """
          from .dense import Matrix
  
          if not isinstance(b, MatrixBase):
              if is_sequence(b):
                  if len(b) != self.cols and len(b) != self.rows:
                      raise ShapeError(
                          "Dimensions incorrect for dot product: %s, %s" % (
                              self.shape, len(b)))
                  return self.dot(Matrix(b))
              else:
                  raise TypeError(
                      "`b` must be an ordered iterable or Matrix, not %s." %
                      type(b))
  
          mat = self
          if mat.cols == b.rows:
              if b.cols != 1:
                  mat = mat.T
                  b = b.T
              prod = flatten((mat * b).tolist())
              return prod
          if mat.cols == b.cols:
              return mat.dot(b.T)
          elif mat.rows == b.rows:
              return mat.T.dot(b)
          else:
              raise ShapeError("Dimensions incorrect for dot product: %s, %s" % (
                  self.shape, b.shape))
  
  
      def berkowitz_charpoly(self, x=Dummy('lambda'), simplify=_simplify):
          return self.charpoly(x=x)
  
      def berkowitz_det(self):
          """Computes determinant using Berkowitz method.
  
          See Also
          ========
  
          det
          berkowitz
          """
          return self.det(method='berkowitz')
  
      def berkowitz_eigenvals(self, **flags):
          """Computes eigenvalues of a Matrix using Berkowitz method.
  
          See Also
          ========
  
          berkowitz
          """
          return self.eigenvals(**flags)
  
      def berkowitz_minors(self):
          """Computes principal minors using Berkowitz method.
  
          See Also
          ========
  
          berkowitz
          """
          sign, minors = self.one, []
  
          for poly in self.berkowitz():
              minors.append(sign * poly[-1])
              sign = -sign
  
          return tuple(minors)
  
      def berkowitz(self):
          from sympy.matrices import zeros
          berk = ((1,),)
          if not self:
              return berk
  
          if not self.is_square:
              raise NonSquareMatrixError()
  
          A, N = self, self.rows
          transforms = [0] * (N - 1)
  
          for n in range(N, 1, -1):
              T, k = zeros(n + 1, n), n - 1
  
              R, C = -A[k, :k], A[:k, k]
              A, a = A[:k, :k], -A[k, k]
  
              items = [C]
  
              for i in range(0, n - 2):
                  items.append(A * items[i])
  
              for i, B in enumerate(items):
                  items[i] = (R * B)[0, 0]
  
              items = [self.one, a] + items
  
              for i in range(n):
                  T[i:, i] = items[:n - i + 1]
  
              transforms[k - 1] = T
  
          polys = [self._new([self.one, -A[0, 0]])]
  
          for i, T in enumerate(transforms):
              polys.append(T * polys[i])
  
          return berk + tuple(map(tuple, polys))
  
      def cofactorMatrix(self, method="berkowitz"):
          return self.cofactor_matrix(method=method)
  
      def det_bareis(self):
          return _det_bareiss(self)
  
      def det_LU_decomposition(self):
          """Compute matrix determinant using LU decomposition
  
  
          Note that this method fails if the LU decomposition itself
          fails. In particular, if the matrix has no inverse this method
          will fail.
  
          TODO: Implement algorithm for sparse matrices (SFF),
          http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps.
  
          See Also
          ========
  
  
          det
          det_bareiss
          berkowitz_det
          """
          return self.det(method='lu')
  
      def jordan_cell(self, eigenval, n):
          return self.jordan_block(size=n, eigenvalue=eigenval)
  
      def jordan_cells(self, calc_transformation=True):
          P, J = self.jordan_form()
          return P, J.get_diag_blocks()
  
      def minorEntry(self, i, j, method="berkowitz"):
          return self.minor(i, j, method=method)
  
      def minorMatrix(self, i, j):
          return self.minor_submatrix(i, j)
  
      def permuteBkwd(self, perm):
          """Permute the rows of the matrix with the given permutation in reverse."""
          return self.permute_rows(perm, direction='backward')
  
      def permuteFwd(self, perm):
          """Permute the rows of the matrix with the given permutation."""
          return self.permute_rows(perm, direction='forward')
  
  
  @Mul._kind_dispatcher.register(_NumberKind, MatrixKind)
  def num_mat_mul(k1, k2):
      """
      Return MatrixKind. The element kind is selected by recursive dispatching.
      Do not need to dispatch in reversed order because KindDispatcher
      searches for this automatically.
      """
      # Deal with Mul._kind_dispatcher's commutativity
      # XXX: this function is called with either k1 or k2 as MatrixKind because
      # the Mul kind dispatcher is commutative. Maybe it shouldn't be. Need to
      # swap the args here because NumberKind does not have an element_kind
      # attribute.
      if not isinstance(k2, MatrixKind):
          k1, k2 = k2, k1
      elemk = Mul._kind_dispatcher(k1, k2.element_kind)
      return MatrixKind(elemk)
  
  
  @Mul._kind_dispatcher.register(MatrixKind, MatrixKind)
  def mat_mat_mul(k1, k2):
      """
      Return MatrixKind. The element kind is selected by recursive dispatching.
      """
      elemk = Mul._kind_dispatcher(k1.element_kind, k2.element_kind)
      return MatrixKind(elemk)
  
  
  class MatrixBase(MatrixDeprecated,
                   MatrixCalculus,
                   MatrixEigen,
                   MatrixCommon,
                   Printable):
      """Base class for matrix objects."""
      # Added just for numpy compatibility
      __array_priority__ = 11
  
      is_Matrix = True
      _class_priority = 3
      _sympify = staticmethod(sympify)
      zero = S.Zero
      one = S.One
  
      @property
      def kind(self) -> MatrixKind:
          elem_kinds = set(e.kind for e in self.flat())
          if len(elem_kinds) == 1:
              elemkind, = elem_kinds
          else:
              elemkind = UndefinedKind
          return MatrixKind(elemkind)
  
      def flat(self):
          return [self[i, j] for i in range(self.rows) for j in range(self.cols)]
  
      def __array__(self, dtype=object):
          from .dense import matrix2numpy
          return matrix2numpy(self, dtype=dtype)
  
      def __len__(self):
          """Return the number of elements of ``self``.
  
          Implemented mainly so bool(Matrix()) == False.
          """
          return self.rows * self.cols
  
      def _matrix_pow_by_jordan_blocks(self, num):
          from sympy.matrices import diag, MutableMatrix
  
          def jordan_cell_power(jc, n):
              N = jc.shape[0]
              l = jc[0,0]
              if l.is_zero:
                  if N == 1 and n.is_nonnegative:
                      jc[0,0] = l**n
                  elif not (n.is_integer and n.is_nonnegative):
                      raise NonInvertibleMatrixError("Non-invertible matrix can only be raised to a nonnegative integer")
                  else:
                      for i in range(N):
                          jc[0,i] = KroneckerDelta(i, n)
              else:
                  for i in range(N):
                      bn = binomial(n, i)
                      if isinstance(bn, binomial):
                          bn = bn._eval_expand_func()
                      jc[0,i] = l**(n-i)*bn
              for i in range(N):
                  for j in range(1, N-i):
                      jc[j,i+j] = jc [j-1,i+j-1]
  
          P, J = self.jordan_form()
          jordan_cells = J.get_diag_blocks()
          # Make sure jordan_cells matrices are mutable:
          jordan_cells = [MutableMatrix(j) for j in jordan_cells]
          for j in jordan_cells:
              jordan_cell_power(j, num)
          return self._new(P.multiply(diag(*jordan_cells))
                  .multiply(P.inv()))
  
      def __str__(self):
          if S.Zero in self.shape:
              return 'Matrix(%s, %s, [])' % (self.rows, self.cols)
          return "Matrix(%s)" % str(self.tolist())
  
      def _format_str(self, printer=None):
          if not printer:
              printer = StrPrinter()
          # Handle zero dimensions:
          if S.Zero in self.shape:
              return 'Matrix(%s, %s, [])' % (self.rows, self.cols)
          if self.rows == 1:
              return "Matrix([%s])" % self.table(printer, rowsep=',\n')
          return "Matrix([\n%s])" % self.table(printer, rowsep=',\n')
  
      @classmethod
      def irregular(cls, ntop, *matrices, **kwargs):
        """Return a matrix filled by the given matrices which
        are listed in order of appearance from left to right, top to
        bottom as they first appear in the matrix. They must fill the
        matrix completely.
  
        Examples
        ========
  
        >>> from sympy import ones, Matrix
        >>> Matrix.irregular(3, ones(2,1), ones(3,3)*2, ones(2,2)*3,
        ...   ones(1,1)*4, ones(2,2)*5, ones(1,2)*6, ones(1,2)*7)
        Matrix([
          [1, 2, 2, 2, 3, 3],
          [1, 2, 2, 2, 3, 3],
          [4, 2, 2, 2, 5, 5],
          [6, 6, 7, 7, 5, 5]])
        """
        ntop = as_int(ntop)
        # make sure we are working with explicit matrices
        b = [i.as_explicit() if hasattr(i, 'as_explicit') else i
            for i in matrices]
        q = list(range(len(b)))
        dat = [i.rows for i in b]
        active = [q.pop(0) for _ in range(ntop)]
        cols = sum([b[i].cols for i in active])
        rows = []
        while any(dat):
            r = []
            for a, j in enumerate(active):
                r.extend(b[j][-dat[j], :])
                dat[j] -= 1
                if dat[j] == 0 and q:
                    active[a] = q.pop(0)
            if len(r) != cols:
              raise ValueError(filldedent('''
                  Matrices provided do not appear to fill
                  the space completely.'''))
            rows.append(r)
        return cls._new(rows)
  
      @classmethod
      def _handle_ndarray(cls, arg):
          # NumPy array or matrix or some other object that implements
          # __array__. So let's first use this method to get a
          # numpy.array() and then make a Python list out of it.
          arr = arg.__array__()
          if len(arr.shape) == 2:
              rows, cols = arr.shape[0], arr.shape[1]
              flat_list = [cls._sympify(i) for i in arr.ravel()]
              return rows, cols, flat_list
          elif len(arr.shape) == 1:
              flat_list = [cls._sympify(i) for i in arr]
              return arr.shape[0], 1, flat_list
          else:
              raise NotImplementedError(
                  "SymPy supports just 1D and 2D matrices")
  
      @classmethod
      def _handle_creation_inputs(cls, *args, **kwargs):
          """Return the number of rows, cols and flat matrix elements.
  
          Examples
          ========
  
          >>> from sympy import Matrix, I
  
          Matrix can be constructed as follows:
  
          * from a nested list of iterables
  
          >>> Matrix( ((1, 2+I), (3, 4)) )
          Matrix([
          [1, 2 + I],
          [3,     4]])
  
          * from un-nested iterable (interpreted as a column)
  
          >>> Matrix( [1, 2] )
          Matrix([
          [1],
          [2]])
  
          * from un-nested iterable with dimensions
  
          >>> Matrix(1, 2, [1, 2] )
          Matrix([[1, 2]])
  
          * from no arguments (a 0 x 0 matrix)
  
          >>> Matrix()
          Matrix(0, 0, [])
  
          * from a rule
  
          >>> Matrix(2, 2, lambda i, j: i/(j + 1) )
          Matrix([
          [0,   0],
          [1, 1/2]])
  
          See Also
          ========
          irregular - filling a matrix with irregular blocks
          """
          from sympy.matrices import SparseMatrix
          from sympy.matrices.expressions.matexpr import MatrixSymbol
          from sympy.matrices.expressions.blockmatrix import BlockMatrix
  
          flat_list = None
  
          if len(args) == 1:
              # Matrix(SparseMatrix(...))
              if isinstance(args[0], SparseMatrix):
                  return args[0].rows, args[0].cols, flatten(args[0].tolist())
  
              # Matrix(Matrix(...))
              elif isinstance(args[0], MatrixBase):
                  return args[0].rows, args[0].cols, args[0].flat()
  
              # Matrix(MatrixSymbol('X', 2, 2))
              elif isinstance(args[0], Basic) and args[0].is_Matrix:
                  return args[0].rows, args[0].cols, args[0].as_explicit().flat()
  
              elif isinstance(args[0], mp.matrix):
                  M = args[0]
                  flat_list = [cls._sympify(x) for x in M]
                  return M.rows, M.cols, flat_list
  
              # Matrix(numpy.ones((2, 2)))
              elif hasattr(args[0], "__array__"):
                  return cls._handle_ndarray(args[0])
  
              # Matrix([1, 2, 3]) or Matrix([[1, 2], [3, 4]])
              elif is_sequence(args[0]) \
                      and not isinstance(args[0], DeferredVector):
                  dat = list(args[0])
                  ismat = lambda i: isinstance(i, MatrixBase) and (
                      evaluate or
                      isinstance(i, BlockMatrix) or
                      isinstance(i, MatrixSymbol))
                  raw = lambda i: is_sequence(i) and not ismat(i)
                  evaluate = kwargs.get('evaluate', True)
  
  
                  if evaluate:
  
                      def make_explicit(x):
                          """make Block and Symbol explicit"""
                          if isinstance(x, BlockMatrix):
                              return x.as_explicit()
                          elif isinstance(x, MatrixSymbol) and all(_.is_Integer for _ in x.shape):
                              return x.as_explicit()
                          else:
                              return x
  
                      def make_explicit_row(row):
                          # Could be list or could be list of lists
                          if isinstance(row, (list, tuple)):
                              return [make_explicit(x) for x in row]
                          else:
                              return make_explicit(row)
  
                      if isinstance(dat, (list, tuple)):
                          dat = [make_explicit_row(row) for row in dat]
  
                  if dat in ([], [[]]):
                      rows = cols = 0
                      flat_list = []
                  elif not any(raw(i) or ismat(i) for i in dat):
                      # a column as a list of values
                      flat_list = [cls._sympify(i) for i in dat]
                      rows = len(flat_list)
                      cols = 1 if rows else 0
                  elif evaluate and all(ismat(i) for i in dat):
                      # a column as a list of matrices
                      ncol = {i.cols for i in dat if any(i.shape)}
                      if ncol:
                          if len(ncol) != 1:
                              raise ValueError('mismatched dimensions')
                          flat_list = [_ for i in dat for r in i.tolist() for _ in r]
                          cols = ncol.pop()
                          rows = len(flat_list)//cols
                      else:
                          rows = cols = 0
                          flat_list = []
                  elif evaluate and any(ismat(i) for i in dat):
                      ncol = set()
                      flat_list = []
                      for i in dat:
                          if ismat(i):
                              flat_list.extend(
                                  [k for j in i.tolist() for k in j])
                              if any(i.shape):
                                  ncol.add(i.cols)
                          elif raw(i):
                              if i:
                                  ncol.add(len(i))
                                  flat_list.extend([cls._sympify(ij) for ij in i])
                          else:
                              ncol.add(1)
                              flat_list.append(i)
                          if len(ncol) > 1:
                              raise ValueError('mismatched dimensions')
                      cols = ncol.pop()
                      rows = len(flat_list)//cols
                  else:
                      # list of lists; each sublist is a logical row
                      # which might consist of many rows if the values in
                      # the row are matrices
                      flat_list = []
                      ncol = set()
                      rows = cols = 0
                      for row in dat:
                          if not is_sequence(row) and \
                                  not getattr(row, 'is_Matrix', False):
                              raise ValueError('expecting list of lists')
  
                          if hasattr(row, '__array__'):
                              if 0 in row.shape:
                                  continue
                          elif not row:
                              continue
  
                          if evaluate and all(ismat(i) for i in row):
                              r, c, flatT = cls._handle_creation_inputs(
                                  [i.T for i in row])
                              T = reshape(flatT, [c])
                              flat = \
                                  [T[i][j] for j in range(c) for i in range(r)]
                              r, c = c, r
                          else:
                              r = 1
                              if getattr(row, 'is_Matrix', False):
                                  c = 1
                                  flat = [row]
                              else:
                                  c = len(row)
                                  flat = [cls._sympify(i) for i in row]
                          ncol.add(c)
                          if len(ncol) > 1:
                              raise ValueError('mismatched dimensions')
                          flat_list.extend(flat)
                          rows += r
                      cols = ncol.pop() if ncol else 0
  
          elif len(args) == 3:
              rows = as_int(args[0])
              cols = as_int(args[1])
  
              if rows < 0 or cols < 0:
                  raise ValueError("Cannot create a {} x {} matrix. "
                                   "Both dimensions must be positive".format(rows, cols))
  
              # Matrix(2, 2, lambda i, j: i+j)
              if len(args) == 3 and isinstance(args[2], Callable):
                  op = args[2]
                  flat_list = []
                  for i in range(rows):
                      flat_list.extend(
                          [cls._sympify(op(cls._sympify(i), cls._sympify(j)))
                           for j in range(cols)])
  
              # Matrix(2, 2, [1, 2, 3, 4])
              elif len(args) == 3 and is_sequence(args[2]):
                  flat_list = args[2]
                  if len(flat_list) != rows * cols:
                      raise ValueError(
                          'List length should be equal to rows*columns')
                  flat_list = [cls._sympify(i) for i in flat_list]
  
  
          # Matrix()
          elif len(args) == 0:
              # Empty Matrix
              rows = cols = 0
              flat_list = []
  
          if flat_list is None:
              raise TypeError(filldedent('''
                  Data type not understood; expecting list of lists
                  or lists of values.'''))
  
          return rows, cols, flat_list
  
      def _setitem(self, key, value):
          """Helper to set value at location given by key.
  
          Examples
          ========
  
          >>> from sympy import Matrix, I, zeros, ones
          >>> m = Matrix(((1, 2+I), (3, 4)))
          >>> m
          Matrix([
          [1, 2 + I],
          [3,     4]])
          >>> m[1, 0] = 9
          >>> m
          Matrix([
          [1, 2 + I],
          [9,     4]])
          >>> m[1, 0] = [[0, 1]]
  
          To replace row r you assign to position r*m where m
          is the number of columns:
  
          >>> M = zeros(4)
          >>> m = M.cols
          >>> M[3*m] = ones(1, m)*2; M
          Matrix([
          [0, 0, 0, 0],
          [0, 0, 0, 0],
          [0, 0, 0, 0],
          [2, 2, 2, 2]])
  
          And to replace column c you can assign to position c:
  
          >>> M[2] = ones(m, 1)*4; M
          Matrix([
          [0, 0, 4, 0],
          [0, 0, 4, 0],
          [0, 0, 4, 0],
          [2, 2, 4, 2]])
          """
          from .dense import Matrix
  
          is_slice = isinstance(key, slice)
          i, j = key = self.key2ij(key)
          is_mat = isinstance(value, MatrixBase)
          if isinstance(i, slice) or isinstance(j, slice):
              if is_mat:
                  self.copyin_matrix(key, value)
                  return
              if not isinstance(value, Expr) and is_sequence(value):
                  self.copyin_list(key, value)
                  return
              raise ValueError('unexpected value: %s' % value)
          else:
              if (not is_mat and
                      not isinstance(value, Basic) and is_sequence(value)):
                  value = Matrix(value)
                  is_mat = True
              if is_mat:
                  if is_slice:
                      key = (slice(*divmod(i, self.cols)),
                             slice(*divmod(j, self.cols)))
                  else:
                      key = (slice(i, i + value.rows),
                             slice(j, j + value.cols))
                  self.copyin_matrix(key, value)
              else:
                  return i, j, self._sympify(value)
              return
  
      def add(self, b):
          """Return self + b """
          return self + b
  
      def condition_number(self):
          """Returns the condition number of a matrix.
  
          This is the maximum singular value divided by the minimum singular value
  
          Examples
          ========
  
          >>> from sympy import Matrix, S
          >>> A = Matrix([[1, 0, 0], [0, 10, 0], [0, 0, S.One/10]])
          >>> A.condition_number()
          100
  
          See Also
          ========
  
          singular_values
          """
  
          if not self:
              return self.zero
          singularvalues = self.singular_values()
          return Max(*singularvalues) / Min(*singularvalues)
  
      def copy(self):
          """
          Returns the copy of a matrix.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> A = Matrix(2, 2, [1, 2, 3, 4])
          >>> A.copy()
          Matrix([
          [1, 2],
          [3, 4]])
  
          """
          return self._new(self.rows, self.cols, self.flat())
  
      def cross(self, b):
          r"""
          Return the cross product of ``self`` and ``b`` relaxing the condition
          of compatible dimensions: if each has 3 elements, a matrix of the
          same type and shape as ``self`` will be returned. If ``b`` has the same
          shape as ``self`` then common identities for the cross product (like
          `a \times b = - b \times a`) will hold.
  
          Parameters
          ==========
              b : 3x1 or 1x3 Matrix
  
          See Also
          ========
  
          dot
          multiply
          multiply_elementwise
          """
          from sympy.matrices.expressions.matexpr import MatrixExpr
  
          if not isinstance(b, (MatrixBase, MatrixExpr)):
              raise TypeError(
                  "{} must be a Matrix, not {}.".format(b, type(b)))
  
          if not (self.rows * self.cols == b.rows * b.cols == 3):
              raise ShapeError("Dimensions incorrect for cross product: %s x %s" %
                               ((self.rows, self.cols), (b.rows, b.cols)))
          else:
              return self._new(self.rows, self.cols, (
                  (self[1] * b[2] - self[2] * b[1]),
                  (self[2] * b[0] - self[0] * b[2]),
                  (self[0] * b[1] - self[1] * b[0])))
  
      @property
      def D(self):
          """Return Dirac conjugate (if ``self.rows == 4``).
  
          Examples
          ========
  
          >>> from sympy import Matrix, I, eye
          >>> m = Matrix((0, 1 + I, 2, 3))
          >>> m.D
          Matrix([[0, 1 - I, -2, -3]])
          >>> m = (eye(4) + I*eye(4))
          >>> m[0, 3] = 2
          >>> m.D
          Matrix([
          [1 - I,     0,      0,      0],
          [    0, 1 - I,      0,      0],
          [    0,     0, -1 + I,      0],
          [    2,     0,      0, -1 + I]])
  
          If the matrix does not have 4 rows an AttributeError will be raised
          because this property is only defined for matrices with 4 rows.
  
          >>> Matrix(eye(2)).D
          Traceback (most recent call last):
          ...
          AttributeError: Matrix has no attribute D.
  
          See Also
          ========
  
          sympy.matrices.common.MatrixCommon.conjugate: By-element conjugation
          sympy.matrices.common.MatrixCommon.H: Hermite conjugation
          """
          from sympy.physics.matrices import mgamma
          if self.rows != 4:
              # In Python 3.2, properties can only return an AttributeError
              # so we can't raise a ShapeError -- see commit which added the
              # first line of this inline comment. Also, there is no need
              # for a message since MatrixBase will raise the AttributeError
              raise AttributeError
          return self.H * mgamma(0)
  
      def dot(self, b, hermitian=None, conjugate_convention=None):
          """Return the dot or inner product of two vectors of equal length.
          Here ``self`` must be a ``Matrix`` of size 1 x n or n x 1, and ``b``
          must be either a matrix of size 1 x n, n x 1, or a list/tuple of length n.
          A scalar is returned.
  
          By default, ``dot`` does not conjugate ``self`` or ``b``, even if there are
          complex entries. Set ``hermitian=True`` (and optionally a ``conjugate_convention``)
          to compute the hermitian inner product.
  
          Possible kwargs are ``hermitian`` and ``conjugate_convention``.
  
          If ``conjugate_convention`` is ``"left"``, ``"math"`` or ``"maths"``,
          the conjugate of the first vector (``self``) is used.  If ``"right"``
          or ``"physics"`` is specified, the conjugate of the second vector ``b`` is used.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
          >>> v = Matrix([1, 1, 1])
          >>> M.row(0).dot(v)
          6
          >>> M.col(0).dot(v)
          12
          >>> v = [3, 2, 1]
          >>> M.row(0).dot(v)
          10
  
          >>> from sympy import I
          >>> q = Matrix([1*I, 1*I, 1*I])
          >>> q.dot(q, hermitian=False)
          -3
  
          >>> q.dot(q, hermitian=True)
          3
  
          >>> q1 = Matrix([1, 1, 1*I])
          >>> q.dot(q1, hermitian=True, conjugate_convention="maths")
          1 - 2*I
          >>> q.dot(q1, hermitian=True, conjugate_convention="physics")
          1 + 2*I
  
  
          See Also
          ========
  
          cross
          multiply
          multiply_elementwise
          """
          from .dense import Matrix
  
          if not isinstance(b, MatrixBase):
              if is_sequence(b):
                  if len(b) != self.cols and len(b) != self.rows:
                      raise ShapeError(
                          "Dimensions incorrect for dot product: %s, %s" % (
                              self.shape, len(b)))
                  return self.dot(Matrix(b))
              else:
                  raise TypeError(
                      "`b` must be an ordered iterable or Matrix, not %s." %
                      type(b))
  
          mat = self
          if (1 not in mat.shape) or (1 not in b.shape) :
              sympy_deprecation_warning(
                  """
                  Using the dot method to multiply non-row/column vectors is
                  deprecated. Use * or @ to perform matrix multiplication.
                  """,
                  deprecated_since_version="1.2",
                  active_deprecations_target="deprecated-matrix-dot-non-vector")
              return mat._legacy_array_dot(b)
          if len(mat) != len(b):
              raise ShapeError("Dimensions incorrect for dot product: %s, %s" % (self.shape, b.shape))
          n = len(mat)
          if mat.shape != (1, n):
              mat = mat.reshape(1, n)
          if b.shape != (n, 1):
              b = b.reshape(n, 1)
  
          # Now ``mat`` is a row vector and ``b`` is a column vector.
  
          # If it so happens that only conjugate_convention is passed
          # then automatically set hermitian to True. If only hermitian
          # is true but no conjugate_convention is not passed then
          # automatically set it to ``"maths"``
  
          if conjugate_convention is not None and hermitian is None:
              hermitian = True
          if hermitian and conjugate_convention is None:
              conjugate_convention = "maths"
  
          if hermitian == True:
              if conjugate_convention in ("maths", "left", "math"):
                  mat = mat.conjugate()
              elif conjugate_convention in ("physics", "right"):
                  b = b.conjugate()
              else:
                  raise ValueError("Unknown conjugate_convention was entered."
                                   " conjugate_convention must be one of the"
                                   " following: math, maths, left, physics or right.")
          return (mat * b)[0]
  
      def dual(self):
          """Returns the dual of a matrix, which is:
  
          ``(1/2)*levicivita(i, j, k, l)*M(k, l)`` summed over indices `k` and `l`
  
          Since the levicivita method is anti_symmetric for any pairwise
          exchange of indices, the dual of a symmetric matrix is the zero
          matrix. Strictly speaking the dual defined here assumes that the
          'matrix' `M` is a contravariant anti_symmetric second rank tensor,
          so that the dual is a covariant second rank tensor.
  
          """
          from sympy.matrices import zeros
  
          M, n = self[:, :], self.rows
          work = zeros(n)
          if self.is_symmetric():
              return work
  
          for i in range(1, n):
              for j in range(1, n):
                  acum = 0
                  for k in range(1, n):
                      acum += LeviCivita(i, j, 0, k) * M[0, k]
                  work[i, j] = acum
                  work[j, i] = -acum
  
          for l in range(1, n):
              acum = 0
              for a in range(1, n):
                  for b in range(1, n):
                      acum += LeviCivita(0, l, a, b) * M[a, b]
              acum /= 2
              work[0, l] = -acum
              work[l, 0] = acum
  
          return work
  
      def _eval_matrix_exp_jblock(self):
          """A helper function to compute an exponential of a Jordan block
          matrix
  
          Examples
          ========
  
          >>> from sympy import Symbol, Matrix
          >>> l = Symbol('lamda')
  
          A trivial example of 1*1 Jordan block:
  
          >>> m = Matrix.jordan_block(1, l)
          >>> m._eval_matrix_exp_jblock()
          Matrix([[exp(lamda)]])
  
          An example of 3*3 Jordan block:
  
          >>> m = Matrix.jordan_block(3, l)
          >>> m._eval_matrix_exp_jblock()
          Matrix([
          [exp(lamda), exp(lamda), exp(lamda)/2],
          [         0, exp(lamda),   exp(lamda)],
          [         0,          0,   exp(lamda)]])
  
          References
          ==========
  
          .. [1] https://en.wikipedia.org/wiki/Matrix_function#Jordan_decomposition
          """
          size = self.rows
          l = self[0, 0]
          exp_l = exp(l)
  
          bands = {i: exp_l / factorial(i) for i in range(size)}
  
          from .sparsetools import banded
          return self.__class__(banded(size, bands))
  
  
      def analytic_func(self, f, x):
          """
          Computes f(A) where A is a Square Matrix
          and f is an analytic function.
  
          Examples
          ========
  
          >>> from sympy import Symbol, Matrix, S, log
  
          >>> x = Symbol('x')
          >>> m = Matrix([[S(5)/4, S(3)/4], [S(3)/4, S(5)/4]])
          >>> f = log(x)
          >>> m.analytic_func(f, x)
          Matrix([
          [     0, log(2)],
          [log(2),      0]])
  
          Parameters
          ==========
  
          f : Expr
              Analytic Function
          x : Symbol
              parameter of f
  
          """
  
          f, x = _sympify(f), _sympify(x)
          if not self.is_square:
              raise NonSquareMatrixError
          if not x.is_symbol:
              raise ValueError("{} must be a symbol.".format(x))
          if x not in f.free_symbols:
              raise ValueError(
                  "{} must be a parameter of {}.".format(x, f))
          if x in self.free_symbols:
              raise ValueError(
                  "{} must not be a parameter of {}.".format(x, self))
  
          eigen = self.eigenvals()
          max_mul = max(eigen.values())
          derivative = {}
          dd = f
          for i in range(max_mul - 1):
              dd = diff(dd, x)
              derivative[i + 1] = dd
          n = self.shape[0]
          r = self.zeros(n)
          f_val = self.zeros(n, 1)
          row = 0
  
          for i in eigen:
              mul = eigen[i]
              f_val[row] = f.subs(x, i)
              if f_val[row].is_number and not f_val[row].is_complex:
                  raise ValueError(
                      "Cannot evaluate the function because the "
                      "function {} is not analytic at the given "
                      "eigenvalue {}".format(f, f_val[row]))
              val = 1
              for a in range(n):
                  r[row, a] = val
                  val *= i
              if mul > 1:
                  coe = [1 for ii in range(n)]
                  deri = 1
                  while mul > 1:
                      row = row + 1
                      mul -= 1
                      d_i = derivative[deri].subs(x, i)
                      if d_i.is_number and not d_i.is_complex:
                          raise ValueError(
                              "Cannot evaluate the function because the "
                              "derivative {} is not analytic at the given "
                              "eigenvalue {}".format(derivative[deri], d_i))
                      f_val[row] = d_i
                      for a in range(n):
                          if a - deri + 1 <= 0:
                              r[row, a] = 0
                              coe[a] = 0
                              continue
                          coe[a] = coe[a]*(a - deri + 1)
                          r[row, a] = coe[a]*pow(i, a - deri)
                      deri += 1
              row += 1
          c = r.solve(f_val)
          ans = self.zeros(n)
          pre = self.eye(n)
          for i in range(n):
              ans = ans + c[i]*pre
              pre *= self
          return ans
  
  
      def exp(self):
  
          """Return the exponential of a square matrix
  
          Examples
          ========
  
          >>> from sympy import Symbol, Matrix
  
          >>> t = Symbol('t')
          >>> m = Matrix([[0, 1], [-1, 0]]) * t
          >>> m.exp()
          Matrix([
          [    exp(I*t)/2 + exp(-I*t)/2, -I*exp(I*t)/2 + I*exp(-I*t)/2],
          [I*exp(I*t)/2 - I*exp(-I*t)/2,      exp(I*t)/2 + exp(-I*t)/2]])
          """
          if not self.is_square:
              raise NonSquareMatrixError(
                  "Exponentiation is valid only for square matrices")
          try:
              P, J = self.jordan_form()
              cells = J.get_diag_blocks()
          except MatrixError:
              raise NotImplementedError(
                  "Exponentiation is implemented only for matrices for which the Jordan normal form can be computed")
  
          blocks = [cell._eval_matrix_exp_jblock() for cell in cells]
          from sympy.matrices import diag
          eJ = diag(*blocks)
          # n = self.rows
          ret = P.multiply(eJ, dotprodsimp=None).multiply(P.inv(), dotprodsimp=None)
          if all(value.is_real for value in self.values()):
              return type(self)(re(ret))
          else:
              return type(self)(ret)
  
      def _eval_matrix_log_jblock(self):
          """Helper function to compute logarithm of a jordan block.
  
          Examples
          ========
  
          >>> from sympy import Symbol, Matrix
          >>> l = Symbol('lamda')
  
          A trivial example of 1*1 Jordan block:
  
          >>> m = Matrix.jordan_block(1, l)
          >>> m._eval_matrix_log_jblock()
          Matrix([[log(lamda)]])
  
          An example of 3*3 Jordan block:
  
          >>> m = Matrix.jordan_block(3, l)
          >>> m._eval_matrix_log_jblock()
          Matrix([
          [log(lamda),    1/lamda, -1/(2*lamda**2)],
          [         0, log(lamda),         1/lamda],
          [         0,          0,      log(lamda)]])
          """
          size = self.rows
          l = self[0, 0]
  
          if l.is_zero:
              raise MatrixError(
                  'Could not take logarithm or reciprocal for the given '
                  'eigenvalue {}'.format(l))
  
          bands = {0: log(l)}
          for i in range(1, size):
              bands[i] = -((-l) ** -i) / i
  
          from .sparsetools import banded
          return self.__class__(banded(size, bands))
  
      def log(self, simplify=cancel):
          """Return the logarithm of a square matrix
  
          Parameters
          ==========
  
          simplify : function, bool
              The function to simplify the result with.
  
              Default is ``cancel``, which is effective to reduce the
              expression growing for taking reciprocals and inverses for
              symbolic matrices.
  
          Examples
          ========
  
          >>> from sympy import S, Matrix
  
          Examples for positive-definite matrices:
  
          >>> m = Matrix([[1, 1], [0, 1]])
          >>> m.log()
          Matrix([
          [0, 1],
          [0, 0]])
  
          >>> m = Matrix([[S(5)/4, S(3)/4], [S(3)/4, S(5)/4]])
          >>> m.log()
          Matrix([
          [     0, log(2)],
          [log(2),      0]])
  
          Examples for non positive-definite matrices:
  
          >>> m = Matrix([[S(3)/4, S(5)/4], [S(5)/4, S(3)/4]])
          >>> m.log()
          Matrix([
          [         I*pi/2, log(2) - I*pi/2],
          [log(2) - I*pi/2,          I*pi/2]])
  
          >>> m = Matrix(
          ...     [[0, 0, 0, 1],
          ...      [0, 0, 1, 0],
          ...      [0, 1, 0, 0],
          ...      [1, 0, 0, 0]])
          >>> m.log()
          Matrix([
          [ I*pi/2,       0,       0, -I*pi/2],
          [      0,  I*pi/2, -I*pi/2,       0],
          [      0, -I*pi/2,  I*pi/2,       0],
          [-I*pi/2,       0,       0,  I*pi/2]])
          """
          if not self.is_square:
              raise NonSquareMatrixError(
                  "Logarithm is valid only for square matrices")
  
          try:
              if simplify:
                  P, J = simplify(self).jordan_form()
              else:
                  P, J = self.jordan_form()
  
              cells = J.get_diag_blocks()
          except MatrixError:
              raise NotImplementedError(
                  "Logarithm is implemented only for matrices for which "
                  "the Jordan normal form can be computed")
  
          blocks = [
              cell._eval_matrix_log_jblock()
              for cell in cells]
          from sympy.matrices import diag
          eJ = diag(*blocks)
  
          if simplify:
              ret = simplify(P * eJ * simplify(P.inv()))
              ret = self.__class__(ret)
          else:
              ret = P * eJ * P.inv()
  
          return ret
  
      def is_nilpotent(self):
          """Checks if a matrix is nilpotent.
  
          A matrix B is nilpotent if for some integer k, B**k is
          a zero matrix.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> a = Matrix([[0, 0, 0], [1, 0, 0], [1, 1, 0]])
          >>> a.is_nilpotent()
          True
  
          >>> a = Matrix([[1, 0, 1], [1, 0, 0], [1, 1, 0]])
          >>> a.is_nilpotent()
          False
          """
          if not self:
              return True
          if not self.is_square:
              raise NonSquareMatrixError(
                  "Nilpotency is valid only for square matrices")
          x = uniquely_named_symbol('x', self, modify=lambda s: '_' + s)
          p = self.charpoly(x)
          if p.args[0] == x ** self.rows:
              return True
          return False
  
      def key2bounds(self, keys):
          """Converts a key with potentially mixed types of keys (integer and slice)
          into a tuple of ranges and raises an error if any index is out of ``self``'s
          range.
  
          See Also
          ========
  
          key2ij
          """
          from sympy.matrices.common import a2idx as a2idx_ # Remove this line after deprecation of a2idx from matrices.py
  
          islice, jslice = [isinstance(k, slice) for k in keys]
          if islice:
              if not self.rows:
                  rlo = rhi = 0
              else:
                  rlo, rhi = keys[0].indices(self.rows)[:2]
          else:
              rlo = a2idx_(keys[0], self.rows)
              rhi = rlo + 1
          if jslice:
              if not self.cols:
                  clo = chi = 0
              else:
                  clo, chi = keys[1].indices(self.cols)[:2]
          else:
              clo = a2idx_(keys[1], self.cols)
              chi = clo + 1
          return rlo, rhi, clo, chi
  
      def key2ij(self, key):
          """Converts key into canonical form, converting integers or indexable
          items into valid integers for ``self``'s range or returning slices
          unchanged.
  
          See Also
          ========
  
          key2bounds
          """
          from sympy.matrices.common import a2idx as a2idx_ # Remove this line after deprecation of a2idx from matrices.py
  
          if is_sequence(key):
              if not len(key) == 2:
                  raise TypeError('key must be a sequence of length 2')
              return [a2idx_(i, n) if not isinstance(i, slice) else i
                      for i, n in zip(key, self.shape)]
          elif isinstance(key, slice):
              return key.indices(len(self))[:2]
          else:
              return divmod(a2idx_(key, len(self)), self.cols)
  
      def normalized(self, iszerofunc=_iszero):
          """Return the normalized version of ``self``.
  
          Parameters
          ==========
  
          iszerofunc : Function, optional
              A function to determine whether ``self`` is a zero vector.
              The default ``_iszero`` tests to see if each element is
              exactly zero.
  
          Returns
          =======
  
          Matrix
              Normalized vector form of ``self``.
              It has the same length as a unit vector. However, a zero vector
              will be returned for a vector with norm 0.
  
          Raises
          ======
  
          ShapeError
              If the matrix is not in a vector form.
  
          See Also
          ========
  
          norm
          """
          if self.rows != 1 and self.cols != 1:
              raise ShapeError("A Matrix must be a vector to normalize.")
          norm = self.norm()
          if iszerofunc(norm):
              out = self.zeros(self.rows, self.cols)
          else:
              out = self.applyfunc(lambda i: i / norm)
          return out
  
      def norm(self, ord=None):
          """Return the Norm of a Matrix or Vector.
          In the simplest case this is the geometric size of the vector
          Other norms can be specified by the ord parameter
  
  
          =====  ============================  ==========================
          ord    norm for matrices             norm for vectors
          =====  ============================  ==========================
          None   Frobenius norm                2-norm
          'fro'  Frobenius norm                - does not exist
          inf    maximum row sum               max(abs(x))
          -inf   --                            min(abs(x))
          1      maximum column sum            as below
          -1     --                            as below
          2      2-norm (largest sing. value)  as below
          -2     smallest singular value       as below
          other  - does not exist              sum(abs(x)**ord)**(1./ord)
          =====  ============================  ==========================
  
          Examples
          ========
  
          >>> from sympy import Matrix, Symbol, trigsimp, cos, sin, oo
          >>> x = Symbol('x', real=True)
          >>> v = Matrix([cos(x), sin(x)])
          >>> trigsimp( v.norm() )
          1
          >>> v.norm(10)
          (sin(x)**10 + cos(x)**10)**(1/10)
          >>> A = Matrix([[1, 1], [1, 1]])
          >>> A.norm(1) # maximum sum of absolute values of A is 2
          2
          >>> A.norm(2) # Spectral norm (max of |Ax|/|x| under 2-vector-norm)
          2
          >>> A.norm(-2) # Inverse spectral norm (smallest singular value)
          0
          >>> A.norm() # Frobenius Norm
          2
          >>> A.norm(oo) # Infinity Norm
          2
          >>> Matrix([1, -2]).norm(oo)
          2
          >>> Matrix([-1, 2]).norm(-oo)
          1
  
          See Also
          ========
  
          normalized
          """
          # Row or Column Vector Norms
          vals = list(self.values()) or [0]
          if S.One in self.shape:
              if ord in (2, None):  # Common case sqrt(<x, x>)
                  return sqrt(Add(*(abs(i) ** 2 for i in vals)))
  
              elif ord == 1:  # sum(abs(x))
                  return Add(*(abs(i) for i in vals))
  
              elif ord is S.Infinity:  # max(abs(x))
                  return Max(*[abs(i) for i in vals])
  
              elif ord is S.NegativeInfinity:  # min(abs(x))
                  return Min(*[abs(i) for i in vals])
  
              # Otherwise generalize the 2-norm, Sum(x_i**ord)**(1/ord)
              # Note that while useful this is not mathematically a norm
              try:
                  return Pow(Add(*(abs(i) ** ord for i in vals)), S.One / ord)
              except (NotImplementedError, TypeError):
                  raise ValueError("Expected order to be Number, Symbol, oo")
  
          # Matrix Norms
          else:
              if ord == 1:  # Maximum column sum
                  m = self.applyfunc(abs)
                  return Max(*[sum(m.col(i)) for i in range(m.cols)])
  
              elif ord == 2:  # Spectral Norm
                  # Maximum singular value
                  return Max(*self.singular_values())
  
              elif ord == -2:
                  # Minimum singular value
                  return Min(*self.singular_values())
  
              elif ord is S.Infinity:   # Infinity Norm - Maximum row sum
                  m = self.applyfunc(abs)
                  return Max(*[sum(m.row(i)) for i in range(m.rows)])
  
              elif (ord is None or isinstance(ord,
                                              str) and ord.lower() in
                  ['f', 'fro', 'frobenius', 'vector']):
                  # Reshape as vector and send back to norm function
                  return self.vec().norm(ord=2)
  
              else:
                  raise NotImplementedError("Matrix Norms under development")
  
      def print_nonzero(self, symb="X"):
          """Shows location of non-zero entries for fast shape lookup.
  
          Examples
          ========
  
          >>> from sympy import Matrix, eye
          >>> m = Matrix(2, 3, lambda i, j: i*3+j)
          >>> m
          Matrix([
          [0, 1, 2],
          [3, 4, 5]])
          >>> m.print_nonzero()
          [ XX]
          [XXX]
          >>> m = eye(4)
          >>> m.print_nonzero("x")
          [x   ]
          [ x  ]
          [  x ]
          [   x]
  
          """
          s = []
          for i in range(self.rows):
              line = []
              for j in range(self.cols):
                  if self[i, j] == 0:
                      line.append(" ")
                  else:
                      line.append(str(symb))
              s.append("[%s]" % ''.join(line))
          print('\n'.join(s))
  
      def project(self, v):
          """Return the projection of ``self`` onto the line containing ``v``.
  
          Examples
          ========
  
          >>> from sympy import Matrix, S, sqrt
          >>> V = Matrix([sqrt(3)/2, S.Half])
          >>> x = Matrix([[1, 0]])
          >>> V.project(x)
          Matrix([[sqrt(3)/2, 0]])
          >>> V.project(-x)
          Matrix([[sqrt(3)/2, 0]])
          """
          return v * (self.dot(v) / v.dot(v))
  
      def table(self, printer, rowstart='[', rowend=']', rowsep='\n',
                colsep=', ', align='right'):
          r"""
          String form of Matrix as a table.
  
          ``printer`` is the printer to use for on the elements (generally
          something like StrPrinter())
  
          ``rowstart`` is the string used to start each row (by default '[').
  
          ``rowend`` is the string used to end each row (by default ']').
  
          ``rowsep`` is the string used to separate rows (by default a newline).
  
          ``colsep`` is the string used to separate columns (by default ', ').
  
          ``align`` defines how the elements are aligned. Must be one of 'left',
          'right', or 'center'.  You can also use '<', '>', and '^' to mean the
          same thing, respectively.
  
          This is used by the string printer for Matrix.
  
          Examples
          ========
  
          >>> from sympy import Matrix, StrPrinter
          >>> M = Matrix([[1, 2], [-33, 4]])
          >>> printer = StrPrinter()
          >>> M.table(printer)
          '[  1, 2]\n[-33, 4]'
          >>> print(M.table(printer))
          [  1, 2]
          [-33, 4]
          >>> print(M.table(printer, rowsep=',\n'))
          [  1, 2],
          [-33, 4]
          >>> print('[%s]' % M.table(printer, rowsep=',\n'))
          [[  1, 2],
          [-33, 4]]
          >>> print(M.table(printer, colsep=' '))
          [  1 2]
          [-33 4]
          >>> print(M.table(printer, align='center'))
          [ 1 , 2]
          [-33, 4]
          >>> print(M.table(printer, rowstart='{', rowend='}'))
          {  1, 2}
          {-33, 4}
          """
          # Handle zero dimensions:
          if S.Zero in self.shape:
              return '[]'
          # Build table of string representations of the elements
          res = []
          # Track per-column max lengths for pretty alignment
          maxlen = [0] * self.cols
          for i in range(self.rows):
              res.append([])
              for j in range(self.cols):
                  s = printer._print(self[i, j])
                  res[-1].append(s)
                  maxlen[j] = max(len(s), maxlen[j])
          # Patch strings together
          align = {
              'left': 'ljust',
              'right': 'rjust',
              'center': 'center',
              '<': 'ljust',
              '>': 'rjust',
              '^': 'center',
          }[align]
          for i, row in enumerate(res):
              for j, elem in enumerate(row):
                  row[j] = getattr(elem, align)(maxlen[j])
              res[i] = rowstart + colsep.join(row) + rowend
          return rowsep.join(res)
  
      def rank_decomposition(self, iszerofunc=_iszero, simplify=False):
          return _rank_decomposition(self, iszerofunc=iszerofunc,
                  simplify=simplify)
  
      def cholesky(self, hermitian=True):
          raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix')
  
      def LDLdecomposition(self, hermitian=True):
          raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix')
  
      def LUdecomposition(self, iszerofunc=_iszero, simpfunc=None,
              rankcheck=False):
          return _LUdecomposition(self, iszerofunc=iszerofunc, simpfunc=simpfunc,
                  rankcheck=rankcheck)
  
      def LUdecomposition_Simple(self, iszerofunc=_iszero, simpfunc=None,
              rankcheck=False):
          return _LUdecomposition_Simple(self, iszerofunc=iszerofunc,
                  simpfunc=simpfunc, rankcheck=rankcheck)
  
      def LUdecompositionFF(self):
          return _LUdecompositionFF(self)
  
      def singular_value_decomposition(self):
          return _singular_value_decomposition(self)
  
      def QRdecomposition(self):
          return _QRdecomposition(self)
  
      def upper_hessenberg_decomposition(self):
          return _upper_hessenberg_decomposition(self)
  
      def diagonal_solve(self, rhs):
          return _diagonal_solve(self, rhs)
  
      def lower_triangular_solve(self, rhs):
          raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix')
  
      def upper_triangular_solve(self, rhs):
          raise NotImplementedError('This function is implemented in DenseMatrix or SparseMatrix')
  
      def cholesky_solve(self, rhs):
          return _cholesky_solve(self, rhs)
  
      def LDLsolve(self, rhs):
          return _LDLsolve(self, rhs)
  
      def LUsolve(self, rhs, iszerofunc=_iszero):
          return _LUsolve(self, rhs, iszerofunc=iszerofunc)
  
      def QRsolve(self, b):
          return _QRsolve(self, b)
  
      def gauss_jordan_solve(self, B, freevar=False):
          return _gauss_jordan_solve(self, B, freevar=freevar)
  
      def pinv_solve(self, B, arbitrary_matrix=None):
          return _pinv_solve(self, B, arbitrary_matrix=arbitrary_matrix)
  
      def solve(self, rhs, method='GJ'):
          return _solve(self, rhs, method=method)
  
      def solve_least_squares(self, rhs, method='CH'):
          return _solve_least_squares(self, rhs, method=method)
  
      def pinv(self, method='RD'):
          return _pinv(self, method=method)
  
      def inv_mod(self, m):
          return _inv_mod(self, m)
  
      def inverse_ADJ(self, iszerofunc=_iszero):
          return _inv_ADJ(self, iszerofunc=iszerofunc)
  
      def inverse_BLOCK(self, iszerofunc=_iszero):
          return _inv_block(self, iszerofunc=iszerofunc)
  
      def inverse_GE(self, iszerofunc=_iszero):
          return _inv_GE(self, iszerofunc=iszerofunc)
  
      def inverse_LU(self, iszerofunc=_iszero):
          return _inv_LU(self, iszerofunc=iszerofunc)
  
      def inverse_CH(self, iszerofunc=_iszero):
          return _inv_CH(self, iszerofunc=iszerofunc)
  
      def inverse_LDL(self, iszerofunc=_iszero):
          return _inv_LDL(self, iszerofunc=iszerofunc)
  
      def inverse_QR(self, iszerofunc=_iszero):
          return _inv_QR(self, iszerofunc=iszerofunc)
  
      def inv(self, method=None, iszerofunc=_iszero, try_block_diag=False):
          return _inv(self, method=method, iszerofunc=iszerofunc,
                  try_block_diag=try_block_diag)
  
      def connected_components(self):
          return _connected_components(self)
  
      def connected_components_decomposition(self):
          return _connected_components_decomposition(self)
  
      def strongly_connected_components(self):
          return _strongly_connected_components(self)
  
      def strongly_connected_components_decomposition(self, lower=True):
          return _strongly_connected_components_decomposition(self, lower=lower)
  
      _sage_ = Basic._sage_
  
      rank_decomposition.__doc__     = _rank_decomposition.__doc__
      cholesky.__doc__               = _cholesky.__doc__
      LDLdecomposition.__doc__       = _LDLdecomposition.__doc__
      LUdecomposition.__doc__        = _LUdecomposition.__doc__
      LUdecomposition_Simple.__doc__ = _LUdecomposition_Simple.__doc__
      LUdecompositionFF.__doc__      = _LUdecompositionFF.__doc__
      singular_value_decomposition.__doc__ = _singular_value_decomposition.__doc__
      QRdecomposition.__doc__        = _QRdecomposition.__doc__
      upper_hessenberg_decomposition.__doc__ = _upper_hessenberg_decomposition.__doc__
  
      diagonal_solve.__doc__         = _diagonal_solve.__doc__
      lower_triangular_solve.__doc__ = _lower_triangular_solve.__doc__
      upper_triangular_solve.__doc__ = _upper_triangular_solve.__doc__
      cholesky_solve.__doc__         = _cholesky_solve.__doc__
      LDLsolve.__doc__               = _LDLsolve.__doc__
      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__
      inverse_BLOCK.__doc__          = _inv_block.__doc__
      inv.__doc__                    = _inv.__doc__
  
      connected_components.__doc__   = _connected_components.__doc__
      connected_components_decomposition.__doc__ = \
          _connected_components_decomposition.__doc__
      strongly_connected_components.__doc__   = \
          _strongly_connected_components.__doc__
      strongly_connected_components_decomposition.__doc__ = \
          _strongly_connected_components_decomposition.__doc__
  
  
  @deprecated(
      """
      sympy.matrices.matrices.classof is deprecated. Use
      sympy.matrices.common.classof instead.
      """,
      deprecated_since_version="1.3",
      active_deprecations_target="deprecated-sympy-matrices-classof-a2idx",
  )
  def classof(A, B):
      from sympy.matrices.common import classof as classof_
      return classof_(A, B)
  
  @deprecated(
      """
      sympy.matrices.matrices.a2idx is deprecated. Use
      sympy.matrices.common.a2idx instead.
      """,
      deprecated_since_version="1.3",
      active_deprecations_target="deprecated-sympy-matrices-classof-a2idx",
  )
  def a2idx(j, n=None):
      from sympy.matrices.common import a2idx as a2idx_
      return a2idx_(j, n)
  