function [X,Y] = sda_affine_mnare(A,B,C,D)
% [X,Y]=SDA_AFFINE_MNARE(A,B,C,D) solves the M-NARE C + XA + DX - XBX = 0
% by means of SDA based on the affine transform
%    A, B, C, D: matrix coefficients
%    X : solution of the M-NARE
%    Y : solution of the dual M-NARE
a = 1/max(diag(A));
n = size(A,1); m = size(D,1);
% Initialization by means of the affine transform
Da = inv(eye(m) + a*D);
F = -Da;
G = a*B*Da;
P = -a*Da*C;
E = -eye(n) + a*A + a*G*C;
% SDA step
[X,Y] = sda(E,F,G,P);