New oct-file for "ordschur": the diffuse filter now works under Octave

matlab/dynare_config.m

@@ -71,12 +71,6 @@ if exist('OCTAVE_VERSION') && octave_ver_less_than('3.2.0')
     addpath([dynareroot '/missing/bicgstab'])
 end
 
-% orschur() is missing in Octave; we don't have a real replacement;
-% the one we provide just exits with an error message
-if exist('OCTAVE_VERSION')
-    addpath([dynareroot '/missing/ordschur'])
-end
-
 % bsxfun is missing in old versions of matlab (octave?)
 if ~exist('OCTAVE_VERSION') && matlab_ver_less_than('7.4')
     addpath([dynareroot '/missing/bsxfun'])

matlab/missing/ordschur/ordschur.m

-## Copyright (C) 2009 Dynare Team
-##
-## This file is part of Dynare.
-##
-## Dynare 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.
-##
-## Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.
-
-function [US,TS] = ordschur(U,T,SELECT)
-  error("The ordschur() function is missing under Octave. For this reason, Dynare can't estimate models with unit roots with Octave.")
-endfunction

mex/build/octave/Makefile.am

@@ -2,7 +2,7 @@ ACLOCAL_AMFLAGS = -I ../../../m4
 
 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
 if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ logposterior qzcomplex
+SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ logposterior qzcomplex ordschur
 if HAVE_GSL
 SUBDIRS += swz
 endif

mex/build/octave/configure.ac

@@ -102,6 +102,7 @@ AC_CONFIG_FILES([Makefile
                  dynare_simul_/Makefile
                  swz/Makefile
                  logposterior/Makefile
-                 qzcomplex/Makefile])
+                 qzcomplex/Makefile
+                 ordschur/Makefile])
 
 AC_OUTPUT

mex/build/octave/ordschur/Makefile.am

+EXEEXT = .oct
+include ../mex.am
+
+vpath %.cc $(top_srcdir)/../../sources/ordschur
+
+noinst_PROGRAMS = ordschur
+
+nodist_ordschur_SOURCES = ordschur.cc

mex/sources/Makefile.am

@@ -8,7 +8,8 @@ EXTRA_DIST = \
 	kronecker \
 	bytecode \
 	qzcomplex \
-	k_order_perturbation
+	k_order_perturbation \
+	ordschur
 
 clean-local:
 	rm -rf `find mex/sources -name *.o`

mex/sources/ordschur/ordschur.cc

+/*
+ * Oct-file for bringing MATLAB's ordschur function to Octave.
+ * Simple wrapper around LAPACK's dtrsen.
+ * Only supports real (double precision) decomposition.
+ * Only selection of eigenvalues with a boolean vector is supported.
+ *
+ * Written by Sebastien Villemot <sebastien.villemot@ens.fr>.
+ */
+
+/*
+ * Copyright (C) 2010 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare 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.
+ *
+ * Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <octave/oct.h>
+#include <octave/f77-fcn.h>
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC(dtrsen, DTRSEN) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
+                            const octave_idx_type *, const octave_idx_type &,
+                            double *, const octave_idx_type &, double *, const octave_idx_type &,
+                            double *, double *, octave_idx_type &, double &, double &, double *,
+                            const octave_idx_type &, octave_idx_type *,
+                            const octave_idx_type &, octave_idx_type &);
+}
+
+DEFUN_DLD (ordschur, args, nargout, "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} [ @var{us}, @var{ts} ] = ordschur (@var{u}, @var{t}, @var{select})\n\
+\n\
+Reorders the real Schur factorization @math{X = U*T*U'} so that selected\n\
+eigenvalues appear in the upper left diagonal blocks of the quasi triangular\n\
+Schur matrix @math{T}. The logical vector @var{select} specifies the selected\n\
+eigenvalues as they appear along @math{T}'s diagonal.\n\
+@end deftypefn\n\
+")
+{
+  int nargin = args.length();
+  octave_value_list retval;
+
+  if (nargin != 3 || nargout != 2)
+    {
+      print_usage();
+      return retval;
+    }
+
+  Matrix U = args(0).matrix_value();
+  Matrix T = args(1).matrix_value();
+  boolNDArray S = args(2).bool_array_value();
+  if (error_state)
+    return retval;
+
+  dim_vector dimU = U.dims();
+  dim_vector dimT = T.dims();
+  octave_idx_type n = dimU(0);
+  if (n != dimU(1) || n != dimT(0) || n != dimT(1))
+    {
+      error("ordschur: input matrices must be square and of same size");
+      return retval;
+    }
+  if (S.nelem() != n)
+    {
+      error("ordschur: selection vector has wrong size");
+      return retval;
+    }
+
+  octave_idx_type lwork = n, liwork = n;
+  OCTAVE_LOCAL_BUFFER(double, wr, n);
+  OCTAVE_LOCAL_BUFFER(double, wi, n);
+  OCTAVE_LOCAL_BUFFER(double, work, lwork);
+  OCTAVE_LOCAL_BUFFER(octave_idx_type, iwork, liwork);
+  octave_idx_type m, info;
+  double cond1, cond2;
+  OCTAVE_LOCAL_BUFFER(octave_idx_type, S2, n);
+  for (int i = 0; i < n; i++)
+    S2[i] = S(i);
+  
+  F77_XFCN (dtrsen, dtrsen, (F77_CONST_CHAR_ARG("N"), F77_CONST_CHAR_ARG("V"),
+                             S2, n, T.fortran_vec(), n, U.fortran_vec(), n,
+                             wr, wi, m, cond1, cond2, work, lwork, 
+                             iwork, liwork, info));
+
+  if (info != 0)
+    {
+      error("qzcomplex: zgges failed");
+      return retval;
+    }
+
+  retval(0) = octave_value(U);
+  retval(1) = octave_value(T);
+  return retval;
+}

tests/Makefile.am

 DYNARE_ROOT = $(abs_top_srcdir)/matlab
 
-# Under Octave we only test a subset of MOD files, because of missing features (models with unit roots, reading Excel files)
+# Under Octave we only test a subset of MOD files, because of missing features
 OCTAVE_MODS = \
 	ramst.mod \
 	ramst_a.mod \
@@ -41,9 +41,14 @@ OCTAVE_MODS = \
 	partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod \
 	arima/mod1.mod \
 	arima/mod1a.mod \
+	arima/mod1b.mod \
+	arima/mod1c.mod \
 	arima/mod2.mod \
+	arima/mod2a.mod \
+	arima/mod2b.mod \
 	arima/mod2c.mod \
 	fs2000/fs2000.mod \
+	fs2000/fs2000a.mod \
 	fs2000/fs2000c.mod \
 	homotopy/homotopy1_test.mod \
 	homotopy/homotopy2_test.mod \
@@ -54,6 +59,12 @@ OCTAVE_MODS = \
 	AIM/fs2000x10L9_L_AIM.mod \
 	AIM/fs2000x10_L9_L.mod \
 	AIM/fs2000x10_L9_L_AIM.mod \
+	AIM/fs2000_b1L1L.mod \
+	AIM/fs2000_b1L1L_AIM.mod \
+	AIM/ls2003_2L0L.mod \
+	AIM/ls2003_2L0L_AIM.mod \
+	AIM/ls2003_2L2L.mod \
+	AIM/ls2003_2L2L_AIM.mod \
 	conditional_variance_decomposition/example1.mod \
 	dsge-var/simul_hybrid.mod \
 	dsge-var/dsgevar_forward_calibrated_lambda.mod \
@@ -66,17 +77,6 @@ OCTAVE_MODS = \
 	external_function/example1_no_deriv_functions_provided_dll.mod
 
 MODS = $(OCTAVE_MODS) \
-	arima/mod1b.mod \
-	arima/mod1c.mod \
-	arima/mod2a.mod \
-	arima/mod2b.mod \
-	fs2000/fs2000a.mod \
-	AIM/fs2000_b1L1L.mod \
-	AIM/fs2000_b1L1L_AIM.mod \
-	AIM/ls2003_2L0L.mod \
-	AIM/ls2003_2L0L_AIM.mod \
-	AIM/ls2003_2L2L.mod \
-	AIM/ls2003_2L2L_AIM.mod \
 	block_bytecode/fs2000_gmres.mod \
 	block_bytecode/ramst_a.mod