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

%MRICCATI_ALG_NEWTON  Solves an M-Riccati equation by Newton's method
%
%   [X iter res] = MRICCATI_ALG_NEWTON(A,B,C,D) uses Newton's 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 singular='y' then the initial matrix 
%   is chosen in order to reduce the number of steps, or else it is the 
%   zero matrix
%
%   Output
%   X     : the minimal nonnegative solution
%   iter  : the number of iterative steps
%
%   Stopping criterion: if the norm of the increment is less than a
%   tolerance or increases stop the iteration

if nargin<6
    singular='n';
end

if nargin<7
    d=1;
end

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

[m n]=size(B);

X=zeros(m,n);
if (singular=='y') % If singularity must be used
    v1=v(1:n);v2=v(n+1:n+m);
    s=sum(v1);
    for k=1:m
        X(k,1:n)=v2(k)/s*ones(1,n);
    end
end

gam=max(max(diag(A)),max(diag(D)));
AA=A;DD=D;CC=C;BB=B;
if (abs(d)<1d-2 &  singular=='y') % If singularity must be used
    eta=gam;
    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

H0=ones(m,n);
iter=0;
res=1;
h0=inf;  

while (iter<maxiter)
    iter=iter+1;
    XC=X*CC;CX=CC*X;
    H0=lyap(AA-XC,DD-CX,-BB-XC*X+AA*X+X*DD);
    h1=norm(H0);
    
    if (h1>h0 | h1<tol)
        break;
    end
    X=X+H0;
    h0=h1;
end