diff --git a/matlab/partial_information/PI_gensys.m b/matlab/partial_information/PI_gensys.m index 589e39dc1b33107ca9512f94a2c7257dc8633cdb..8e905b8ff18dbd27ac30cac2fad5fcc8b6692e20 100644 --- a/matlab/partial_information/PI_gensys.m +++ b/matlab/partial_information/PI_gensys.m @@ -176,13 +176,9 @@ end G0pi=eye(n+FL_RANK+NX); try - % In Matlab: [aa bb q z v w] = qz(a,b) s.t. qaz = aa, qbz = bb % - % In Octave: [aa bb q z v w] = qz(a,b) s.t. q'az = aa, q'bz=bb % - % and qzcomplex() extension based on lapack zgges produces same - % qz output for Octave as Matlab qz() does for Matlab thus: if isoctave - [a b q z]=qzcomplex(G0pi,G1pi); - q=q'; + % Need to force QZ complex on Octave (otherwise it returns the real one) + [a b q z]=qz(complex(G0pi),complex(G1pi)); else [a b q z]=qz(G0pi,G1pi); end diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index 9544e0e52367ea56d15e4cf82ab201336150b92b..64bb633ba3ac4b526ab9c14bc66f398fb5664548 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 qzcomplex block_kalman_filter sobol local_state_space_iterations +SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv block_kalman_filter sobol local_state_space_iterations if HAVE_MATIO SUBDIRS += k_order_perturbation dynare_simul_ diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac index 8c664c6e2d1471a62a138aae0d23c3c1c846c41d..d3661feb6db9f77772ad3c56eddefa4fe8fc7c41 100755 --- a/mex/build/octave/configure.ac +++ b/mex/build/octave/configure.ac @@ -132,7 +132,6 @@ AC_CONFIG_FILES([Makefile gensylv/Makefile k_order_perturbation/Makefile dynare_simul_/Makefile - qzcomplex/Makefile kalman_steady_state/Makefile ms_sbvar/Makefile block_kalman_filter/Makefile diff --git a/mex/build/octave/qzcomplex/Makefile.am b/mex/build/octave/qzcomplex/Makefile.am deleted file mode 100644 index 252686edcf45920945792d308404f19cc44d9c8e..0000000000000000000000000000000000000000 --- a/mex/build/octave/qzcomplex/Makefile.am +++ /dev/null @@ -1,6 +0,0 @@ -EXEEXT = .oct -include ../mex.am - -mex_PROGRAMS = qzcomplex - -nodist_qzcomplex_SOURCES = $(top_srcdir)/../../sources/qzcomplex/qzcomplex.cc diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am index 6f201c506d7d509192691f98abdb36c425f5adb4..4d7007ee151a0d1633fc1468c82feb8b798ff3ec 100644 --- a/mex/sources/Makefile.am +++ b/mex/sources/Makefile.am @@ -6,7 +6,6 @@ EXTRA_DIST = \ mjdgges \ kronecker \ bytecode \ - qzcomplex \ k_order_perturbation \ kalman_steady_state \ ms-sbvar \ diff --git a/mex/sources/qzcomplex/qzcomplex.cc b/mex/sources/qzcomplex/qzcomplex.cc deleted file mode 100644 index 74773fcdc59dbdb8e9a8e85aac8742570090882b..0000000000000000000000000000000000000000 --- a/mex/sources/qzcomplex/qzcomplex.cc +++ /dev/null @@ -1,122 +0,0 @@ -/* - * Oct-file for bringing the complex QZ decomposition to Octave. - * Simple wrapper around LAPACK's zgges. - * - * Written by Sébastien Villemot <sebastien@dynare.org>. - */ - -/* - * Copyright (C) 2010-2017 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/>. - */ - -#ifdef __MINGW32__ -# define __CROSS_COMPILATION__ -#endif - -#ifdef __MINGW64__ -# define __CROSS_COMPILATION__ -#endif - -#ifdef __CROSS_COMPILATION__ -# define M_PI 3.14159265358979323846 -#endif - -#include <octave/oct.h> -#include <octave/f77-fcn.h> - -extern "C" -{ - F77_RET_T - F77_FUNC(zgges, ZGGES) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - octave_idx_type (*)(Complex *, Complex *), const octave_idx_type &, - Complex *, const octave_idx_type &, Complex *, const octave_idx_type &, - octave_idx_type &, Complex *, Complex *, Complex *, const octave_idx_type &, - Complex *, const octave_idx_type &, Complex *, const octave_idx_type &, - double *, octave_idx_type *, octave_idx_type &); -} - -DEFUN_DLD(qzcomplex, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} [ @var{aa}, @var{bb}, @var{q}, @var{z} ] = qzcomplex (@var{a}, @var{b})\n\ -\n\ -Computes the complex QZ decomposition of @math{(A, B)}, satisfying:\n\ -@example\n\ - AA = Q'*A*Z, BB = Q'*B*Z\n\ -@end example\n\ -@end deftypefn\n\ -") -{ - int nargin = args.length(); - octave_value_list retval; - - if (nargin != 2 || nargout != 4) - { - print_usage(); - return retval; - } - - ComplexMatrix A = args(0).complex_matrix_value(); - ComplexMatrix B = args(1).complex_matrix_value(); - if (error_state) - return retval; - - dim_vector dimA = A.dims(); - dim_vector dimB = B.dims(); - octave_idx_type n = dimA(0); - if (n != dimA(1) || n != dimB(0) || n != dimB(1)) - { - error("qzcomplex: input matrices must be square and of same size"); - return retval; - } - - octave_idx_type lwork = 2*n; - OCTAVE_LOCAL_BUFFER(Complex, alpha, n); - OCTAVE_LOCAL_BUFFER(Complex, beta, n); - OCTAVE_LOCAL_BUFFER(Complex, work, lwork); - OCTAVE_LOCAL_BUFFER(double, rwork, 8*n); - ComplexMatrix vsl(n, n), vsr(n, n); - octave_idx_type sdim, info; - - F77_XFCN(zgges, ZGGES, (F77_CONST_CHAR_ARG("V"), F77_CONST_CHAR_ARG("V"), - F77_CONST_CHAR_ARG("N"), NULL, - n, A.fortran_vec(), n, B.fortran_vec(), n, sdim, - alpha, beta, vsl.fortran_vec(), n, vsr.fortran_vec(), n, - work, lwork, rwork, NULL, info)); - - if (info != 0) - { - error("qzcomplex: zgges failed"); - return retval; - } - - retval(0) = octave_value(A); - retval(1) = octave_value(B); - retval(2) = octave_value(vsl); - retval(3) = octave_value(vsr); - return retval; -} - -/* - - %!test - %! A = [ 1 2 3+1i; 4 5-1i 0; -1 -5 3]; - %! B = [ -2 -8i 4; 1 5+3i 5; 7 -10 -2]; - %! [AA,BB,Q,Z] = qzcomplex(A,B); - %! assert(Q'*A*Z, AA, sqrt(eps)) - %! assert(Q'*B*Z, BB, sqrt(eps)) - -*/