'The D operator, mapping P into
.. math::
D(P) := P + PC(\theta I - C\'PC)^{-1} C\'P.
Parameters
P : array_like(float, ndim=2)
A matrix that should be n x n
Returns
dP : array_like(float, ndim=2)
The matrix P after applying the D operator'

def d_operator(self, P):

    (C, theta) = (self.C, self.theta)
    I = np.identity(self.j)
    S1 = dot(P, C)
    S2 = dot(C.T, S1)
    dP = (P + dot(S1, solve(((theta * I) - S2), S1.T)))
    return dP
