From 81823ad0353a810c0fa500c62ddb3b335e27ea76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 21 Oct 2010 15:43:13 +0200 Subject: [PATCH] New oct-file for "ordschur": the diffuse filter now works under Octave --- matlab/dynare_config.m | 6 -- matlab/missing/ordschur/ordschur.m | 20 ----- mex/build/octave/Makefile.am | 2 +- mex/build/octave/configure.ac | 3 +- mex/build/octave/ordschur/Makefile.am | 8 ++ mex/sources/Makefile.am | 3 +- mex/sources/ordschur/ordschur.cc | 107 ++++++++++++++++++++++++++ tests/Makefile.am | 24 +++--- 8 files changed, 132 insertions(+), 41 deletions(-) delete mode 100644 matlab/missing/ordschur/ordschur.m create mode 100644 mex/build/octave/ordschur/Makefile.am create mode 100644 mex/sources/ordschur/ordschur.cc diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index 878160ab9..eee684ef9 100644 --- a/matlab/dynare_config.m +++ b/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']) diff --git a/matlab/missing/ordschur/ordschur.m b/matlab/missing/ordschur/ordschur.m deleted file mode 100644 index c9ba12c8f..000000000 --- a/matlab/missing/ordschur/ordschur.m +++ /dev/null @@ -1,20 +0,0 @@ -## 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 . - -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 diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index 6f7a3add2..0471b4409 100644 --- a/mex/build/octave/Makefile.am +++ b/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 diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac index e54aa85b0..a30d348fa 100644 --- a/mex/build/octave/configure.ac +++ b/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 diff --git a/mex/build/octave/ordschur/Makefile.am b/mex/build/octave/ordschur/Makefile.am new file mode 100644 index 000000000..e9225d7ec --- /dev/null +++ b/mex/build/octave/ordschur/Makefile.am @@ -0,0 +1,8 @@ +EXEEXT = .oct +include ../mex.am + +vpath %.cc $(top_srcdir)/../../sources/ordschur + +noinst_PROGRAMS = ordschur + +nodist_ordschur_SOURCES = ordschur.cc diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am index ba2152c08..c6fdccfcb 100644 --- a/mex/sources/Makefile.am +++ b/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` diff --git a/mex/sources/ordschur/ordschur.cc b/mex/sources/ordschur/ordschur.cc new file mode 100644 index 000000000..5f7767c1b --- /dev/null +++ b/mex/sources/ordschur/ordschur.cc @@ -0,0 +1,107 @@ +/* + * 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 . + */ + +/* + * 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 . + */ + +#include +#include + +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; +} diff --git a/tests/Makefile.am b/tests/Makefile.am index 1b7514d67..4b5f0781b 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,6 @@ 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