'Return (P, D), where D is diagonal and
D = P^-1 * M * P
where M is current matrix.
Parameters
reals_only : bool. Whether to throw an error if complex numbers are need
to diagonalize. (Default: False)
sort : bool. Sort the eigenvalues along the diagonal. (Default: False)
normalize : bool. If True, normalize the columns of P. (Default: False)
Examples
>>> from sympy import Matrix
>>> m = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
>>> m
Matrix([
[1,  2, 0],
[0,  3, 0],
[2, -4, 2]])
>>> (P, D) = m.diagonalize()
>>> D
Matrix([
[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
>>> P
Matrix([
[-1, 0, -1],
[ 0, 0, -1],
[ 2, 1,  2]])
>>> P.inv() * m * P
Matrix([
[1, 0, 0],
[0, 2, 0],
[0, 0, 3]])
See Also
is_diagonal
is_diagonalizable'

def diagonalize(self, reals_only=False, sort=False, normalize=False):

    if (not self.is_square):
        raise NonSquareMatrixError()
    if (not self.is_diagonalizable(reals_only=reals_only, clear_cache=False)):
        raise MatrixError('Matrix    is    not    diagonalizable')
    eigenvecs = self._cache_eigenvects
    if (eigenvecs is None):
        eigenvecs = self.eigenvects(simplify=True)
    if sort:
        eigenvecs = sorted(eigenvecs, key=default_sort_key)
    (p_cols, diag) = ([], [])
    for (val, mult, basis) in eigenvecs:
        diag += ([val] * mult)
        p_cols += basis
    if normalize:
        p_cols = [(v / v.norm()) for v in p_cols]
    return (self.hstack(*p_cols), self.diag(*diag))
