function [p1,p2,lamp,perr,wellposed,p] = daresolv(A,B,Q,R,Type)
%
% [P1,P2,LAMP,PERR,WELLPOSED,P] = DARESOLV(A,B,Q,R,TYPE) computes the 
% solution to the discrete algebraic Riccati equation
%                                   -1
%            A'PA - P - A'PB(R+B'PB)  B'PA + Q = 0
%
%   based on the "Type" one selects:
% 
%                Type = 'eigen' ---- Eigenstructure approach
%
%                Type = 'Schur' ---- Schur vector approach
%
%  Outputs:   P = P2/P1 riccati solution
%             [P1;P2] stable eigenspace of hamiltonian matrix
%             LAMP closed loop eigenvalues
%             PERR residual error matrix 
%             WELLPOSED = 'TRUE ' if hamiltonian has no imaginary-axis
%                  eigenvalues, otherwise 'FALSE'

% R. Y. Chiang & M. G. Safonov 2/88
% Copyright (c) 1988 by the MathWorks, Inc.
% All Rights Reserved.

wellposed='TRUE ';  %initialize 
if nargin == 3
   Type = 'eigen';
end
%
[n,n] = size(A);
%
% ------ If A is stable and Q is zero, p2 = zero(n), p1 = eye(n)
%
lamA = eig(A);
ind = find(lamA>0);
if (sum(ind) == 0) & (norm(q,'fro') == 0)
   p2 = zeros(n); p1 = eye(n);
   p = p2; perr = p2;
   return
end
%   
RR = B/R*B';
ham = [A+RR/A'*Q -RR/A';-A'\Q inv(A')];
%
% ------ Eigenstructure approach:
%
if Type == 'eigen'
   [v,d] = reig(ham,2);
   if rank(v) < 2*n
      Type = 'Schur';
   end
   p1 = v(1:n,1:n);
   p2 = v((n+1):(2*n),1:n);
end
%
% ------ Schur vector approach:
%
if Type == 'Schur'
   [v,t,m,swap] = cschur(ham,6);
   vv = [real(v(:,1:m)) imag(v(:,1:m))];
   [vh,rh] = qr(vv);
   p1 = vh(1:n,1:n);
   p2 = vh(n+1:2*n,1:n);
end
%
P = real(p2/p1);
%
perr = A'*P*A-P-A'*P*B/(R+B'*P*B)*B'*P*A+Q;
%
% ------ Check jw-axis roots:
%
lamp = eig(A-B/(R+B'*P*B)*B'*P'*A);
%
if max(abs(lamp)) > 1-1.e-8
     wellposed='FALSE';
end
%
% ------- End of DARESOLV.M ---- % RYC/MGS %