Commit 69fc0e32 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

New oct-file "qzcomplex" for bringing the complex QZ decomposition to Octave....

New oct-file "qzcomplex" for bringing the complex QZ decomposition to Octave. Fixes issues with partial information under Octave.
parent 8f3d84ed
......@@ -84,6 +84,7 @@ ylwrap
# DLL rules
*.mex
*.dll
*.oct
*.mexglx
*.mexa64
*.mexw32
......
......@@ -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
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ logposterior qzcomplex
if HAVE_GSL
SUBDIRS += swz
endif
......
EXEXET = .mex
include ../mex.am
include ../../bytecode.am
......@@ -101,6 +101,7 @@ AC_CONFIG_FILES([Makefile
k_order_perturbation/Makefile
dynare_simul_/Makefile
swz/Makefile
logposterior/Makefile])
logposterior/Makefile
qzcomplex/Makefile])
AC_OUTPUT
EXEEXT = .mex
include ../mex.am
include ../../dynare_simul_.am
EXEEXT = .mex
include ../mex.am
include ../../gensylv.am
EXEEXT = .mex
include ../mex.am
include ../../k_order_perturbation.am
EXEEXT = .mex
include ../mex.am
include ../../kronecker.am
EXEEXT = .mex
include ../mex.am
include ../../libdynare++.am
EXEEXT = .mex
include ../mex.am
include ../../logposterior.am
EXEEXT = .mex
CPPFLAGS += $(shell $(MKOCTFILE) -p CPPFLAGS)
CPPFLAGS += $(shell $(MKOCTFILE) -p INCFLAGS)
CPPFLAGS += -I$(top_srcdir)/../../sources
......
EXEEXT = .mex
include ../mex.am
include ../../mjdgges.am
EXEEXT = .oct
include ../mex.am
vpath %.cc $(top_srcdir)/../../sources/qzcomplex
noinst_PROGRAMS = qzcomplex
nodist_qzcomplex_SOURCES = qzcomplex.cc
EXEEXT = .mex
include ../mex.am
include ../../swz.am
......@@ -7,6 +7,7 @@ EXTRA_DIST = \
mjdgges \
kronecker \
bytecode \
qzcomplex \
k_order_perturbation
clean-local:
......
/*
* Oct-file for bringing the complex QZ decomposition to Octave.
* Simple wrapper around LAPACK's zgges.
*
* 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(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)
{
error("qzcomplex: needs two input arguments");
return retval;
}
if (nargout != 4)
{
error("qzcomplex: needs four output arguments");
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;
}
......@@ -110,7 +110,7 @@ SectionGroupEnd
Section "MEX files for Octave 3.2.4 (MinGW build)"
SetOutPath $INSTDIR\mex\octave
File ..\mex\octave\*.mex
File ..\mex\octave\*.mex ..\mex\octave\*.oct
SectionEnd
Section "Dynare++ (standalone executable)"
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment