text: """
  Basic methods common to all matrices to be used
  when creating more advanced matrices (e.g., matrices over rings,
  etc.).
  """
  
  from collections import defaultdict
  from collections.abc import Iterable
  from inspect import isfunction
  from functools import reduce
  
  from sympy.core.logic import FuzzyBool
  from sympy.assumptions.refine import refine
  from sympy.core import SympifyError, Add
  from sympy.core.basic import Atom
  from sympy.core.compatibility import as_int, is_sequence
  from sympy.core.decorators import call_highest_priority
  from sympy.core.kind import Kind, NumberKind
  from sympy.core.logic import fuzzy_and
  from sympy.core.singleton import S
  from sympy.core.symbol import Symbol
  from sympy.core.sympify import sympify
  from sympy.functions import Abs
  from sympy.polys.polytools import Poly
  from sympy.simplify import simplify as _simplify
  from sympy.simplify.simplify import dotprodsimp as _dotprodsimp
  from sympy.utilities.exceptions import SymPyDeprecationWarning
  from sympy.utilities.iterables import flatten
  from sympy.utilities.misc import filldedent
  from sympy.tensor.array import NDimArray
  
  from .utilities import _get_intermediate_simp_bool
  
  
  class MatrixError(Exception):
      pass
  
  
  class ShapeError(ValueError, MatrixError):
      """Wrong matrix shape"""
      pass
  
  
  class NonSquareMatrixError(ShapeError):
      pass
  
  
  class NonInvertibleMatrixError(ValueError, MatrixError):
      """The matrix in not invertible (division by multidimensional zero error)."""
      pass
  
  
  class NonPositiveDefiniteMatrixError(ValueError, MatrixError):
      """The matrix is not a positive-definite matrix."""
      pass
  
  
  class MatrixRequired:
      """All subclasses of matrix objects must implement the
      required matrix properties listed here."""
      rows = None  # type: int
      cols = None  # type: int
      _simplify = None
  
      @classmethod
      def _new(cls, *args, **kwargs):
          """`_new` must, at minimum, be callable as
          `_new(rows, cols, mat) where mat is a flat list of the
          elements of the matrix."""
          raise NotImplementedError("Subclasses must implement this.")
  
      def __eq__(self, other):
          raise NotImplementedError("Subclasses must implement this.")
  
      def __getitem__(self, key):
          """Implementations of __getitem__ should accept ints, in which
          case the matrix is indexed as a flat list, tuples (i,j) in which
          case the (i,j) entry is returned, slices, or mixed tuples (a,b)
          where a and b are any combintion of slices and integers."""
          raise NotImplementedError("Subclasses must implement this.")
  
      def __len__(self):
          """The total number of entries in the matrix."""
          raise NotImplementedError("Subclasses must implement this.")
  
      @property
      def shape(self):
          raise NotImplementedError("Subclasses must implement this.")
  
  
  class MatrixShaping(MatrixRequired):
      """Provides basic matrix shaping and extracting of submatrices"""
  
      def _eval_col_del(self, col):
          def entry(i, j):
              return self[i, j] if j < col else self[i, j + 1]
          return self._new(self.rows, self.cols - 1, entry)
  
      def _eval_col_insert(self, pos, other):
  
          def entry(i, j):
              if j < pos:
                  return self[i, j]
              elif pos <= j < pos + other.cols:
                  return other[i, j - pos]
              return self[i, j - other.cols]
  
          return self._new(self.rows, self.cols + other.cols,
                           lambda i, j: entry(i, j))
  
      def _eval_col_join(self, other):
          rows = self.rows
  
          def entry(i, j):
              if i < rows:
                  return self[i, j]
              return other[i - rows, j]
  
          return classof(self, other)._new(self.rows + other.rows, self.cols,
                                           lambda i, j: entry(i, j))
  
      def _eval_extract(self, rowsList, colsList):
          mat = list(self)
          cols = self.cols
          indices = (i * cols + j for i in rowsList for j in colsList)
          return self._new(len(rowsList), len(colsList),
                           list(mat[i] for i in indices))
  
      def _eval_get_diag_blocks(self):
          sub_blocks = []
  
          def recurse_sub_blocks(M):
              i = 1
              while i <= M.shape[0]:
                  if i == 1:
                      to_the_right = M[0, i:]
                      to_the_bottom = M[i:, 0]
                  else:
                      to_the_right = M[:i, i:]
                      to_the_bottom = M[i:, :i]
                  if any(to_the_right) or any(to_the_bottom):
                      i += 1
                      continue
                  else:
                      sub_blocks.append(M[:i, :i])
                      if M.shape == M[:i, :i].shape:
                          return
                      else:
                          recurse_sub_blocks(M[i:, i:])
                          return
  
          recurse_sub_blocks(self)
          return sub_blocks
  
      def _eval_row_del(self, row):
          def entry(i, j):
              return self[i, j] if i < row else self[i + 1, j]
          return self._new(self.rows - 1, self.cols, entry)
  
      def _eval_row_insert(self, pos, other):
          entries = list(self)
          insert_pos = pos * self.cols
          entries[insert_pos:insert_pos] = list(other)
          return self._new(self.rows + other.rows, self.cols, entries)
  
      def _eval_row_join(self, other):
          cols = self.cols
  
          def entry(i, j):
              if j < cols:
                  return self[i, j]
              return other[i, j - cols]
  
          return classof(self, other)._new(self.rows, self.cols + other.cols,
                                           lambda i, j: entry(i, j))
  
      def _eval_tolist(self):
          return [list(self[i,:]) for i in range(self.rows)]
  
      def _eval_todok(self):
          dok = {}
          rows, cols = self.shape
          for i in range(rows):
              for j in range(cols):
                  val = self[i, j]
                  if val != self.zero:
                      dok[i, j] = val
          return dok
  
      def _eval_vec(self):
          rows = self.rows
  
          def entry(n, _):
              # we want to read off the columns first
              j = n // rows
              i = n - j * rows
              return self[i, j]
  
          return self._new(len(self), 1, entry)
  
      def _eval_vech(self, diagonal):
          c = self.cols
          v = []
          if diagonal:
              for j in range(c):
                  for i in range(j, c):
                      v.append(self[i, j])
          else:
              for j in range(c):
                  for i in range(j + 1, c):
                      v.append(self[i, j])
          return self._new(len(v), 1, v)
  
      def col_del(self, col):
          """Delete the specified column."""
          if col < 0:
              col += self.cols
          if not 0 <= col < self.cols:
              raise IndexError("Column {} is out of range.".format(col))
          return self._eval_col_del(col)
  
      def col_insert(self, pos, other):
          """Insert one or more columns at the given column position.
  
          Examples
          ========
  
          >>> from sympy import zeros, ones
          >>> M = zeros(3)
          >>> V = ones(3, 1)
          >>> M.col_insert(1, V)
          Matrix([
          [0, 1, 0, 0],
          [0, 1, 0, 0],
          [0, 1, 0, 0]])
  
          See Also
          ========
  
          col
          row_insert
          """
          # Allows you to build a matrix even if it is null matrix
          if not self:
              return type(self)(other)
  
          pos = as_int(pos)
  
          if pos < 0:
              pos = self.cols + pos
          if pos < 0:
              pos = 0
          elif pos > self.cols:
              pos = self.cols
  
          if self.rows != other.rows:
              raise ShapeError(
                  "`self` and `other` must have the same number of rows.")
  
          return self._eval_col_insert(pos, other)
  
      def col_join(self, other):
          """Concatenates two matrices along self's last and other's first row.
  
          Examples
          ========
  
          >>> from sympy import zeros, ones
          >>> M = zeros(3)
          >>> V = ones(1, 3)
          >>> M.col_join(V)
          Matrix([
          [0, 0, 0],
          [0, 0, 0],
          [0, 0, 0],
          [1, 1, 1]])
  
          See Also
          ========
  
          col
          row_join
          """
          # A null matrix can always be stacked (see  #10770)
          if self.rows == 0 and self.cols != other.cols:
              return self._new(0, other.cols, []).col_join(other)
  
          if self.cols != other.cols:
              raise ShapeError(
                  "`self` and `other` must have the same number of columns.")
          return self._eval_col_join(other)
  
      def col(self, j):
          """Elementary column selector.
  
          Examples
          ========
  
          >>> from sympy import eye
          >>> eye(2).col(0)
          Matrix([
          [1],
          [0]])
  
          See Also
          ========
  
          row
          sympy.matrices.dense.MutableDenseMatrix.col_op
          sympy.matrices.dense.MutableDenseMatrix.col_swap
          col_del
          col_join
          col_insert
          """
          return self[:, j]
  
      def extract(self, rowsList, colsList):
          """Return a submatrix by specifying a list of rows and columns.
          Negative indices can be given. All indices must be in the range
          -n <= i < n where n is the number of rows or columns.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(4, 3, range(12))
          >>> m
          Matrix([
          [0,  1,  2],
          [3,  4,  5],
          [6,  7,  8],
          [9, 10, 11]])
          >>> m.extract([0, 1, 3], [0, 1])
          Matrix([
          [0,  1],
          [3,  4],
          [9, 10]])
  
          Rows or columns can be repeated:
  
          >>> m.extract([0, 0, 1], [-1])
          Matrix([
          [2],
          [2],
          [5]])
  
          Every other row can be taken by using range to provide the indices:
  
          >>> m.extract(range(0, m.rows, 2), [-1])
          Matrix([
          [2],
          [8]])
  
          RowsList or colsList can also be a list of booleans, in which case
          the rows or columns corresponding to the True values will be selected:
  
          >>> m.extract([0, 1, 2, 3], [True, False, True])
          Matrix([
          [0,  2],
          [3,  5],
          [6,  8],
          [9, 11]])
          """
  
          if not is_sequence(rowsList) or not is_sequence(colsList):
              raise TypeError("rowsList and colsList must be iterable")
          # ensure rowsList and colsList are lists of integers
          if rowsList and all(isinstance(i, bool) for i in rowsList):
              rowsList = [index for index, item in enumerate(rowsList) if item]
          if colsList and all(isinstance(i, bool) for i in colsList):
              colsList = [index for index, item in enumerate(colsList) if item]
  
          # ensure everything is in range
          rowsList = [a2idx(k, self.rows) for k in rowsList]
          colsList = [a2idx(k, self.cols) for k in colsList]
  
          return self._eval_extract(rowsList, colsList)
  
      def get_diag_blocks(self):
          """Obtains the square sub-matrices on the main diagonal of a square matrix.
  
          Useful for inverting symbolic matrices or solving systems of
          linear equations which may be decoupled by having a block diagonal
          structure.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> from sympy.abc import x, y, z
          >>> A = Matrix([[1, 3, 0, 0], [y, z*z, 0, 0], [0, 0, x, 0], [0, 0, 0, 0]])
          >>> a1, a2, a3 = A.get_diag_blocks()
          >>> a1
          Matrix([
          [1,    3],
          [y, z**2]])
          >>> a2
          Matrix([[x]])
          >>> a3
          Matrix([[0]])
  
          """
          return self._eval_get_diag_blocks()
  
      @classmethod
      def hstack(cls, *args):
          """Return a matrix formed by joining args horizontally (i.e.
          by repeated application of row_join).
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix, eye
          >>> Matrix.hstack(eye(2), 2*eye(2))
          Matrix([
          [1, 0, 2, 0],
          [0, 1, 0, 2]])
          """
          if len(args) == 0:
              return cls._new()
  
          kls = type(args[0])
          return reduce(kls.row_join, args)
  
      def reshape(self, rows, cols):
          """Reshape the matrix. Total number of elements must remain the same.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(2, 3, lambda i, j: 1)
          >>> m
          Matrix([
          [1, 1, 1],
          [1, 1, 1]])
          >>> m.reshape(1, 6)
          Matrix([[1, 1, 1, 1, 1, 1]])
          >>> m.reshape(3, 2)
          Matrix([
          [1, 1],
          [1, 1],
          [1, 1]])
  
          """
          if self.rows * self.cols != rows * cols:
              raise ValueError("Invalid reshape parameters %d %d" % (rows, cols))
          return self._new(rows, cols, lambda i, j: self[i * cols + j])
  
      def row_del(self, row):
          """Delete the specified row."""
          if row < 0:
              row += self.rows
          if not 0 <= row < self.rows:
              raise IndexError("Row {} is out of range.".format(row))
  
          return self._eval_row_del(row)
  
      def row_insert(self, pos, other):
          """Insert one or more rows at the given row position.
  
          Examples
          ========
  
          >>> from sympy import zeros, ones
          >>> M = zeros(3)
          >>> V = ones(1, 3)
          >>> M.row_insert(1, V)
          Matrix([
          [0, 0, 0],
          [1, 1, 1],
          [0, 0, 0],
          [0, 0, 0]])
  
          See Also
          ========
  
          row
          col_insert
          """
          # Allows you to build a matrix even if it is null matrix
          if not self:
              return self._new(other)
  
          pos = as_int(pos)
  
          if pos < 0:
              pos = self.rows + pos
          if pos < 0:
              pos = 0
          elif pos > self.rows:
              pos = self.rows
  
          if self.cols != other.cols:
              raise ShapeError(
                  "`self` and `other` must have the same number of columns.")
  
          return self._eval_row_insert(pos, other)
  
      def row_join(self, other):
          """Concatenates two matrices along self's last and rhs's first column
  
          Examples
          ========
  
          >>> from sympy import zeros, ones
          >>> M = zeros(3)
          >>> V = ones(3, 1)
          >>> M.row_join(V)
          Matrix([
          [0, 0, 0, 1],
          [0, 0, 0, 1],
          [0, 0, 0, 1]])
  
          See Also
          ========
  
          row
          col_join
          """
          # A null matrix can always be stacked (see  #10770)
          if self.cols == 0 and self.rows != other.rows:
              return self._new(other.rows, 0, []).row_join(other)
  
          if self.rows != other.rows:
              raise ShapeError(
                  "`self` and `rhs` must have the same number of rows.")
          return self._eval_row_join(other)
  
      def diagonal(self, k=0):
          """Returns the kth diagonal of self. The main diagonal
          corresponds to `k=0`; diagonals above and below correspond to
          `k > 0` and `k < 0`, respectively. The values of `self[i, j]`
          for which `j - i = k`, are returned in order of increasing
          `i + j`, starting with `i + j = |k|`.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(3, 3, lambda i, j: j - i); m
          Matrix([
          [ 0,  1, 2],
          [-1,  0, 1],
          [-2, -1, 0]])
          >>> _.diagonal()
          Matrix([[0, 0, 0]])
          >>> m.diagonal(1)
          Matrix([[1, 1]])
          >>> m.diagonal(-2)
          Matrix([[-2]])
  
          Even though the diagonal is returned as a Matrix, the element
          retrieval can be done with a single index:
  
          >>> Matrix.diag(1, 2, 3).diagonal()[1]  # instead of [0, 1]
          2
  
          See Also
          ========
          diag - to create a diagonal matrix
          """
          rv = []
          k = as_int(k)
          r = 0 if k > 0 else -k
          c = 0 if r else k
          while True:
              if r == self.rows or c == self.cols:
                  break
              rv.append(self[r, c])
              r += 1
              c += 1
          if not rv:
              raise ValueError(filldedent('''
              The %s diagonal is out of range [%s, %s]''' % (
              k, 1 - self.rows, self.cols - 1)))
          return self._new(1, len(rv), rv)
  
      def row(self, i):
          """Elementary row selector.
  
          Examples
          ========
  
          >>> from sympy import eye
          >>> eye(2).row(0)
          Matrix([[1, 0]])
  
          See Also
          ========
  
          col
          sympy.matrices.dense.MutableDenseMatrix.row_op
          sympy.matrices.dense.MutableDenseMatrix.row_swap
          row_del
          row_join
          row_insert
          """
          return self[i, :]
  
      @property
      def shape(self):
          """The shape (dimensions) of the matrix as the 2-tuple (rows, cols).
  
          Examples
          ========
  
          >>> from sympy.matrices import zeros
          >>> M = zeros(2, 3)
          >>> M.shape
          (2, 3)
          >>> M.rows
          2
          >>> M.cols
          3
          """
          return (self.rows, self.cols)
  
      def todok(self):
          """Return the matrix as dictionary of keys.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> M = Matrix.eye(3)
          >>> M.todok()
          {(0, 0): 1, (1, 1): 1, (2, 2): 1}
          """
          return self._eval_todok()
  
      def tolist(self):
          """Return the Matrix as a nested Python list.
  
          Examples
          ========
  
          >>> from sympy import Matrix, ones
          >>> m = Matrix(3, 3, range(9))
          >>> m
          Matrix([
          [0, 1, 2],
          [3, 4, 5],
          [6, 7, 8]])
          >>> m.tolist()
          [[0, 1, 2], [3, 4, 5], [6, 7, 8]]
          >>> ones(3, 0).tolist()
          [[], [], []]
  
          When there are no rows then it will not be possible to tell how
          many columns were in the original matrix:
  
          >>> ones(0, 3).tolist()
          []
  
          """
          if not self.rows:
              return []
          if not self.cols:
              return [[] for i in range(self.rows)]
          return self._eval_tolist()
  
      def vec(self):
          """Return the Matrix converted into a one column matrix by stacking columns
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m=Matrix([[1, 3], [2, 4]])
          >>> m
          Matrix([
          [1, 3],
          [2, 4]])
          >>> m.vec()
          Matrix([
          [1],
          [2],
          [3],
          [4]])
  
          See Also
          ========
  
          vech
          """
          return self._eval_vec()
  
      def vech(self, diagonal=True, check_symmetry=True):
          """Reshapes the matrix into a column vector by stacking the
          elements in the lower triangle.
  
          Parameters
          ==========
  
          diagonal : bool, optional
              If ``True``, it includes the diagonal elements.
  
          check_symmetry : bool, optional
              If ``True``, it checks whether the matrix is symmetric.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m=Matrix([[1, 2], [2, 3]])
          >>> m
          Matrix([
          [1, 2],
          [2, 3]])
          >>> m.vech()
          Matrix([
          [1],
          [2],
          [3]])
          >>> m.vech(diagonal=False)
          Matrix([[2]])
  
          Notes
          =====
  
          This should work for symmetric matrices and ``vech`` can
          represent symmetric matrices in vector form with less size than
          ``vec``.
  
          See Also
          ========
  
          vec
          """
          if not self.is_square:
              raise NonSquareMatrixError
  
          if check_symmetry and not self.is_symmetric():
              raise ValueError("The matrix is not symmetric.")
  
          return self._eval_vech(diagonal)
  
      @classmethod
      def vstack(cls, *args):
          """Return a matrix formed by joining args vertically (i.e.
          by repeated application of col_join).
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix, eye
          >>> Matrix.vstack(eye(2), 2*eye(2))
          Matrix([
          [1, 0],
          [0, 1],
          [2, 0],
          [0, 2]])
          """
          if len(args) == 0:
              return cls._new()
  
          kls = type(args[0])
          return reduce(kls.col_join, args)
  
  
  class MatrixSpecial(MatrixRequired):
      """Construction of special matrices"""
  
      @classmethod
      def _eval_diag(cls, rows, cols, diag_dict):
          """diag_dict is a defaultdict containing
          all the entries of the diagonal matrix."""
          def entry(i, j):
              return diag_dict[(i, j)]
          return cls._new(rows, cols, entry)
  
      @classmethod
      def _eval_eye(cls, rows, cols):
          vals = [cls.zero]*(rows*cols)
          vals[::cols+1] = [cls.one]*min(rows, cols)
          return cls._new(rows, cols, vals, copy=False)
  
      @classmethod
      def _eval_jordan_block(cls, rows, cols, eigenvalue, band='upper'):
          if band == 'lower':
              def entry(i, j):
                  if i == j:
                      return eigenvalue
                  elif j + 1 == i:
                      return cls.one
                  return cls.zero
          else:
              def entry(i, j):
                  if i == j:
                      return eigenvalue
                  elif i + 1 == j:
                      return cls.one
                  return cls.zero
          return cls._new(rows, cols, entry)
  
      @classmethod
      def _eval_ones(cls, rows, cols):
          def entry(i, j):
              return cls.one
          return cls._new(rows, cols, entry)
  
      @classmethod
      def _eval_zeros(cls, rows, cols):
          return cls._new(rows, cols, [cls.zero]*(rows*cols), copy=False)
  
      @classmethod
      def _eval_wilkinson(cls, n):
          def entry(i, j):
              return cls.one if i + 1 == j else cls.zero
  
          D = cls._new(2*n + 1, 2*n + 1, entry)
  
          wminus = cls.diag([i for i in range(-n, n + 1)], unpack=True) + D + D.T
          wplus = abs(cls.diag([i for i in range(-n, n + 1)], unpack=True)) + D + D.T
  
          return wminus, wplus
  
      @classmethod
      def diag(kls, *args, strict=False, unpack=True, rows=None, cols=None, **kwargs):
          """Returns a matrix with the specified diagonal.
          If matrices are passed, a block-diagonal matrix
          is created (i.e. the "direct sum" of the matrices).
  
          kwargs
          ======
  
          rows : rows of the resulting matrix; computed if
                 not given.
  
          cols : columns of the resulting matrix; computed if
                 not given.
  
          cls : class for the resulting matrix
  
          unpack : bool which, when True (default), unpacks a single
          sequence rather than interpreting it as a Matrix.
  
          strict : bool which, when False (default), allows Matrices to
          have variable-length rows.
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> Matrix.diag(1, 2, 3)
          Matrix([
          [1, 0, 0],
          [0, 2, 0],
          [0, 0, 3]])
  
          The current default is to unpack a single sequence. If this is
          not desired, set `unpack=False` and it will be interpreted as
          a matrix.
  
          >>> Matrix.diag([1, 2, 3]) == Matrix.diag(1, 2, 3)
          True
  
          When more than one element is passed, each is interpreted as
          something to put on the diagonal. Lists are converted to
          matrices. Filling of the diagonal always continues from
          the bottom right hand corner of the previous item: this
          will create a block-diagonal matrix whether the matrices
          are square or not.
  
          >>> col = [1, 2, 3]
          >>> row = [[4, 5]]
          >>> Matrix.diag(col, row)
          Matrix([
          [1, 0, 0],
          [2, 0, 0],
          [3, 0, 0],
          [0, 4, 5]])
  
          When `unpack` is False, elements within a list need not all be
          of the same length. Setting `strict` to True would raise a
          ValueError for the following:
  
          >>> Matrix.diag([[1, 2, 3], [4, 5], [6]], unpack=False)
          Matrix([
          [1, 2, 3],
          [4, 5, 0],
          [6, 0, 0]])
  
          The type of the returned matrix can be set with the ``cls``
          keyword.
  
          >>> from sympy.matrices import ImmutableMatrix
          >>> from sympy.utilities.misc import func_name
          >>> func_name(Matrix.diag(1, cls=ImmutableMatrix))
          'ImmutableDenseMatrix'
  
          A zero dimension matrix can be used to position the start of
          the filling at the start of an arbitrary row or column:
  
          >>> from sympy import ones
          >>> r2 = ones(0, 2)
          >>> Matrix.diag(r2, 1, 2)
          Matrix([
          [0, 0, 1, 0],
          [0, 0, 0, 2]])
  
          See Also
          ========
          eye
          diagonal - to extract a diagonal
          .dense.diag
          .expressions.blockmatrix.BlockMatrix
          .sparsetools.banded - to create multi-diagonal matrices
         """
          from sympy.matrices.matrices import MatrixBase
          from sympy.matrices.dense import Matrix
          from sympy.matrices.sparse import SparseMatrix
          klass = kwargs.get('cls', kls)
          if unpack and len(args) == 1 and is_sequence(args[0]) and \
                  not isinstance(args[0], MatrixBase):
              args = args[0]
  
          # fill a default dict with the diagonal entries
          diag_entries = defaultdict(int)
          rmax = cmax = 0  # keep track of the biggest index seen
          for m in args:
              if isinstance(m, list):
                  if strict:
                      # if malformed, Matrix will raise an error
                      _ = Matrix(m)
                      r, c = _.shape
                      m = _.tolist()
                  else:
                      r, c, smat = SparseMatrix._handle_creation_inputs(m)
                      for (i, j), _ in smat.items():
                          diag_entries[(i + rmax, j + cmax)] = _
                      m = []  # to skip process below
              elif hasattr(m, 'shape'):  # a Matrix
                  # convert to list of lists
                  r, c = m.shape
                  m = m.tolist()
              else:  # in this case, we're a single value
                  diag_entries[(rmax, cmax)] = m
                  rmax += 1
                  cmax += 1
                  continue
              # process list of lists
              for i in range(len(m)):
                  for j, _ in enumerate(m[i]):
                      diag_entries[(i + rmax, j + cmax)] = _
              rmax += r
              cmax += c
          if rows is None:
              rows, cols = cols, rows
          if rows is None:
              rows, cols = rmax, cmax
          else:
              cols = rows if cols is None else cols
          if rows < rmax or cols < cmax:
              raise ValueError(filldedent('''
                  The constructed matrix is {} x {} but a size of {} x {}
                  was specified.'''.format(rmax, cmax, rows, cols)))
          return klass._eval_diag(rows, cols, diag_entries)
  
      @classmethod
      def eye(kls, rows, cols=None, **kwargs):
          """Returns an identity matrix.
  
          Args
          ====
  
          rows : rows of the matrix
          cols : cols of the matrix (if None, cols=rows)
  
          kwargs
          ======
          cls : class of the returned matrix
          """
          if cols is None:
              cols = rows
          if rows < 0 or cols < 0:
              raise ValueError("Cannot create a {} x {} matrix. "
                               "Both dimensions must be positive".format(rows, cols))
          klass = kwargs.get('cls', kls)
          rows, cols = as_int(rows), as_int(cols)
  
          return klass._eval_eye(rows, cols)
  
      @classmethod
      def jordan_block(kls, size=None, eigenvalue=None, *, band='upper', **kwargs):
          """Returns a Jordan block
  
          Parameters
          ==========
  
          size : Integer, optional
              Specifies the shape of the Jordan block matrix.
  
          eigenvalue : Number or Symbol
              Specifies the value for the main diagonal of the matrix.
  
              .. note::
                  The keyword ``eigenval`` is also specified as an alias
                  of this keyword, but it is not recommended to use.
  
                  We may deprecate the alias in later release.
  
          band : 'upper' or 'lower', optional
              Specifies the position of the off-diagonal to put `1` s on.
  
          cls : Matrix, optional
              Specifies the matrix class of the output form.
  
              If it is not specified, the class type where the method is
              being executed on will be returned.
  
          rows, cols : Integer, optional
              Specifies the shape of the Jordan block matrix. See Notes
              section for the details of how these key works.
  
              .. note::
                  This feature will be deprecated in the future.
  
  
          Returns
          =======
  
          Matrix
              A Jordan block matrix.
  
          Raises
          ======
  
          ValueError
              If insufficient arguments are given for matrix size
              specification, or no eigenvalue is given.
  
          Examples
          ========
  
          Creating a default Jordan block:
  
          >>> from sympy import Matrix
          >>> from sympy.abc import x
          >>> Matrix.jordan_block(4, x)
          Matrix([
          [x, 1, 0, 0],
          [0, x, 1, 0],
          [0, 0, x, 1],
          [0, 0, 0, x]])
  
          Creating an alternative Jordan block matrix where `1` is on
          lower off-diagonal:
  
          >>> Matrix.jordan_block(4, x, band='lower')
          Matrix([
          [x, 0, 0, 0],
          [1, x, 0, 0],
          [0, 1, x, 0],
          [0, 0, 1, x]])
  
          Creating a Jordan block with keyword arguments
  
          >>> Matrix.jordan_block(size=4, eigenvalue=x)
          Matrix([
          [x, 1, 0, 0],
          [0, x, 1, 0],
          [0, 0, x, 1],
          [0, 0, 0, x]])
  
          Notes
          =====
  
          .. note::
              This feature will be deprecated in the future.
  
          The keyword arguments ``size``, ``rows``, ``cols`` relates to
          the Jordan block size specifications.
  
          If you want to create a square Jordan block, specify either
          one of the three arguments.
  
          If you want to create a rectangular Jordan block, specify
          ``rows`` and ``cols`` individually.
  
          +--------------------------------+---------------------+
          |        Arguments Given         |     Matrix Shape    |
          +----------+----------+----------+----------+----------+
          |   size   |   rows   |   cols   |   rows   |   cols   |
          +==========+==========+==========+==========+==========+
          |   size   |         Any         |   size   |   size   |
          +----------+----------+----------+----------+----------+
          |          |        None         |     ValueError      |
          |          +----------+----------+----------+----------+
          |   None   |   rows   |   None   |   rows   |   rows   |
          |          +----------+----------+----------+----------+
          |          |   None   |   cols   |   cols   |   cols   |
          +          +----------+----------+----------+----------+
          |          |   rows   |   cols   |   rows   |   cols   |
          +----------+----------+----------+----------+----------+
  
          References
          ==========
  
          .. [1] https://en.wikipedia.org/wiki/Jordan_matrix
          """
          if 'rows' in kwargs or 'cols' in kwargs:
              SymPyDeprecationWarning(
                  feature="Keyword arguments 'rows' or 'cols'",
                  issue=16102,
                  useinstead="a more generic banded matrix constructor",
                  deprecated_since_version="1.4"
              ).warn()
  
          klass = kwargs.pop('cls', kls)
          rows = kwargs.pop('rows', None)
          cols = kwargs.pop('cols', None)
  
          eigenval = kwargs.get('eigenval', None)
          if eigenvalue is None and eigenval is None:
              raise ValueError("Must supply an eigenvalue")
          elif eigenvalue != eigenval and None not in (eigenval, eigenvalue):
              raise ValueError(
                  "Inconsistent values are given: 'eigenval'={}, "
                  "'eigenvalue'={}".format(eigenval, eigenvalue))
          else:
              if eigenval is not None:
                  eigenvalue = eigenval
  
          if (size, rows, cols) == (None, None, None):
              raise ValueError("Must supply a matrix size")
  
          if size is not None:
              rows, cols = size, size
          elif rows is not None and cols is None:
              cols = rows
          elif cols is not None and rows is None:
              rows = cols
  
          rows, cols = as_int(rows), as_int(cols)
  
          return klass._eval_jordan_block(rows, cols, eigenvalue, band)
  
      @classmethod
      def ones(kls, rows, cols=None, **kwargs):
          """Returns a matrix of ones.
  
          Args
          ====
  
          rows : rows of the matrix
          cols : cols of the matrix (if None, cols=rows)
  
          kwargs
          ======
          cls : class of the returned matrix
          """
          if cols is None:
              cols = rows
          klass = kwargs.get('cls', kls)
          rows, cols = as_int(rows), as_int(cols)
  
          return klass._eval_ones(rows, cols)
  
      @classmethod
      def zeros(kls, rows, cols=None, **kwargs):
          """Returns a matrix of zeros.
  
          Args
          ====
  
          rows : rows of the matrix
          cols : cols of the matrix (if None, cols=rows)
  
          kwargs
          ======
          cls : class of the returned matrix
          """
          if cols is None:
              cols = rows
          if rows < 0 or cols < 0:
              raise ValueError("Cannot create a {} x {} matrix. "
                               "Both dimensions must be positive".format(rows, cols))
          klass = kwargs.get('cls', kls)
          rows, cols = as_int(rows), as_int(cols)
  
          return klass._eval_zeros(rows, cols)
  
      @classmethod
      def companion(kls, poly):
          """Returns a companion matrix of a polynomial.
  
          Examples
          ========
  
          >>> from sympy import Matrix, Poly, Symbol, symbols
          >>> x = Symbol('x')
          >>> c0, c1, c2, c3, c4 = symbols('c0:5')
          >>> p = Poly(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + x**5, x)
          >>> Matrix.companion(p)
          Matrix([
          [0, 0, 0, 0, -c0],
          [1, 0, 0, 0, -c1],
          [0, 1, 0, 0, -c2],
          [0, 0, 1, 0, -c3],
          [0, 0, 0, 1, -c4]])
          """
          poly = kls._sympify(poly)
          if not isinstance(poly, Poly):
              raise ValueError("{} must be a Poly instance.".format(poly))
          if not poly.is_monic:
              raise ValueError("{} must be a monic polynomial.".format(poly))
          if not poly.is_univariate:
              raise ValueError(
                  "{} must be a univariate polynomial.".format(poly))
  
          size = poly.degree()
          if not size >= 1:
              raise ValueError(
                  "{} must have degree not less than 1.".format(poly))
  
          coeffs = poly.all_coeffs()
          def entry(i, j):
              if j == size - 1:
                  return -coeffs[-1 - i]
              elif i == j + 1:
                  return kls.one
              return kls.zero
          return kls._new(size, size, entry)
  
  
      @classmethod
      def wilkinson(kls, n, **kwargs):
          """Returns two square Wilkinson Matrix of size 2*n + 1
          $W_{2n + 1}^-, W_{2n + 1}^+ =$ Wilkinson(n)
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> wminus, wplus = Matrix.wilkinson(3)
          >>> wminus
          Matrix([
          [-3,  1,  0, 0, 0, 0, 0],
          [ 1, -2,  1, 0, 0, 0, 0],
          [ 0,  1, -1, 1, 0, 0, 0],
          [ 0,  0,  1, 0, 1, 0, 0],
          [ 0,  0,  0, 1, 1, 1, 0],
          [ 0,  0,  0, 0, 1, 2, 1],
          [ 0,  0,  0, 0, 0, 1, 3]])
          >>> wplus
          Matrix([
          [3, 1, 0, 0, 0, 0, 0],
          [1, 2, 1, 0, 0, 0, 0],
          [0, 1, 1, 1, 0, 0, 0],
          [0, 0, 1, 0, 1, 0, 0],
          [0, 0, 0, 1, 1, 1, 0],
          [0, 0, 0, 0, 1, 2, 1],
          [0, 0, 0, 0, 0, 1, 3]])
  
          References
          ==========
  
          .. [1] https://blogs.mathworks.com/cleve/2013/04/15/wilkinsons-matrices-2/
          .. [2] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Claredon Press, Oxford, 1965, 662 pp.
  
          """
          klass = kwargs.get('cls', kls)
          n = as_int(n)
          return klass._eval_wilkinson(n)
  
  class MatrixProperties(MatrixRequired):
      """Provides basic properties of a matrix."""
  
      def _eval_atoms(self, *types):
          result = set()
          for i in self:
              result.update(i.atoms(*types))
          return result
  
      def _eval_free_symbols(self):
          return set().union(*(i.free_symbols for i in self if i))
  
      def _eval_has(self, *patterns):
          return any(a.has(*patterns) for a in self)
  
      def _eval_is_anti_symmetric(self, simpfunc):
          if not all(simpfunc(self[i, j] + self[j, i]).is_zero for i in range(self.rows) for j in range(self.cols)):
              return False
          return True
  
      def _eval_is_diagonal(self):
          for i in range(self.rows):
              for j in range(self.cols):
                  if i != j and self[i, j]:
                      return False
          return True
  
      # _eval_is_hermitian is called by some general sympy
      # routines and has a different *args signature.  Make
      # sure the names don't clash by adding `_matrix_` in name.
      def _eval_is_matrix_hermitian(self, simpfunc):
          mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i].conjugate()))
          return mat.is_zero_matrix
  
      def _eval_is_Identity(self) -> FuzzyBool:
          def dirac(i, j):
              if i == j:
                  return 1
              return 0
  
          return all(self[i, j] == dirac(i, j)
                  for i in range(self.rows)
                  for j in range(self.cols))
  
      def _eval_is_lower_hessenberg(self):
          return all(self[i, j].is_zero
                     for i in range(self.rows)
                     for j in range(i + 2, self.cols))
  
      def _eval_is_lower(self):
          return all(self[i, j].is_zero
                     for i in range(self.rows)
                     for j in range(i + 1, self.cols))
  
      def _eval_is_symbolic(self):
          return self.has(Symbol)
  
      def _eval_is_symmetric(self, simpfunc):
          mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i]))
          return mat.is_zero_matrix
  
      def _eval_is_zero_matrix(self):
          if any(i.is_zero == False for i in self):
              return False
          if any(i.is_zero is None for i in self):
              return None
          return True
  
      def _eval_is_upper_hessenberg(self):
          return all(self[i, j].is_zero
                     for i in range(2, self.rows)
                     for j in range(min(self.cols, (i - 1))))
  
      def _eval_values(self):
          return [i for i in self if not i.is_zero]
  
      def _has_positive_diagonals(self):
          diagonal_entries = (self[i, i] for i in range(self.rows))
          return fuzzy_and(x.is_positive for x in diagonal_entries)
  
      def _has_nonnegative_diagonals(self):
          diagonal_entries = (self[i, i] for i in range(self.rows))
          return fuzzy_and(x.is_nonnegative for x in diagonal_entries)
  
      def atoms(self, *types):
          """Returns the atoms that form the current object.
  
          Examples
          ========
  
          >>> from sympy.abc import x, y
          >>> from sympy.matrices import Matrix
          >>> Matrix([[x]])
          Matrix([[x]])
          >>> _.atoms()
          {x}
          >>> Matrix([[x, y], [y, x]])
          Matrix([
          [x, y],
          [y, x]])
          >>> _.atoms()
          {x, y}
          """
  
          types = tuple(t if isinstance(t, type) else type(t) for t in types)
          if not types:
              types = (Atom,)
          return self._eval_atoms(*types)
  
      @property
      def free_symbols(self):
          """Returns the free symbols within the matrix.
  
          Examples
          ========
  
          >>> from sympy.abc import x
          >>> from sympy.matrices import Matrix
          >>> Matrix([[x], [1]]).free_symbols
          {x}
          """
          return self._eval_free_symbols()
  
      def has(self, *patterns):
          """Test whether any subexpression matches any of the patterns.
  
          Examples
          ========
  
          >>> from sympy import Matrix, SparseMatrix, Float
          >>> from sympy.abc import x, y
          >>> A = Matrix(((1, x), (0.2, 3)))
          >>> B = SparseMatrix(((1, x), (0.2, 3)))
          >>> A.has(x)
          True
          >>> A.has(y)
          False
          >>> A.has(Float)
          True
          >>> B.has(x)
          True
          >>> B.has(y)
          False
          >>> B.has(Float)
          True
          """
          return self._eval_has(*patterns)
  
      def is_anti_symmetric(self, simplify=True):
          """Check if matrix M is an antisymmetric matrix,
          that is, M is a square matrix with all M[i, j] == -M[j, i].
  
          When ``simplify=True`` (default), the sum M[i, j] + M[j, i] is
          simplified before testing to see if it is zero. By default,
          the SymPy simplify function is used. To use a custom function
          set simplify to a function that accepts a single argument which
          returns a simplified expression. To skip simplification, set
          simplify to False but note that although this will be faster,
          it may induce false negatives.
  
          Examples
          ========
  
          >>> from sympy import Matrix, symbols
          >>> m = Matrix(2, 2, [0, 1, -1, 0])
          >>> m
          Matrix([
          [ 0, 1],
          [-1, 0]])
          >>> m.is_anti_symmetric()
          True
          >>> x, y = symbols('x y')
          >>> m = Matrix(2, 3, [0, 0, x, -y, 0, 0])
          >>> m
          Matrix([
          [ 0, 0, x],
          [-y, 0, 0]])
          >>> m.is_anti_symmetric()
          False
  
          >>> from sympy.abc import x, y
          >>> m = Matrix(3, 3, [0, x**2 + 2*x + 1, y,
          ...                   -(x + 1)**2 , 0, x*y,
          ...                   -y, -x*y, 0])
  
          Simplification of matrix elements is done by default so even
          though two elements which should be equal and opposite wouldn't
          pass an equality test, the matrix is still reported as
          anti-symmetric:
  
          >>> m[0, 1] == -m[1, 0]
          False
          >>> m.is_anti_symmetric()
          True
  
          If 'simplify=False' is used for the case when a Matrix is already
          simplified, this will speed things up. Here, we see that without
          simplification the matrix does not appear anti-symmetric:
  
          >>> m.is_anti_symmetric(simplify=False)
          False
  
          But if the matrix were already expanded, then it would appear
          anti-symmetric and simplification in the is_anti_symmetric routine
          is not needed:
  
          >>> m = m.expand()
          >>> m.is_anti_symmetric(simplify=False)
          True
          """
          # accept custom simplification
          simpfunc = simplify
          if not isfunction(simplify):
              simpfunc = _simplify if simplify else lambda x: x
  
          if not self.is_square:
              return False
          return self._eval_is_anti_symmetric(simpfunc)
  
      def is_diagonal(self):
          """Check if matrix is diagonal,
          that is matrix in which the entries outside the main diagonal are all zero.
  
          Examples
          ========
  
          >>> from sympy import Matrix, diag
          >>> m = Matrix(2, 2, [1, 0, 0, 2])
          >>> m
          Matrix([
          [1, 0],
          [0, 2]])
          >>> m.is_diagonal()
          True
  
          >>> m = Matrix(2, 2, [1, 1, 0, 2])
          >>> m
          Matrix([
          [1, 1],
          [0, 2]])
          >>> m.is_diagonal()
          False
  
          >>> m = diag(1, 2, 3)
          >>> m
          Matrix([
          [1, 0, 0],
          [0, 2, 0],
          [0, 0, 3]])
          >>> m.is_diagonal()
          True
  
          See Also
          ========
  
          is_lower
          is_upper
          sympy.matrices.matrices.MatrixEigen.is_diagonalizable
          diagonalize
          """
          return self._eval_is_diagonal()
  
      @property
      def is_weakly_diagonally_dominant(self):
          r"""Tests if the matrix is row weakly diagonally dominant.
  
          Explanation
          ===========
  
          A $n, n$ matrix $A$ is row weakly diagonally dominant if
  
          .. math::
              \left|A_{i, i}\right| \ge \sum_{j = 0, j \neq i}^{n-1}
              \left|A_{i, j}\right| \quad {\text{for all }}
              i \in \{ 0, ..., n-1 \}
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> A = Matrix([[3, -2, 1], [1, -3, 2], [-1, 2, 4]])
          >>> A.is_weakly_diagonally_dominant
          True
  
          >>> A = Matrix([[-2, 2, 1], [1, 3, 2], [1, -2, 0]])
          >>> A.is_weakly_diagonally_dominant
          False
  
          >>> A = Matrix([[-4, 2, 1], [1, 6, 2], [1, -2, 5]])
          >>> A.is_weakly_diagonally_dominant
          True
  
          Notes
          =====
  
          If you want to test whether a matrix is column diagonally
          dominant, you can apply the test after transposing the matrix.
          """
          if not self.is_square:
              return False
  
          rows, cols = self.shape
  
          def test_row(i):
              summation = self.zero
              for j in range(cols):
                  if i != j:
                      summation += Abs(self[i, j])
              return (Abs(self[i, i]) - summation).is_nonnegative
  
          return fuzzy_and(test_row(i) for i in range(rows))
  
      @property
      def is_strongly_diagonally_dominant(self):
          r"""Tests if the matrix is row strongly diagonally dominant.
  
          Explanation
          ===========
  
          A $n, n$ matrix $A$ is row strongly diagonally dominant if
  
          .. math::
              \left|A_{i, i}\right| > \sum_{j = 0, j \neq i}^{n-1}
              \left|A_{i, j}\right| \quad {\text{for all }}
              i \in \{ 0, ..., n-1 \}
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> A = Matrix([[3, -2, 1], [1, -3, 2], [-1, 2, 4]])
          >>> A.is_strongly_diagonally_dominant
          False
  
          >>> A = Matrix([[-2, 2, 1], [1, 3, 2], [1, -2, 0]])
          >>> A.is_strongly_diagonally_dominant
          False
  
          >>> A = Matrix([[-4, 2, 1], [1, 6, 2], [1, -2, 5]])
          >>> A.is_strongly_diagonally_dominant
          True
  
          Notes
          =====
  
          If you want to test whether a matrix is column diagonally
          dominant, you can apply the test after transposing the matrix.
          """
          if not self.is_square:
              return False
  
          rows, cols = self.shape
  
          def test_row(i):
              summation = self.zero
              for j in range(cols):
                  if i != j:
                      summation += Abs(self[i, j])
              return (Abs(self[i, i]) - summation).is_positive
  
          return fuzzy_and(test_row(i) for i in range(rows))
  
      @property
      def is_hermitian(self):
          """Checks if the matrix is Hermitian.
  
          In a Hermitian matrix element i,j is the complex conjugate of
          element j,i.
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> from sympy import I
          >>> from sympy.abc import x
          >>> a = Matrix([[1, I], [-I, 1]])
          >>> a
          Matrix([
          [ 1, I],
          [-I, 1]])
          >>> a.is_hermitian
          True
          >>> a[0, 0] = 2*I
          >>> a.is_hermitian
          False
          >>> a[0, 0] = x
          >>> a.is_hermitian
          >>> a[0, 1] = a[1, 0]*I
          >>> a.is_hermitian
          False
          """
          if not self.is_square:
              return False
  
          return self._eval_is_matrix_hermitian(_simplify)
  
      @property
      def is_Identity(self) -> FuzzyBool:
          if not self.is_square:
              return False
          return self._eval_is_Identity()
  
      @property
      def is_lower_hessenberg(self):
          r"""Checks if the matrix is in the lower-Hessenberg form.
  
          The lower hessenberg matrix has zero entries
          above the first superdiagonal.
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> a = Matrix([[1, 2, 0, 0], [5, 2, 3, 0], [3, 4, 3, 7], [5, 6, 1, 1]])
          >>> a
          Matrix([
          [1, 2, 0, 0],
          [5, 2, 3, 0],
          [3, 4, 3, 7],
          [5, 6, 1, 1]])
          >>> a.is_lower_hessenberg
          True
  
          See Also
          ========
  
          is_upper_hessenberg
          is_lower
          """
          return self._eval_is_lower_hessenberg()
  
      @property
      def is_lower(self):
          """Check if matrix is a lower triangular matrix. True can be returned
          even if the matrix is not square.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(2, 2, [1, 0, 0, 1])
          >>> m
          Matrix([
          [1, 0],
          [0, 1]])
          >>> m.is_lower
          True
  
          >>> m = Matrix(4, 3, [0, 0, 0, 2, 0, 0, 1, 4 , 0, 6, 6, 5])
          >>> m
          Matrix([
          [0, 0, 0],
          [2, 0, 0],
          [1, 4, 0],
          [6, 6, 5]])
          >>> m.is_lower
          True
  
          >>> from sympy.abc import x, y
          >>> m = Matrix(2, 2, [x**2 + y, y**2 + x, 0, x + y])
          >>> m
          Matrix([
          [x**2 + y, x + y**2],
          [       0,    x + y]])
          >>> m.is_lower
          False
  
          See Also
          ========
  
          is_upper
          is_diagonal
          is_lower_hessenberg
          """
          return self._eval_is_lower()
  
      @property
      def is_square(self):
          """Checks if a matrix is square.
  
          A matrix is square if the number of rows equals the number of columns.
          The empty matrix is square by definition, since the number of rows and
          the number of columns are both zero.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> a = Matrix([[1, 2, 3], [4, 5, 6]])
          >>> b = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
          >>> c = Matrix([])
          >>> a.is_square
          False
          >>> b.is_square
          True
          >>> c.is_square
          True
          """
          return self.rows == self.cols
  
      def is_symbolic(self):
          """Checks if any elements contain Symbols.
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> from sympy.abc import x, y
          >>> M = Matrix([[x, y], [1, 0]])
          >>> M.is_symbolic()
          True
  
          """
          return self._eval_is_symbolic()
  
      def is_symmetric(self, simplify=True):
          """Check if matrix is symmetric matrix,
          that is square matrix and is equal to its transpose.
  
          By default, simplifications occur before testing symmetry.
          They can be skipped using 'simplify=False'; while speeding things a bit,
          this may however induce false negatives.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(2, 2, [0, 1, 1, 2])
          >>> m
          Matrix([
          [0, 1],
          [1, 2]])
          >>> m.is_symmetric()
          True
  
          >>> m = Matrix(2, 2, [0, 1, 2, 0])
          >>> m
          Matrix([
          [0, 1],
          [2, 0]])
          >>> m.is_symmetric()
          False
  
          >>> m = Matrix(2, 3, [0, 0, 0, 0, 0, 0])
          >>> m
          Matrix([
          [0, 0, 0],
          [0, 0, 0]])
          >>> m.is_symmetric()
          False
  
          >>> from sympy.abc import x, y
          >>> m = Matrix(3, 3, [1, x**2 + 2*x + 1, y, (x + 1)**2 , 2, 0, y, 0, 3])
          >>> m
          Matrix([
          [         1, x**2 + 2*x + 1, y],
          [(x + 1)**2,              2, 0],
          [         y,              0, 3]])
          >>> m.is_symmetric()
          True
  
          If the matrix is already simplified, you may speed-up is_symmetric()
          test by using 'simplify=False'.
  
          >>> bool(m.is_symmetric(simplify=False))
          False
          >>> m1 = m.expand()
          >>> m1.is_symmetric(simplify=False)
          True
          """
          simpfunc = simplify
          if not isfunction(simplify):
              simpfunc = _simplify if simplify else lambda x: x
  
          if not self.is_square:
              return False
  
          return self._eval_is_symmetric(simpfunc)
  
      @property
      def is_upper_hessenberg(self):
          """Checks if the matrix is the upper-Hessenberg form.
  
          The upper hessenberg matrix has zero entries
          below the first subdiagonal.
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> a = Matrix([[1, 4, 2, 3], [3, 4, 1, 7], [0, 2, 3, 4], [0, 0, 1, 3]])
          >>> a
          Matrix([
          [1, 4, 2, 3],
          [3, 4, 1, 7],
          [0, 2, 3, 4],
          [0, 0, 1, 3]])
          >>> a.is_upper_hessenberg
          True
  
          See Also
          ========
  
          is_lower_hessenberg
          is_upper
          """
          return self._eval_is_upper_hessenberg()
  
      @property
      def is_upper(self):
          """Check if matrix is an upper triangular matrix. True can be returned
          even if the matrix is not square.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(2, 2, [1, 0, 0, 1])
          >>> m
          Matrix([
          [1, 0],
          [0, 1]])
          >>> m.is_upper
          True
  
          >>> m = Matrix(4, 3, [5, 1, 9, 0, 4 , 6, 0, 0, 5, 0, 0, 0])
          >>> m
          Matrix([
          [5, 1, 9],
          [0, 4, 6],
          [0, 0, 5],
          [0, 0, 0]])
          >>> m.is_upper
          True
  
          >>> m = Matrix(2, 3, [4, 2, 5, 6, 1, 1])
          >>> m
          Matrix([
          [4, 2, 5],
          [6, 1, 1]])
          >>> m.is_upper
          False
  
          See Also
          ========
  
          is_lower
          is_diagonal
          is_upper_hessenberg
          """
          return all(self[i, j].is_zero
                     for i in range(1, self.rows)
                     for j in range(min(i, self.cols)))
  
      @property
      def is_zero_matrix(self):
          """Checks if a matrix is a zero matrix.
  
          A matrix is zero if every element is zero.  A matrix need not be square
          to be considered zero.  The empty matrix is zero by the principle of
          vacuous truth.  For a matrix that may or may not be zero (e.g.
          contains a symbol), this will be None
  
          Examples
          ========
  
          >>> from sympy import Matrix, zeros
          >>> from sympy.abc import x
          >>> a = Matrix([[0, 0], [0, 0]])
          >>> b = zeros(3, 4)
          >>> c = Matrix([[0, 1], [0, 0]])
          >>> d = Matrix([])
          >>> e = Matrix([[x, 0], [0, 0]])
          >>> a.is_zero_matrix
          True
          >>> b.is_zero_matrix
          True
          >>> c.is_zero_matrix
          False
          >>> d.is_zero_matrix
          True
          >>> e.is_zero_matrix
          """
          return self._eval_is_zero_matrix()
  
      def values(self):
          """Return non-zero values of self."""
          return self._eval_values()
  
  
  class MatrixOperations(MatrixRequired):
      """Provides basic matrix shape and elementwise
      operations.  Should not be instantiated directly."""
  
      def _eval_adjoint(self):
          return self.transpose().conjugate()
  
      def _eval_applyfunc(self, f):
          out = self._new(self.rows, self.cols, [f(x) for x in self])
          return out
  
      def _eval_as_real_imag(self):  # type: ignore
          from sympy.functions.elementary.complexes import re, im
  
          return (self.applyfunc(re), self.applyfunc(im))
  
      def _eval_conjugate(self):
          return self.applyfunc(lambda x: x.conjugate())
  
      def _eval_permute_cols(self, perm):
          # apply the permutation to a list
          mapping = list(perm)
  
          def entry(i, j):
              return self[i, mapping[j]]
  
          return self._new(self.rows, self.cols, entry)
  
      def _eval_permute_rows(self, perm):
          # apply the permutation to a list
          mapping = list(perm)
  
          def entry(i, j):
              return self[mapping[i], j]
  
          return self._new(self.rows, self.cols, entry)
  
      def _eval_trace(self):
          return sum(self[i, i] for i in range(self.rows))
  
      def _eval_transpose(self):
          return self._new(self.cols, self.rows, lambda i, j: self[j, i])
  
      def adjoint(self):
          """Conjugate transpose or Hermitian conjugation."""
          return self._eval_adjoint()
  
      def applyfunc(self, f):
          """Apply a function to each element of the matrix.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> m = Matrix(2, 2, lambda i, j: i*2+j)
          >>> m
          Matrix([
          [0, 1],
          [2, 3]])
          >>> m.applyfunc(lambda i: 2*i)
          Matrix([
          [0, 2],
          [4, 6]])
  
          """
          if not callable(f):
              raise TypeError("`f` must be callable.")
  
          return self._eval_applyfunc(f)
  
      def as_real_imag(self, deep=True, **hints):
          """Returns a tuple containing the (real, imaginary) part of matrix."""
          # XXX: Ignoring deep and hints...
          return self._eval_as_real_imag()
  
      def conjugate(self):
          """Return the by-element conjugation.
  
          Examples
          ========
  
          >>> from sympy.matrices import SparseMatrix
          >>> from sympy import I
          >>> a = SparseMatrix(((1, 2 + I), (3, 4), (I, -I)))
          >>> a
          Matrix([
          [1, 2 + I],
          [3,     4],
          [I,    -I]])
          >>> a.C
          Matrix([
          [ 1, 2 - I],
          [ 3,     4],
          [-I,     I]])
  
          See Also
          ========
  
          transpose: Matrix transposition
          H: Hermite conjugation
          sympy.matrices.matrices.MatrixBase.D: Dirac conjugation
          """
          return self._eval_conjugate()
  
      def doit(self, **kwargs):
          return self.applyfunc(lambda x: x.doit())
  
      def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None, verbose=False):
          """Apply evalf() to each element of self."""
          options = {'subs':subs, 'maxn':maxn, 'chop':chop, 'strict':strict,
                  'quad':quad, 'verbose':verbose}
          return self.applyfunc(lambda i: i.evalf(n, **options))
  
      def expand(self, deep=True, modulus=None, power_base=True, power_exp=True,
                 mul=True, log=True, multinomial=True, basic=True, **hints):
          """Apply core.function.expand to each entry of the matrix.
  
          Examples
          ========
  
          >>> from sympy.abc import x
          >>> from sympy.matrices import Matrix
          >>> Matrix(1, 1, [x*(x+1)])
          Matrix([[x*(x + 1)]])
          >>> _.expand()
          Matrix([[x**2 + x]])
  
          """
          return self.applyfunc(lambda x: x.expand(
              deep, modulus, power_base, power_exp, mul, log, multinomial, basic,
              **hints))
  
      @property
      def H(self):
          """Return Hermite conjugate.
  
          Examples
          ========
  
          >>> from sympy import Matrix, I
          >>> m = Matrix((0, 1 + I, 2, 3))
          >>> m
          Matrix([
          [    0],
          [1 + I],
          [    2],
          [    3]])
          >>> m.H
          Matrix([[0, 1 - I, 2, 3]])
  
          See Also
          ========
  
          conjugate: By-element conjugation
          sympy.matrices.matrices.MatrixBase.D: Dirac conjugation
          """
          return self.T.C
  
      def permute(self, perm, orientation='rows', direction='forward'):
          r"""Permute the rows or columns of a matrix by the given list of
          swaps.
  
          Parameters
          ==========
  
          perm : Permutation, list, or list of lists
              A representation for the permutation.
  
              If it is ``Permutation``, it is used directly with some
              resizing with respect to the matrix size.
  
              If it is specified as list of lists,
              (e.g., ``[[0, 1], [0, 2]]``), then the permutation is formed
              from applying the product of cycles. The direction how the
              cyclic product is applied is described in below.
  
              If it is specified as a list, the list should represent
              an array form of a permutation. (e.g., ``[1, 2, 0]``) which
              would would form the swapping function
              `0 \mapsto 1, 1 \mapsto 2, 2\mapsto 0`.
  
          orientation : 'rows', 'cols'
              A flag to control whether to permute the rows or the columns
  
          direction : 'forward', 'backward'
              A flag to control whether to apply the permutations from
              the start of the list first, or from the back of the list
              first.
  
              For example, if the permutation specification is
              ``[[0, 1], [0, 2]]``,
  
              If the flag is set to ``'forward'``, the cycle would be
              formed as `0 \mapsto 2, 2 \mapsto 1, 1 \mapsto 0`.
  
              If the flag is set to ``'backward'``, the cycle would be
              formed as `0 \mapsto 1, 1 \mapsto 2, 2 \mapsto 0`.
  
              If the argument ``perm`` is not in a form of list of lists,
              this flag takes no effect.
  
          Examples
          ========
  
          >>> from sympy.matrices import eye
          >>> M = eye(3)
          >>> M.permute([[0, 1], [0, 2]], orientation='rows', direction='forward')
          Matrix([
          [0, 0, 1],
          [1, 0, 0],
          [0, 1, 0]])
  
          >>> from sympy.matrices import eye
          >>> M = eye(3)
          >>> M.permute([[0, 1], [0, 2]], orientation='rows', direction='backward')
          Matrix([
          [0, 1, 0],
          [0, 0, 1],
          [1, 0, 0]])
  
          Notes
          =====
  
          If a bijective function
          `\sigma : \mathbb{N}_0 \rightarrow \mathbb{N}_0` denotes the
          permutation.
  
          If the matrix `A` is the matrix to permute, represented as
          a horizontal or a vertical stack of vectors:
  
          .. math::
              A =
              \begin{bmatrix}
              a_0 \\ a_1 \\ \vdots \\ a_{n-1}
              \end{bmatrix} =
              \begin{bmatrix}
              \alpha_0 & \alpha_1 & \cdots & \alpha_{n-1}
              \end{bmatrix}
  
          If the matrix `B` is the result, the permutation of matrix rows
          is defined as:
  
          .. math::
              B := \begin{bmatrix}
              a_{\sigma(0)} \\ a_{\sigma(1)} \\ \vdots \\ a_{\sigma(n-1)}
              \end{bmatrix}
  
          And the permutation of matrix columns is defined as:
  
          .. math::
              B := \begin{bmatrix}
              \alpha_{\sigma(0)} & \alpha_{\sigma(1)} &
              \cdots & \alpha_{\sigma(n-1)}
              \end{bmatrix}
          """
          from sympy.combinatorics import Permutation
  
          # allow british variants and `columns`
          if direction == 'forwards':
              direction = 'forward'
          if direction == 'backwards':
              direction = 'backward'
          if orientation == 'columns':
              orientation = 'cols'
  
          if direction not in ('forward', 'backward'):
              raise TypeError("direction='{}' is an invalid kwarg. "
                              "Try 'forward' or 'backward'".format(direction))
          if orientation not in ('rows', 'cols'):
              raise TypeError("orientation='{}' is an invalid kwarg. "
                              "Try 'rows' or 'cols'".format(orientation))
  
          if not isinstance(perm, (Permutation, Iterable)):
              raise ValueError(
                  "{} must be a list, a list of lists, "
                  "or a SymPy permutation object.".format(perm))
  
          # ensure all swaps are in range
          max_index = self.rows if orientation == 'rows' else self.cols
          if not all(0 <= t <= max_index for t in flatten(list(perm))):
              raise IndexError("`swap` indices out of range.")
  
          if perm and not isinstance(perm, Permutation) and \
              isinstance(perm[0], Iterable):
              if direction == 'forward':
                  perm = list(reversed(perm))
              perm = Permutation(perm, size=max_index+1)
          else:
              perm = Permutation(perm, size=max_index+1)
  
          if orientation == 'rows':
              return self._eval_permute_rows(perm)
          if orientation == 'cols':
              return self._eval_permute_cols(perm)
  
      def permute_cols(self, swaps, direction='forward'):
          """Alias for
          ``self.permute(swaps, orientation='cols', direction=direction)``
  
          See Also
          ========
  
          permute
          """
          return self.permute(swaps, orientation='cols', direction=direction)
  
      def permute_rows(self, swaps, direction='forward'):
          """Alias for
          ``self.permute(swaps, orientation='rows', direction=direction)``
  
          See Also
          ========
  
          permute
          """
          return self.permute(swaps, orientation='rows', direction=direction)
  
      def refine(self, assumptions=True):
          """Apply refine to each element of the matrix.
  
          Examples
          ========
  
          >>> from sympy import Symbol, Matrix, Abs, sqrt, Q
          >>> x = Symbol('x')
          >>> Matrix([[Abs(x)**2, sqrt(x**2)],[sqrt(x**2), Abs(x)**2]])
          Matrix([
          [ Abs(x)**2, sqrt(x**2)],
          [sqrt(x**2),  Abs(x)**2]])
          >>> _.refine(Q.real(x))
          Matrix([
          [  x**2, Abs(x)],
          [Abs(x),   x**2]])
  
          """
          return self.applyfunc(lambda x: refine(x, assumptions))
  
      def replace(self, F, G, map=False, simultaneous=True, exact=None):
          """Replaces Function F in Matrix entries with Function G.
  
          Examples
          ========
  
          >>> from sympy import symbols, Function, Matrix
          >>> F, G = symbols('F, G', cls=Function)
          >>> M = Matrix(2, 2, lambda i, j: F(i+j)) ; M
          Matrix([
          [F(0), F(1)],
          [F(1), F(2)]])
          >>> N = M.replace(F,G)
          >>> N
          Matrix([
          [G(0), G(1)],
          [G(1), G(2)]])
          """
          return self.applyfunc(
              lambda x: x.replace(F, G, map=map, simultaneous=simultaneous, exact=exact))
  
      def rot90(self, k=1):
          """Rotates Matrix by 90 degrees
  
          Parameters
          ==========
  
          k : int
              Specifies how many times the matrix is rotated by 90 degrees
              (clockwise when positive, counter-clockwise when negative).
  
          Examples
          ========
  
          >>> from sympy import Matrix, symbols
          >>> A = Matrix(2, 2, symbols('a:d'))
          >>> A
          Matrix([
          [a, b],
          [c, d]])
  
          Rotating the matrix clockwise one time:
  
          >>> A.rot90(1)
          Matrix([
          [c, a],
          [d, b]])
  
          Rotating the matrix anticlockwise two times:
  
          >>> A.rot90(-2)
          Matrix([
          [d, c],
          [b, a]])
          """
  
          mod = k%4
          if mod == 0:
              return self
          if mod == 1:
              return self[::-1, ::].T
          if mod == 2:
              return self[::-1, ::-1]
          if mod == 3:
              return self[::, ::-1].T
  
      def simplify(self, **kwargs):
          """Apply simplify to each element of the matrix.
  
          Examples
          ========
  
          >>> from sympy.abc import x, y
          >>> from sympy import sin, cos
          >>> from sympy.matrices import SparseMatrix
          >>> SparseMatrix(1, 1, [x*sin(y)**2 + x*cos(y)**2])
          Matrix([[x*sin(y)**2 + x*cos(y)**2]])
          >>> _.simplify()
          Matrix([[x]])
          """
          return self.applyfunc(lambda x: x.simplify(**kwargs))
  
      def subs(self, *args, **kwargs):  # should mirror core.basic.subs
          """Return a new matrix with subs applied to each entry.
  
          Examples
          ========
  
          >>> from sympy.abc import x, y
          >>> from sympy.matrices import SparseMatrix, Matrix
          >>> SparseMatrix(1, 1, [x])
          Matrix([[x]])
          >>> _.subs(x, y)
          Matrix([[y]])
          >>> Matrix(_).subs(y, x)
          Matrix([[x]])
          """
  
          if len(args) == 1 and  not isinstance(args[0], (dict, set)) and iter(args[0]) and not is_sequence(args[0]):
              args = (list(args[0]),)
  
          return self.applyfunc(lambda x: x.subs(*args, **kwargs))
  
      def trace(self):
          """
          Returns the trace of a square matrix i.e. the sum of the
          diagonal elements.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> A = Matrix(2, 2, [1, 2, 3, 4])
          >>> A.trace()
          5
  
          """
          if self.rows != self.cols:
              raise NonSquareMatrixError()
          return self._eval_trace()
  
      def transpose(self):
          """
          Returns the transpose of the matrix.
  
          Examples
          ========
  
          >>> from sympy import Matrix
          >>> A = Matrix(2, 2, [1, 2, 3, 4])
          >>> A.transpose()
          Matrix([
          [1, 3],
          [2, 4]])
  
          >>> from sympy import Matrix, I
          >>> m=Matrix(((1, 2+I), (3, 4)))
          >>> m
          Matrix([
          [1, 2 + I],
          [3,     4]])
          >>> m.transpose()
          Matrix([
          [    1, 3],
          [2 + I, 4]])
          >>> m.T == m.transpose()
          True
  
          See Also
          ========
  
          conjugate: By-element conjugation
  
          """
          return self._eval_transpose()
  
      @property
      def T(self):
          '''Matrix transposition'''
          return self.transpose()
  
      @property
      def C(self):
          '''By-element conjugation'''
          return self.conjugate()
  
      def n(self, *args, **kwargs):
          """Apply evalf() to each element of self."""
          return self.evalf(*args, **kwargs)
  
      def xreplace(self, rule):  # should mirror core.basic.xreplace
          """Return a new matrix with xreplace applied to each entry.
  
          Examples
          ========
  
          >>> from sympy.abc import x, y
          >>> from sympy.matrices import SparseMatrix, Matrix
          >>> SparseMatrix(1, 1, [x])
          Matrix([[x]])
          >>> _.xreplace({x: y})
          Matrix([[y]])
          >>> Matrix(_).xreplace({y: x})
          Matrix([[x]])
          """
          return self.applyfunc(lambda x: x.xreplace(rule))
  
      def _eval_simplify(self, **kwargs):
          # XXX: We can't use self.simplify here as mutable subclasses will
          # override simplify and have it return None
          return MatrixOperations.simplify(self, **kwargs)
  
      def _eval_trigsimp(self, **opts):
          from sympy.simplify import trigsimp
          return self.applyfunc(lambda x: trigsimp(x, **opts))
  
      def upper_triangular(self, k=0):
          """returns the elements on and above the kth diagonal of a matrix.
          If k is not specified then simply returns upper-triangular portion
          of a matrix
  
          Examples
          ========
  
          >>> from sympy import ones
          >>> A = ones(4)
          >>> A.upper_triangular()
          Matrix([
          [1, 1, 1, 1],
          [0, 1, 1, 1],
          [0, 0, 1, 1],
          [0, 0, 0, 1]])
  
          >>> A.upper_triangular(2)
          Matrix([
          [0, 0, 1, 1],
          [0, 0, 0, 1],
          [0, 0, 0, 0],
          [0, 0, 0, 0]])
  
          >>> A.upper_triangular(-1)
          Matrix([
          [1, 1, 1, 1],
          [1, 1, 1, 1],
          [0, 1, 1, 1],
          [0, 0, 1, 1]])
  
          """
  
          def entry(i, j):
              return self[i, j] if i + k <= j else self.zero
  
          return self._new(self.rows, self.cols, entry)
  
  
      def lower_triangular(self, k=0):
          """returns the elements on and below the kth diagonal of a matrix.
          If k is not specified then simply returns lower-triangular portion
          of a matrix
  
          Examples
          ========
  
          >>> from sympy import ones
          >>> A = ones(4)
          >>> A.lower_triangular()
          Matrix([
          [1, 0, 0, 0],
          [1, 1, 0, 0],
          [1, 1, 1, 0],
          [1, 1, 1, 1]])
  
          >>> A.lower_triangular(-2)
          Matrix([
          [0, 0, 0, 0],
          [0, 0, 0, 0],
          [1, 0, 0, 0],
          [1, 1, 0, 0]])
  
          >>> A.lower_triangular(1)
          Matrix([
          [1, 1, 0, 0],
          [1, 1, 1, 0],
          [1, 1, 1, 1],
          [1, 1, 1, 1]])
  
          """
  
          def entry(i, j):
              return self[i, j] if i + k >= j else self.zero
  
          return self._new(self.rows, self.cols, entry)
  
  
  
  class MatrixArithmetic(MatrixRequired):
      """Provides basic matrix arithmetic operations.
      Should not be instantiated directly."""
  
      _op_priority = 10.01
  
      def _eval_Abs(self):
          return self._new(self.rows, self.cols, lambda i, j: Abs(self[i, j]))
  
      def _eval_add(self, other):
          return self._new(self.rows, self.cols,
                           lambda i, j: self[i, j] + other[i, j])
  
      def _eval_matrix_mul(self, other):
          def entry(i, j):
              vec = [self[i,k]*other[k,j] for k in range(self.cols)]
              try:
                  return Add(*vec)
              except (TypeError, SympifyError):
                  # Some matrices don't work with `sum` or `Add`
                  # They don't work with `sum` because `sum` tries to add `0`
                  # Fall back to a safe way to multiply if the `Add` fails.
                  return reduce(lambda a, b: a + b, vec)
  
          return self._new(self.rows, other.cols, entry)
  
      def _eval_matrix_mul_elementwise(self, other):
          return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other[i,j])
  
      def _eval_matrix_rmul(self, other):
          def entry(i, j):
              return sum(other[i,k]*self[k,j] for k in range(other.cols))
          return self._new(other.rows, self.cols, entry)
  
      def _eval_pow_by_recursion(self, num):
          if num == 1:
              return self
  
          if num % 2 == 1:
              a, b = self, self._eval_pow_by_recursion(num - 1)
          else:
              a = b = self._eval_pow_by_recursion(num // 2)
  
          return a.multiply(b)
  
      def _eval_pow_by_cayley(self, exp):
          from sympy.discrete.recurrences import linrec_coeffs
          row = self.shape[0]
          p = self.charpoly()
  
          coeffs = (-p).all_coeffs()[1:]
          coeffs = linrec_coeffs(coeffs, exp)
          new_mat = self.eye(row)
          ans = self.zeros(row)
  
          for i in range(row):
              ans += coeffs[i]*new_mat
              new_mat *= self
  
          return ans
  
      def _eval_pow_by_recursion_dotprodsimp(self, num, prevsimp=None):
          if prevsimp is None:
              prevsimp = [True]*len(self)
  
          if num == 1:
              return self
  
          if num % 2 == 1:
              a, b = self, self._eval_pow_by_recursion_dotprodsimp(num - 1,
                      prevsimp=prevsimp)
          else:
              a = b = self._eval_pow_by_recursion_dotprodsimp(num // 2,
                      prevsimp=prevsimp)
  
          m     = a.multiply(b, dotprodsimp=False)
          lenm  = len(m)
          elems = [None]*lenm
  
          for i in range(lenm):
              if prevsimp[i]:
                  elems[i], prevsimp[i] = _dotprodsimp(m[i], withsimp=True)
              else:
                  elems[i] = m[i]
  
          return m._new(m.rows, m.cols, elems)
  
      def _eval_scalar_mul(self, other):
          return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other)
  
      def _eval_scalar_rmul(self, other):
          return self._new(self.rows, self.cols, lambda i, j: other*self[i,j])
  
      def _eval_Mod(self, other):
          from sympy import Mod
          return self._new(self.rows, self.cols, lambda i, j: Mod(self[i, j], other))
  
      # python arithmetic functions
      def __abs__(self):
          """Returns a new matrix with entry-wise absolute values."""
          return self._eval_Abs()
  
      @call_highest_priority('__radd__')
      def __add__(self, other):
          """Return self + other, raising ShapeError if shapes don't match."""
          if isinstance(other, NDimArray): # Matrix and array addition is currently not implemented
              return NotImplemented
          other = _matrixify(other)
          # matrix-like objects can have shapes.  This is
          # our first sanity check.
          if hasattr(other, 'shape'):
              if self.shape != other.shape:
                  raise ShapeError("Matrix size mismatch: %s + %s" % (
                      self.shape, other.shape))
  
          # honest sympy matrices defer to their class's routine
          if getattr(other, 'is_Matrix', False):
              # call the highest-priority class's _eval_add
              a, b = self, other
              if a.__class__ != classof(a, b):
                  b, a = a, b
              return a._eval_add(b)
          # Matrix-like objects can be passed to CommonMatrix routines directly.
          if getattr(other, 'is_MatrixLike', False):
              return MatrixArithmetic._eval_add(self, other)
  
          raise TypeError('cannot add %s and %s' % (type(self), type(other)))
  
      @call_highest_priority('__rtruediv__')
      def __truediv__(self, other):
          return self * (self.one / other)
  
      @call_highest_priority('__rmatmul__')
      def __matmul__(self, other):
          other = _matrixify(other)
          if not getattr(other, 'is_Matrix', False) and not getattr(other, 'is_MatrixLike', False):
              return NotImplemented
  
          return self.__mul__(other)
  
      def __mod__(self, other):
          return self.applyfunc(lambda x: x % other)
  
      @call_highest_priority('__rmul__')
      def __mul__(self, other):
          """Return self*other where other is either a scalar or a matrix
          of compatible dimensions.
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
          >>> 2*A == A*2 == Matrix([[2, 4, 6], [8, 10, 12]])
          True
          >>> B = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
          >>> A*B
          Matrix([
          [30, 36, 42],
          [66, 81, 96]])
          >>> B*A
          Traceback (most recent call last):
          ...
          ShapeError: Matrices size mismatch.
          >>>
  
          See Also
          ========
  
          matrix_multiply_elementwise
          """
  
          return self.multiply(other)
  
      def multiply(self, other, dotprodsimp=None):
          """Same as __mul__() but with optional simplification.
  
          Parameters
          ==========
  
          dotprodsimp : bool, optional
              Specifies whether intermediate term algebraic simplification is used
              during matrix multiplications to control expression blowup and thus
              speed up calculation. Default is off.
          """
  
          isimpbool = _get_intermediate_simp_bool(False, dotprodsimp)
          other = _matrixify(other)
          # matrix-like objects can have shapes.  This is
          # our first sanity check. Double check other is not explicitly not a Matrix.
          if (hasattr(other, 'shape') and len(other.shape) == 2 and
              (getattr(other, 'is_Matrix', True) or
               getattr(other, 'is_MatrixLike', True))):
              if self.shape[1] != other.shape[0]:
                  raise ShapeError("Matrix size mismatch: %s * %s." % (
                      self.shape, other.shape))
  
          # honest sympy matrices defer to their class's routine
          if getattr(other, 'is_Matrix', False):
              m = self._eval_matrix_mul(other)
              if isimpbool:
                  return m._new(m.rows, m.cols, [_dotprodsimp(e) for e in m])
              return m
  
          # Matrix-like objects can be passed to CommonMatrix routines directly.
          if getattr(other, 'is_MatrixLike', False):
              return MatrixArithmetic._eval_matrix_mul(self, other)
  
          # if 'other' is not iterable then scalar multiplication.
          if not isinstance(other, Iterable):
              try:
                  return self._eval_scalar_mul(other)
              except TypeError:
                  pass
  
          return NotImplemented
  
      def multiply_elementwise(self, other):
          """Return the Hadamard product (elementwise product) of A and B
  
          Examples
          ========
  
          >>> from sympy.matrices import Matrix
          >>> A = Matrix([[0, 1, 2], [3, 4, 5]])
          >>> B = Matrix([[1, 10, 100], [100, 10, 1]])
          >>> A.multiply_elementwise(B)
          Matrix([
          [  0, 10, 200],
          [300, 40,   5]])
  
          See Also
          ========
  
          sympy.matrices.matrices.MatrixBase.cross
          sympy.matrices.matrices.MatrixBase.dot
          multiply
          """
          if self.shape != other.shape:
              raise ShapeError("Matrix shapes must agree {} != {}".format(self.shape, other.shape))
  
          return self._eval_matrix_mul_elementwise(other)
  
      def __neg__(self):
          return self._eval_scalar_mul(-1)
  
      @call_highest_priority('__rpow__')
      def __pow__(self, exp):
          """Return self**exp a scalar or symbol."""
  
          return self.pow(exp)
  
  
      def pow(self, exp, method=None):
          r"""Return self**exp a scalar or symbol.
  
          Parameters
          ==========
  
          method : multiply, mulsimp, jordan, cayley
              If multiply then it returns exponentiation using recursion.
              If jordan then Jordan form exponentiation will be used.
              If cayley then the exponentiation is done using Cayley-Hamilton
              theorem.
              If mulsimp then the exponentiation is done using recursion
              with dotprodsimp. This specifies whether intermediate term
              algebraic simplification is used during naive matrix power to
              control expression blowup and thus speed up calculation.
              If None, then it heuristically decides which method to use.
  
          """
  
          if method is not None and method not in ['multiply', 'mulsimp', 'jordan', 'cayley']:
              raise TypeError('No such method')
          if self.rows != self.cols:
              raise NonSquareMatrixError()
          a = self
          jordan_pow = getattr(a, '_matrix_pow_by_jordan_blocks', None)
          exp = sympify(exp)
  
          if exp.is_zero:
              return a._new(a.rows, a.cols, lambda i, j: int(i == j))
          if exp == 1:
              return a
  
          diagonal = getattr(a, 'is_diagonal', None)
          if diagonal is not None and diagonal():
              return a._new(a.rows, a.cols, lambda i, j: a[i,j]**exp if i == j else 0)
  
          if exp.is_Number and exp % 1 == 0:
              if a.rows == 1:
                  return a._new([[a[0]**exp]])
              if exp < 0:
                  exp = -exp
                  a = a.inv()
          # When certain conditions are met,
          # Jordan block algorithm is faster than
          # computation by recursion.
          if method == 'jordan':
              try:
                  return jordan_pow(exp)
              except MatrixError:
                  if method == 'jordan':
                      raise
  
          elif method == 'cayley':
              if not exp.is_Number or exp % 1 != 0:
                  raise ValueError("cayley method is only valid for integer powers")
              return a._eval_pow_by_cayley(exp)
  
          elif method == "mulsimp":
              if not exp.is_Number or exp % 1 != 0:
                  raise ValueError("mulsimp method is only valid for integer powers")
              return a._eval_pow_by_recursion_dotprodsimp(exp)
  
          elif method == "multiply":
              if not exp.is_Number or exp % 1 != 0:
                  raise ValueError("multiply method is only valid for integer powers")
              return a._eval_pow_by_recursion(exp)
  
          elif method is None and exp.is_Number and exp % 1 == 0:
              # Decide heuristically which method to apply
              if a.rows == 2 and exp > 100000:
                  return jordan_pow(exp)
              elif _get_intermediate_simp_bool(True, None):
                  return a._eval_pow_by_recursion_dotprodsimp(exp)
              elif exp > 10000:
                  return a._eval_pow_by_cayley(exp)
              else:
                  return a._eval_pow_by_recursion(exp)
  
          if jordan_pow:
              try:
                  return jordan_pow(exp)
              except NonInvertibleMatrixError:
                  # Raised by jordan_pow on zero determinant matrix unless exp is
                  # definitely known to be a non-negative integer.
                  # Here we raise if n is definitely not a non-negative integer
                  # but otherwise we can leave this as an unevaluated MatPow.
                  if exp.is_integer is False or exp.is_nonnegative is False:
                      raise
  
          from sympy.matrices.expressions import MatPow
          return MatPow(a, exp)
  
      @call_highest_priority('__add__')
      def __radd__(self, other):
          return self + other
  
      @call_highest_priority('__matmul__')
      def __rmatmul__(self, other):
          other = _matrixify(other)
          if not getattr(other, 'is_Matrix', False) and not getattr(other, 'is_MatrixLike', False):
              return NotImplemented
  
          return self.__rmul__(other)
  
      @call_highest_priority('__mul__')
      def __rmul__(self, other):
          return self.rmultiply(other)
  
      def rmultiply(self, other, dotprodsimp=None):
          """Same as __rmul__() but with optional simplification.
  
          Parameters
          ==========
  
          dotprodsimp : bool, optional
              Specifies whether intermediate term algebraic simplification is used
              during matrix multiplications to control expression blowup and thus
              speed up calculation. Default is off.
          """
          isimpbool = _get_intermediate_simp_bool(False, dotprodsimp)
          other = _matrixify(other)
          # matrix-like objects can have shapes.  This is
          # our first sanity check. Double check other is not explicitly not a Matrix.
          if (hasattr(other, 'shape') and len(other.shape) == 2 and
              (getattr(other, 'is_Matrix', True) or
               getattr(other, 'is_MatrixLike', True))):
              if self.shape[0] != other.shape[1]:
                  raise ShapeError("Matrix size mismatch.")
  
          # honest sympy matrices defer to their class's routine
          if getattr(other, 'is_Matrix', False):
              m = self._eval_matrix_rmul(other)
              if isimpbool:
                  return m._new(m.rows, m.cols, [_dotprodsimp(e) for e in m])
              return m
          # Matrix-like objects can be passed to CommonMatrix routines directly.
          if getattr(other, 'is_MatrixLike', False):
              return MatrixArithmetic._eval_matrix_rmul(self, other)
  
          # if 'other' is not iterable then scalar multiplication.
          if not isinstance(other, Iterable):
              try:
                  return self._eval_scalar_rmul(other)
              except TypeError:
                  pass
  
          return NotImplemented
  
      @call_highest_priority('__sub__')
      def __rsub__(self, a):
          return (-self) + a
  
      @call_highest_priority('__rsub__')
      def __sub__(self, a):
          return self + (-a)
  
  class MatrixCommon(MatrixArithmetic, MatrixOperations, MatrixProperties,
                    MatrixSpecial, MatrixShaping):
      """All common matrix operations including basic arithmetic, shaping,
      and special matrices like `zeros`, and `eye`."""
      _diff_wrt = True  # type: bool
  
  
  class _MinimalMatrix:
      """Class providing the minimum functionality
      for a matrix-like object and implementing every method
      required for a `MatrixRequired`.  This class does not have everything
      needed to become a full-fledged SymPy object, but it will satisfy the
      requirements of anything inheriting from `MatrixRequired`.  If you wish
      to make a specialized matrix type, make sure to implement these
      methods and properties with the exception of `__init__` and `__repr__`
      which are included for convenience."""
  
      is_MatrixLike = True
      _sympify = staticmethod(sympify)
      _class_priority = 3
      zero = S.Zero
      one = S.One
  
      is_Matrix = True
      is_MatrixExpr = False
  
      @classmethod
      def _new(cls, *args, **kwargs):
          return cls(*args, **kwargs)
  
      def __init__(self, rows, cols=None, mat=None, copy=False):
          if isfunction(mat):
              # if we passed in a function, use that to populate the indices
              mat = list(mat(i, j) for i in range(rows) for j in range(cols))
          if cols is None and mat is None:
              mat = rows
          rows, cols = getattr(mat, 'shape', (rows, cols))
          try:
              # if we passed in a list of lists, flatten it and set the size
              if cols is None and mat is None:
                  mat = rows
              cols = len(mat[0])
              rows = len(mat)
              mat = [x for l in mat for x in l]
          except (IndexError, TypeError):
              pass
          self.mat = tuple(self._sympify(x) for x in mat)
          self.rows, self.cols = rows, cols
          if self.rows is None or self.cols is None:
              raise NotImplementedError("Cannot initialize matrix with given parameters")
  
      def __getitem__(self, key):
          def _normalize_slices(row_slice, col_slice):
              """Ensure that row_slice and col_slice don't have
              `None` in their arguments.  Any integers are converted
              to slices of length 1"""
              if not isinstance(row_slice, slice):
                  row_slice = slice(row_slice, row_slice + 1, None)
              row_slice = slice(*row_slice.indices(self.rows))
  
              if not isinstance(col_slice, slice):
                  col_slice = slice(col_slice, col_slice + 1, None)
              col_slice = slice(*col_slice.indices(self.cols))
  
              return (row_slice, col_slice)
  
          def _coord_to_index(i, j):
              """Return the index in _mat corresponding
              to the (i,j) position in the matrix. """
              return i * self.cols + j
  
          if isinstance(key, tuple):
              i, j = key
              if isinstance(i, slice) or isinstance(j, slice):
                  # if the coordinates are not slices, make them so
                  # and expand the slices so they don't contain `None`
                  i, j = _normalize_slices(i, j)
  
                  rowsList, colsList = list(range(self.rows))[i], \
                                       list(range(self.cols))[j]
                  indices = (i * self.cols + j for i in rowsList for j in
                             colsList)
                  return self._new(len(rowsList), len(colsList),
                                   list(self.mat[i] for i in indices))
  
              # if the key is a tuple of ints, change
              # it to an array index
              key = _coord_to_index(i, j)
          return self.mat[key]
  
      def __eq__(self, other):
          try:
              classof(self, other)
          except TypeError:
              return False
          return (
              self.shape == other.shape and list(self) == list(other))
  
      def __len__(self):
          return self.rows*self.cols
  
      def __repr__(self):
          return "_MinimalMatrix({}, {}, {})".format(self.rows, self.cols,
                                                     self.mat)
  
      @property
      def shape(self):
          return (self.rows, self.cols)
  
  
  class _CastableMatrix: # this is needed here ONLY FOR TESTS.
      def as_mutable(self):
          return self
  
      def as_immutable(self):
          return self
  
  
  class _MatrixWrapper:
      """Wrapper class providing the minimum functionality for a matrix-like
      object: .rows, .cols, .shape, indexability, and iterability. CommonMatrix
      math operations should work on matrix-like objects. This one is intended for
      matrix-like objects which use the same indexing format as SymPy with respect
      to returning matrix elements instead of rows for non-tuple indexes.
      """
  
      is_Matrix     = False # needs to be here because of __getattr__
      is_MatrixLike = True
  
      def __init__(self, mat, shape):
          self.mat = mat
          self.shape = shape
          self.rows, self.cols = shape
  
      def __getitem__(self, key):
          if isinstance(key, tuple):
              return sympify(self.mat.__getitem__(key))
  
          return sympify(self.mat.__getitem__((key // self.rows, key % self.cols)))
  
      def __iter__(self): # supports numpy.matrix and numpy.array
          mat = self.mat
          cols = self.cols
  
          return iter(sympify(mat[r, c]) for r in range(self.rows) for c in range(cols))
  
  
  class MatrixKind(Kind):
      """
      Kind for all matrices in SymPy.
  
      Basic class for this kind is ``MatrixBase`` and ``MatrixExpr``,
      but any expression representing the matrix can have this.
  
      Parameters
      ==========
  
      element_kind : Kind
          Kind of the element. Default is :obj:NumberKind `<sympy.core.kind.NumberKind>`,
          which means that the matrix contains only numbers.
  
      Examples
      ========
  
      Any instance of matrix class has ``MatrixKind``.
  
      >>> from sympy import MatrixSymbol
      >>> A = MatrixSymbol('A', 2,2)
      >>> A.kind
      MatrixKind(NumberKind)
  
      Although expression representing a matrix may be not instance of
      matrix class, it will have ``MatrixKind`` as well.
  
      >>> from sympy import Integral
      >>> from sympy.matrices.expressions import MatrixExpr
      >>> from sympy.abc import x
      >>> intM = Integral(A, x)
      >>> isinstance(intM, MatrixExpr)
      False
      >>> intM.kind
      MatrixKind(NumberKind)
  
      Use ``isinstance()`` to check for ``MatrixKind` without specifying
      the element kind. Use ``is`` with specifying the element kind.
  
      >>> from sympy import Matrix
      >>> from sympy.matrices import MatrixKind
      >>> from sympy.core.kind import NumberKind
      >>> boolM = Matrix([True, False])
      >>> isinstance(boolM.kind, MatrixKind)
      True
      >>> boolM.kind is MatrixKind(NumberKind)
      False
  
      See Also
      ========
  
      shape : Function to return the shape of objects with ``MatrixKind``.
  
      """
      def __new__(cls, element_kind=NumberKind):
          obj = super().__new__(cls, element_kind)
          obj.element_kind = element_kind
          return obj
  
      def __repr__(self):
          return "MatrixKind(%s)" % self.element_kind
  
  
  def _matrixify(mat):
      """If `mat` is a Matrix or is matrix-like,
      return a Matrix or MatrixWrapper object.  Otherwise
      `mat` is passed through without modification."""
  
      if getattr(mat, 'is_Matrix', False) or getattr(mat, 'is_MatrixLike', False):
          return mat
  
      if not(getattr(mat, 'is_Matrix', True) or getattr(mat, 'is_MatrixLike', True)):
          return mat
  
      shape = None
  
      if hasattr(mat, 'shape'): # numpy, scipy.sparse
          if len(mat.shape) == 2:
              shape = mat.shape
      elif hasattr(mat, 'rows') and hasattr(mat, 'cols'): # mpmath
          shape = (mat.rows, mat.cols)
  
      if shape:
          return _MatrixWrapper(mat, shape)
  
      return mat
  
  
  def a2idx(j, n=None):
      """Return integer after making positive and validating against n."""
      if type(j) is not int:
          jindex = getattr(j, '__index__', None)
          if jindex is not None:
              j = jindex()
          else:
              raise IndexError("Invalid index a[%r]" % (j,))
      if n is not None:
          if j < 0:
              j += n
          if not (j >= 0 and j < n):
              raise IndexError("Index out of range: a[%s]" % (j,))
      return int(j)
  
  
  def classof(A, B):
      """
      Get the type of the result when combining matrices of different types.
  
      Currently the strategy is that immutability is contagious.
  
      Examples
      ========
  
      >>> from sympy import Matrix, ImmutableMatrix
      >>> from sympy.matrices.common import classof
      >>> M = Matrix([[1, 2], [3, 4]]) # a Mutable Matrix
      >>> IM = ImmutableMatrix([[1, 2], [3, 4]])
      >>> classof(M, IM)
      <class 'sympy.matrices.immutable.ImmutableDenseMatrix'>
      """
      priority_A = getattr(A, '_class_priority', None)
      priority_B = getattr(B, '_class_priority', None)
      if None not in (priority_A, priority_B):
          if A._class_priority > B._class_priority:
              return A.__class__
          else:
              return B.__class__
  
      try:
          import numpy
      except ImportError:
          pass
      else:
          if isinstance(A, numpy.ndarray):
              return B.__class__
          if isinstance(B, numpy.ndarray):
              return A.__class__
  
      raise TypeError("Incompatible classes %s, %s" % (A.__class__, B.__class__))
  