def d_operator(self, P):
        r"""
        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

        """
        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
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
r"""
        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