code: def _LDLdecomposition(M, hermitian=True):
      """Returns the LDL Decomposition (L, D) of matrix A,
      such that L * D * L.H == A if hermitian flag is True, or
      L * D * L.T == A if hermitian is False.
      This method eliminates the use of square root.
      Further this ensures that all the diagonal entries of L are 1.
      A must be a Hermitian positive-definite matrix if hermitian is True,
      or a symmetric matrix otherwise.
  
      Examples
      ========
  
      >>> from sympy import Matrix, eye
      >>> A = Matrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))
      >>> L, D = A.LDLdecomposition()
      >>> L
      Matrix([
      [   1,   0, 0],
      [ 3/5,   1, 0],
      [-1/5, 1/3, 1]])
      >>> D
      Matrix([
      [25, 0, 0],
      [ 0, 9, 0],
      [ 0, 0, 9]])
      >>> L * D * L.T * A.inv() == eye(A.rows)
      True
  
      The matrix can have complex entries:
  
      >>> from sympy import I
      >>> A = Matrix(((9, 3*I), (-3*I, 5)))
      >>> L, D = A.LDLdecomposition()
      >>> L
      Matrix([
      [   1, 0],
      [-I/3, 1]])
      >>> D
      Matrix([
      [9, 0],
      [0, 4]])
      >>> L*D*L.H == A
      True
  
      See Also
      ========
  
      sympy.matrices.dense.DenseMatrix.cholesky
      sympy.matrices.matrixbase.MatrixBase.LUdecomposition
      QRdecomposition
      """
  
      from .dense import MutableDenseMatrix
  
      if not M.is_square:
          raise NonSquareMatrixError("Matrix must be square.")
      if hermitian and not M.is_hermitian:
          raise ValueError("Matrix must be Hermitian.")
      if not hermitian and not M.is_symmetric():
          raise ValueError("Matrix must be symmetric.")
  
      D = MutableDenseMatrix.zeros(M.rows, M.rows)
      L = MutableDenseMatrix.eye(M.rows)
  
      if hermitian:
          for i in range(M.rows):
              for j in range(i):
                  L[i, j] = (1 / D[j, j]) * (
                      M[i, j]
                      - sum(
                          L[i, k] * L[j, k].conjugate() * D[k, k]
                          for k in range(j)
                      )
                  )
  
              D[i, i] = M[i, i] - sum(
                  L[i, k] * L[i, k].conjugate() * D[k, k]
                  for k in range(i)
              )
  
              if D[i, i].is_positive is False:
                  raise NonPositiveDefiniteMatrixError(
                      "Matrix must be positive-definite"
                  )
  
      else:
          for i in range(M.rows):
              for j in range(i):
                  L[i, j] = (1 / D[j, j]) * (
                      M[i, j]
                      - sum(
                          L[i, k] * L[j, k] * D[k, k] for k in range(j)
                      )
                  )
  
              D[i, i] = M[i, i] - sum(
                  L[i, k] ** 2 * D[k, k] for k in range(i)
              )
  
      return M._new(L), M._new(D)
  
imports: import copy
  from sympy.core import S
  from sympy.core.function import expand_mul
  from sympy.functions.elementary.miscellaneous import Min, sqrt
  from sympy.functions.elementary.complexes import sign
  from .exceptions import NonSquareMatrixError, NonPositiveDefiniteMatrixError
  from .utilities import _get_intermediate_simp, _iszero
  from .determinant import _find_reasonable_pivot_naive
label_desc: Use image node_id 6 for calling a global function with example usage: _LDLdecomposition(M, hermitian) and returns: M,  M