function  [f,v] = dlqreigg(a,b,q,r)
%
%DLQREIGG  Discrete linear quadratic regulator solution
%  [f,v] = dlqreigg(a,b,q,r) solves the Riccati equation
%  associated with the real discrete time optimal linear regulator
%  problem of the form
%
%        min     sum {x'qx + u'ru}
%        st      x(t+1) = ax(t) + bu(t)
%
%  Returned are the matrices f and v such that the optimal decision rule
%  is -fx and the value function is x'vx
%
%  ALGORITHM: Generalized Eigenvector
%  WRITTEN BY:  Evan W. Anderson   December 1, 1994



%Get the ordered generalized eigenvector decomposition  of (L,N) 

tol = 1e-15; 
n=size(a,1);
[V,D] = eig([a zeros(n); -q eye(n)],[eye(n) (b/r)*b'; zeros(n) a']);

%Calculate the solution to the Riccati equation and the optimal
% decision rule. 

[E,I] = sort(abs(diag(D)));
v = V(n+1:2*n,I(1:n))/V(1:n,I(1:n));
f = (b'*v*b + r)\b'*v*a;