From a22d1d415adac0ce3301540777ee46446d7db232 Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@mjui.fr>
Date: Mon, 10 Sep 2012 13:25:31 +0200
Subject: [PATCH] replaced rank() by rcond() in evaluating whether Z22 is full
 rank in checking Blanchard and Kahn conditions with CHECK

---
 matlab/check.m                  |  4 ++--
 matlab/dr_block.m               | 23 +++++++++++++++++------
 matlab/dyn_first_order_solver.m |  6 +++++-
 3 files changed, 24 insertions(+), 9 deletions(-)

diff --git a/matlab/check.m b/matlab/check.m
index 6014bfdf1..1eeb7c788 100644
--- a/matlab/check.m
+++ b/matlab/check.m
@@ -86,7 +86,7 @@ end;
 n_explod = nnz(abs(eigenvalues_) > options.qz_criterium);
 
 result = 0;
-if (nyf== n_explod) && (dr.rank == nyf)
+if (nyf== n_explod) && (dr.full_rank)
     result = 1;
 end
 
@@ -99,7 +99,7 @@ if options.noprint == 0
     disp(sprintf('\nThere are %d eigenvalue(s) larger than 1 in modulus ', n_explod));
     disp(sprintf('for %d forward-looking variable(s)',nyf));
     disp(' ')
-    if dr.rank == nyf && nyf == n_explod
+    if result
         disp('The rank condition is verified.')
     else
         disp('The rank conditions ISN''T verified!')
diff --git a/matlab/dr_block.m b/matlab/dr_block.m
index 44b95efcb..d6dca0a1c 100644
--- a/matlab/dr_block.m
+++ b/matlab/dr_block.m
@@ -72,7 +72,7 @@ else
     chck = 0;
 end;
 mexErrCheck('bytecode', chck);
-dr.rank = 0;
+dr.full_rank = 1;
 dr.eigval = [];
 dr.nstatic = 0;
 dr.nfwrd = 0;
@@ -151,7 +151,6 @@ for i = 1:Size;
             data(i).rank = 0;
         end
         dr.eigval = [dr.eigval ; data(i).eigval];
-        dr.rank = dr.rank + data(i).rank;
         %First order approximation
         if task ~= 1
             [tmp1, tmp2, indx_c] = find(M_.block_structure.block(i).lead_lag_incidence(2,:));
@@ -221,12 +220,14 @@ for i = 1:Size;
             indx_c = M_.block_structure.block(i).lead_lag_incidence(3,indx_r);
             data(i).eigval = 1 ./ diag(jacob(indx_r, indx_c));
             data(i).rank = sum(abs(data(i).eigval) > 0);
+            full_rank = (rcond(jacob(indx_r, indx_c)) > 1e-9);
         else
             data(i).eigval = [];
             data(i).rank = 0;
+            full_rank = 1;
         end
         dr.eigval = [dr.eigval ; data(i).eigval];
-        dr.rank = dr.rank + data(i).rank;
+        dr.full_rank = dr.full_rank && full_rank;
         %First order approximation
         if task ~= 1
             if (maximum_lag > 0)
@@ -310,11 +311,13 @@ for i = 1:Size;
         if maximum_lead > 0 && n_fwrd > 0
             data(i).eigval = - jacob(1 , n_pred + n - n_fwrd + 1 : n_pred + n) / jacob(1 , n_pred + n + 1 : n_pred + n + n_fwrd) ;
             data(i).rank = sum(abs(data(i).eigval) > 0);
+            full_rank = (abs(jacob(1,n_pred+n+1: n_pred_n+n_fwrd)) > 1e-9);
         else
             data(i).eigval = [];
             data(i).rank = 0;
+            full_rank = 1;
         end;
-        dr.rank = dr.rank + data(i).rank;
+        dr.full_rank = dr.full_rank && full_rank;
         dr.eigval = [dr.eigval ; data(i).eigval];
       case 6
       %% ------------------------------------------------------------------
@@ -323,11 +326,15 @@ for i = 1:Size;
             data(i).eigval = eig(- jacob(: , 1 : n_pred) / ...
                                  jacob(: , (n_pred + n_static + 1 : n_pred + n_static + n_pred )));
             data(i).rank = 0;
+            full_rank = (rcond(jacob(: , (n_pred + n_static + 1 : n_pred ...
+                                          + n_static + n_pred ))) > 1e-9);
         else
             data(i).eigval = [];
             data(i).rank = 0;
+            full_rank = 1;
         end;
         dr.eigval = [dr.eigval ; data(i).eigval];
+        dr.full_rank = dr.full_rank && full_rank;
         if task ~= 1
             if (maximum_lag > 0)
                  ghx = - jacob(: , 1 : n_pred) / jacob(: , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both);
@@ -390,11 +397,14 @@ for i = 1:Size;
             data(i).eigval = eig(- jacob(: , n_pred + n - n_fwrd + 1: n_pred + n))/ ...
                 jacob(: , n_pred + n + 1 : n_pred + n + n_fwrd);
             data(i).rank = sum(abs(data(i).eigval) > 0);
+            full_rank = (rcond(jacob(: , n_pred + n + 1 : n_pred + n + ...
+                                     n_fwrd)) > 1e-9);
         else
             data(i).eigval = [];
             data(i).rank = 0;
+            full_rank = 1;
         end;
-        dr.rank = dr.rank + data(i).rank;
+        dr.full_rank = dr.full_rank && full_rank;
         dr.eigval = [dr.eigval ; data(i).eigval];
       case {5,8}
       %% ------------------------------------------------------------------
@@ -450,7 +460,8 @@ for i = 1:Size;
             nba = nd-sdim;
             if task == 1
                 data(i).rank = rank(w(nd-nyf+1:end,nd-nyf+1:end));
-                dr.rank = dr.rank + data(i).rank;
+                dr.full_rank = dr.full_rank && (rcond(w(nd-nyf+1:end,nd- ...
+                                                        nyf+1:end)) > 1e-9);
                 if ~exist('OCTAVE_VERSION','builtin')
                     data(i).eigval = eig(E,D);
                 end
diff --git a/matlab/dyn_first_order_solver.m b/matlab/dyn_first_order_solver.m
index 8ad529c9e..b0d9a81c0 100644
--- a/matlab/dyn_first_order_solver.m
+++ b/matlab/dyn_first_order_solver.m
@@ -233,7 +233,11 @@ else
     nba = nd-sdim;
 
     if task == 1
-        dr.rank = rank(w(npred+nboth+1:end,npred+nboth+1:end));
+        if rcond(w(npred+nboth+1:end,npred+nboth+1:end)) < 1e-9
+            dr.full_rank = 0;
+        else
+            dr.full_rank = 1;
+        end            
         % Under Octave, eig(A,B) doesn't exist, and
         % lambda = qz(A,B) won't return infinite eigenvalues
         if ~exist('OCTAVE_VERSION')