Added use of G Anderson AIM - see Readme.txt in the new AIM subdirectory.

DR1 tries to invoke AIM if options_.useAIM == 1 is set.
for start, options_.useAIM = 0 is set by default in global_initialization.m so that system uses mjdgges by default. 
If AIM is to be used, options_.useAIM = 1 needs to be set either in the model <>.mod file, before invoking, estimate and/or stoch_simul, or by issuing appropriate command for estimate and/or stoch_simul. 

However, if other exceptions occur and aimcode is higher than 5, the system resets options_.useAIM = 0 tries to use mjdgges instead.

GP July '08

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1967 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/dr1.m

@@ -209,6 +209,24 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
     b = jacobia_(:,k0);
     
     if M_.maximum_endo_lead == 0;  % backward models
+        % If required, try Gary Anderson and G Moore AIM solver if not
+        % check only and if 1st order (added by GP July'08)
+        if (options_.useAIM == 1) && (task == 0) && (options_.order == 1) 
+            try
+                [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr);
+                if aimcode ~=1
+                    info(1) = aimcode;
+                    info(2) = 1.0e+8;
+                    return
+                end
+            catch
+                %warning('Problem with using AIM solver - Using Dynare solver instead');
+                disp('Problem with using AIM solver - Using Dynare solver instead');
+                options_.useAIM = 0; % and try mjdgges instead
+            end
+        end % end if useAIM and...
+        %else use original Dynare solver
+        if ~((options_.useAIM == 1) && (task == 0) && (options_.order == 1)) 
             a = jacobia_(:,nonzeros(k1'));
             dr.ghx = zeros(size(a));
             m = 0;
@@ -220,6 +238,7 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
             if M_.exo_nbr
                 dr.ghu = -b\jacobia_(:,nz+1:end);
             end
+        end % if not use AIM or not...
         dr.eigval = eig(transition_matrix(dr));
         dr.rank = 0;
         if any(abs(dr.eigval) > options_.qz_criterium)
@@ -239,6 +258,54 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
     else
         aa = jacobia_;
     end
+
+    % If required, try Gary Anderson and G Moore AIM solver, that is, if not check
+    % only and if 1st order (added by GP July'08)
+    if (options_.useAIM == 1) && (task == 0) && (options_.order == 1) 
+        try
+            [dr,aimcode]=dynAIMsolver1(aa,M_,dr);
+
+            % reuse some of the bypassed code and tests that may be needed 
+            
+            if aimcode ~=1
+                info(1) = aimcode;
+                info(2) = 1.0e+8;
+                return
+            end
+            [A,B] =transition_matrix(dr);
+            dr.eigval = eig(A);
+%            if any(abs(dr.eigval) > options_.qz_criterium)
+%                temp = sort(abs(dr.eigval));
+%                nba = nnz(abs(dr.eigval) > options_.qz_criterium);
+%                temp = temp(nd-nba+1:nd)-1-options_.qz_criterium;
+%                info(1) = 3;
+%                info(2) = temp'*temp;
+%                return
+%            end
+            sdim = sum( abs(dr.eigval) < options_.qz_criterium );
+            nba = nd-sdim;
+
+            nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1);
+            if nba ~= nyf
+                temp = sort(abs(dr.eigval));
+                if nba > nyf
+                    temp = temp(nd-nba+1:nd-nyf)-1-options_.qz_criterium;
+                    info(1) = 3
+                elseif nba < nyf;
+                    temp = temp(nd-nyf+1:nd-nba)-1-options_.qz_criterium;
+                    info(1) = 4
+                end
+                info(2) = temp'*temp;
+                return
+            end
+        catch
+            %warning('Problem with using AIM solver - Using Dynare solver instead');
+            disp('Problem with using AIM solver - Using Dynare solver instead');
+            options_.useAIM = 0; % and then try mjdgges instead
+        end
+    end % end if useAIM and...
+    %else  % use original Dynare solver
+    if  ~((options_.useAIM == 1)&& (task == 0) && (options_.order == 1)) % || isempty(options_.useAIM)  
         a = aa(:,nonzeros(k1'));
         b = aa(:,k0);
         b10 = b(1:nstatic,1:nstatic);
@@ -270,6 +337,12 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
             end
         end
 
+        % 1) if mjdgges.dll (or .mexw32 or ....) doesn't exit, 
+        % matlab/qz is added to the path. There exists now qz/mjdgges.m that 
+        % contains the calls to the old Sims code 
+        % 2) In  global_initialization.m, if mjdgges.m is visible exist(...)==2, 
+        % this means that the DLL isn't avaiable and use_qzdiv is set to 1
+
         [ss,tt,w,sdim,dr.eigval,info1] = mjdgges(e,d,options_.qz_criterium);
 
         if info1
@@ -296,7 +369,7 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
             temp = sort(abs(dr.eigval));
             if nba > nyf
                 temp = temp(nd-nba+1:nd-nyf)-1-options_.qz_criterium;
-            info(1) = 3;
+                info(1) = 3
             elseif nba < nyf;
                 temp = temp(nd-nyf+1:nd-nba)-1-options_.qz_criterium;
                 info(1) = 4;
@@ -354,6 +427,8 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
             dr.ghx = [temp; dr.ghx];
             temp = [];
         end
+    end % if not use AIM and ....
+    % End of if... and if not... main AIM Blocks, continue as per usual...
     
     if options_.loglinear == 1
         k = find(dr.kstate(:,2) <= M_.maximum_endo_lag+1);
@@ -365,6 +440,7 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
         dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
     end
     
+    if  ~((options_.useAIM == 1) && (options_.order == 1)) % if not use AIM ...
         %% Necessary when using Sims' routines for QZ
         if options_.use_qzdiv
             gx = real(gx);
@@ -372,6 +448,7 @@ function [dr,info,M_,options_,oo_] = dr1(dr,task,M_,options_,oo_)
             dr.ghx = real(dr.ghx);
             dr.ghu = real(dr.ghu);
         end
+    end % if not use AIM 
     
     %exogenous deterministic variables
     if M_.exo_det_nbr > 0