function [b,rts,ia,nexact,nnumeric,lgroots,mcode] = AiMEigenvalues(h,neq,nlag,nlead,condn,uprbnd)
% AiMEigenvalues: Solve a linear perfect foresight model using the matlab eig
%                 function to find the invariant subspace associated with the big
%                 roots.  This procedure will fail if the companion matrix is
%                 defective and does not have a linearly independent set of
%                 eigenvectors associated with the big roots.
%
% USAGE:
% 
%       [b,rts,ia,nexact,nnumeric,lgroots,mcode] = AiMEigenvalues(h,neq,nlag,nlead,condn,uprbnd)
%
% REQUIRED INPUT:  h (matrix) with structural coefficients - dimension (neq,neq*(nlag+1+nlead)).
%
%                  neq (integer) with the number of equations.
%
%                  nlag (integer) with the number of lags.
%
%                  nlead (integer) with the number of leads.
%
%                  condn (scalar) tolerance used as a condition number test by AiMNumericShift and AiMReducedForm
%
%                  uprbnd (scalar) upper bound for the modulus of roots allowed in the reduced form.
%
% REQUIRED OUTPUT: b (matrix) with reduced form coefficients - dimension (neq,neq*nlag).
%
%                  rts (vector) with roots returned by eig.
%
%                  ia (integer) Dimension of companion matrix (number of non-trivial elements in rts).
%
%                  nexact (integer) with number of exact shiftrights.
%
%                  nnumeric (integer) with number of numeric shiftrights.
%
%                  lgroots (integer) with number of roots greater in modulus than uprbnd.
%
%                  mcode (integer) with return code.
%
% NOTE: This function is based on aim_eig.m
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2009 European Central Bank.
%
%                       First version: August 25, 2006.
%                        This version: March 3, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 19-09-2006: Initialized the ia integer to 0.
%
% * 30-10-2007: Updated the documentation.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: Added the mcode value 7 when a has infinite or NaN entries. This avoids
%               running the eig function on such a matrix.
%
% * 05-12-2008: Updated the documentation.
%
% * 03-03-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% Initialization.
%
nexact   = 0;
nnumeric = 0;
lgroots  = 0;
iq       = 0;
mcode    = 0;
ia       = 0;
%
qrows = neq*nlead;
qcols = neq*(nlag+nlead);
bcols = neq*nlag;
%
q        = zeros(qrows,qcols);
rts      = zeros(qcols,1);
%
b = zeros(neq,bcols);
%
% Compute the auxiliary initial conditions and store them in q.
%
[h,q,iq,nexact] = AiMExactShift(h,q,iq,qrows,qcols,neq);
%
if (iq>qrows);
   mcode = 61;
   return;
end;
%
[h,q,iq,nnumeric] = AiMNumericShift(h,q,iq,qrows,qcols,neq,condn);
%
if (iq>qrows);
   mcode = 62;
   return;
end;
%
%  Build the companion matrix.  Compute the stability conditions, and
%  combine them with the auxiliary initial conditions in q.  
%
[a,ia,js] = AiMBuildCompanion(h,qcols,neq);
%
if (ia~=0);
   if (max(max(isinf(a)))==1)|(max(max(isnan(a)))==1);
      mcode = 7;
      return;
   end;
   [w,rts,lgroots] = AiMEigenSystem(a,uprbnd);
   q = AiMCopyEigenvectors(q,w,js,iq,qrows);
end;
%
test = nexact+nnumeric+lgroots;
%
if (test>qrows);
   mcode = 3;
elseif (test<qrows);
   mcode = 4;
end;
%
% If the right-hand block of q is invertible, compute the reduced form.
%
[nonsing,b] = AiMReducedForm(q,qrows,qcols,bcols,neq,b,condn);
%
if (nonsing&mcode==0);
   mcode = 1;
elseif (~nonsing&mcode==0);
   mcode = 5;
elseif (~nonsing&mcode==3);
   mcode = 35;
elseif (~nonsing&mcode==4);
   mcode = 45;
end

%
% end of AiMEigenvalues.m
%
