function [X,Y] = cr_affine_mnare(A,B,C,D)
% [X,Y]=CR_AFFINE_MNARE(A,B,C,D) solves the M-NARE C + XA + DX - XBX = 0
% by means of CR based on the affine transform
%    A, B, C, D: matrix coefficients
%    X : solution of the M-NARE
%    Y : solution of the dual M-NARE
alpha = 1/max(diag(A));
n = size(A,1); m = size(D,1);
R1 = -eye(n) + alpha*A;
R2 = -alpha*C;
R3 = -alpha*B;
R4 = zeros(m,n);
R5 = -eye(m) - alpha*D;
R6 = -eye(m);
% CR iteration
[X,Y] = cr_struct(R1,R2,R3,R4,R5,R6);