function [X,Y] = cr_cayley_mnare(A,B,C,D,g)
% [X,Y]=CR_CAYLEY_MNARE(A,B,C,D,g) solves the M-NARE C + XA + DX - XBX = 0
% by means of CR based on the Cayley transform
%    A, B, C, D: matrix coefficients
%    g: scalar gamma used by the Cayley transform
%    X : solution of the M-NARE
%    Y : solution of the dual M-NARE
n = size(A,1); m = size(D,1);
S1 = (A+g*eye(n))\B;
S2 = (g*eye(m)-D)\C;
V = -D + g*eye(m) - C*S1;
W = A + g*eye(n) - B*S2;
VI = inv(V);
WI = inv(W);
R1 = eye(n) - 2*g*WI;
R2 = -2*g*S2*WI;
R3 = -2*g*S1*VI;
R4 = zeros(m,n);
R5 = eye(m) - 2*g*VI;
R6 = -eye(m);
% CR iteration
[X,Y] = cr_struct(R1,R2,R3,R4,R5,R6);