function [X,res,iter] = mriccati (A,B,C,D,alg,singular,shift)

%MRICCATI  Solves an M-Riccati equation by various algorithm
%
%   [X res iter] = MRICCATI(A,B,C,D,alg) uses various algorithms
%   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 ]
%   alg: the algorithm of choice
%
%   Implemented algorithms
%       'newt' Newton's method (requires sylvsol.m)
%       'sdas' Structure-preserving doubling algorithm
%       'bilm' CR applied to the UQME
%
% Output: X     : the minimal nonnegative solution
%         res   : a relative residual
%         iter  : the number of iterative steps

if nargin<5
    alg='newt';
end

if nargin<6
    singular='y';
end

[m n]=size(B);
[s d u v]=mriccati_detect(A,B,C,D);
if (size(u,1)>0)
    u(n+1:n+m)=-u(n+1:n+m);
end

if (s=='n') % If M is nonsingular don't use singularity
    singular='n';
end

if (s=='t' & singular=='y') % Reverse the coefficients in the transient case
    T=A;A=D';B=B';C=C';D=T';
    v(1:m)=u(n+1:n+m);v(m+1:m+n)=u(1:n);
end    
    
switch alg
    
case 'newt' % Newton's method
    [X iter]=mriccati_alg_newton(A,B,C,D,v,singular,d);

case 'sdas' % Structure-preserving Doubling Algorithm (SDA)
    [X iter]=mriccati_alg_sda(A,B,C,D,v,singular);

case 'bilm' % Transformation to UQME and Cyclic Reduction
    [X iter]=mriccati_alg_bilm(A,B,C,D,v,singular);
    
end

if (s=='t' & singular=='y') % Reverse the coefficients in the transient case
    T=A;A=D';B=B';C=C';D=T';
    X=X';
end    

if (nargout>1) 
    XC=X*C;AX=A*X;XD=X*D;XCX=XC*X;
    res=norm(XCX-AX-XD+B,1)/(norm(XCX,1)+norm(AX,1)+norm(XD,1)+norm(B,1));
end