function [X,iter] = mriccati_alg_bilm (A,B,C,D,v,singular)

%MRICCATI_ALG_BILM  Solves an M-Riccati equation by transforming it into
%                   a unilateral quadratic matrix equation an the applying
%                   the Cyclic Reduction algorithm
%
%   [X iter res] = MRICCATI_ALG_BILM(A,B,C,D) use the CR method
%   to compute the minimal nonnegative solution of the nonsymmetric 
%   algebraic Riccati equation
%               XCX-AX-XD+B=0
%   whose coefficients form an M-matrix
%               [ D  -C ]
%               [ -B  A ]
%   If the optional parameter is different from 's' the shift technique is avoided
%
%   Output
%   X     : the minimal nonnegative solution
%   iter  : the number of iterative steps
%
%   Stopping criterion: if the matrices of the CR converging to zero have
%   norm smaller than a fixed tolerance tol, then convergence is declared

if nargin<6
    singular='n';
end

maxiter=30;
tol=1.d-15;

[m n]=size(B);
In=eye(n);Im=eye(m);Onm=zeros(n,m);Om=zeros(m);On=zeros(n);

gamma=max(max(diag(A)),max(diag(D)));

AA=A;DD=D;CC=C;BB=B;
if (singular=='y') % If singularity must be used
    eta=gamma;
    a1=v(1:n); % a is the vector v
    a2=v(n+1:n+m);
    b1=ones(n,1)/sum(v); % b is a vector p such that p^Tv=1
    b2=ones(m,1)/sum(v);
    AA=A-eta*a2*b2';
    BB=B+eta*a2*b1';
    CC=C-eta*a1*b2';
    DD=D+eta*a1*b1';
end

t=1/gamma;
A0=[Im-t*A t*B;Onm On];
B0=[-Im-t*A  2*t*B+2/n*ones(m,n)-1/n*t*B*ones(n,n);t*C -In-t*D];BB0=B0;
C0=[Om t*B-1/n*t*B*ones(n,n);Onm In-t*D-1/n*ones(n)+1/n*t*D*ones(n)];
K=[Om t*B;Onm In-t*D];

iter=0;
while (norm(C0,1)*norm(A0,1)>tol^2)
    iter=iter+1;
    Binv=inv(B0);
    
    % block decomposition 
    B11=Binv(1:m,1:m);B12=Binv(1:m,m+1:m+n);
    B21=Binv(m+1:m+n,1:m);B22=Binv(m+1:m+n,m+1:m+n);
    C12=C0(1:m,m+1:m+n);C22=C0(m+1:m+n,m+1:m+n);
    A11=A0(1:m,1:m);A12=A0(1:m,m+1:m+n);
    
    % updates computation
    A1B1=A11*B11;A1B2=A11*B12;A2B3=A12*B21;A2B4=A12*B22;
    C1B3=C12*B21;C2B3=C22*B21;C1B4=C12*B22;C2B4=C22*B22;
    
    ABA=[A1B1*A11+A2B3*A11 A1B1*A12+A2B3*A12;Onm On];
    CBC=[Om C1B3*C12+C1B4*C22;Onm C2B3*C12+C2B4*C22];
    CBA=[C1B3*A11 C1B3*A12;C2B3*A11 C2B3*A12];
    ABC=[Om A1B1*C12+A2B3*C12+A1B2*C22+A2B4*C22;Onm On];
    
    % iteration
    A1=-ABA;B1=B0-ABC-CBA;BB1=BB0-ABC;C1=-CBC;
    A0=A1;B0=B1;BB0=BB1;C0=C1;
    
    % to see convergence norm(C0,1);
end

T=-inv(BB0)*K;

X=T(1:m,m+1:m+n)*inv(In+T(m+1:m+n,m+1:m+n));