'Solves Ax = B using LDL decomposition,
for a general square and non-singular matrix.
For a non-square matrix with rows > cols,
the least squares solution is returned.
Examples
>>> from sympy.matrices import Matrix, eye
>>> A = eye(2)*2
>>> B = Matrix([[1, 2], [3, 4]])
>>> A.LDLsolve(B) == B/2
True
See Also
LDLdecomposition
lower_triangular_solve
upper_triangular_solve
gauss_jordan_solve
cholesky_solve
diagonal_solve
LUsolve
QRsolve
pinv_solve'

def LDLsolve(self, rhs):

    if self.is_symmetric():
        (L, D) = self.LDLdecomposition()
    elif (self.rows >= self.cols):
        (L, D) = (self.T * self).LDLdecomposition()
        rhs = (self.T * rhs)
    else:
        raise NotImplementedError('Under-determined    System.    Try    M.gauss_jordan_solve(rhs)')
    Y = L._lower_triangular_solve(rhs)
    Z = D._diagonal_solve(Y)
    return L.T._upper_triangular_solve(Z)
