From 2e65a9ab96fdd0e3e67145f06e852a65a34ae742 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Mon, 6 Aug 2012 18:20:29 +0200 Subject: [PATCH] Provide a better implementation of linsolve for Octave Closes: #273 --- matlab/dynare_config.m | 5 -- matlab/missing/linsolve/linsolve.m | 33 ------- mex/build/octave/Makefile.am | 2 +- mex/build/octave/configure.ac | 3 +- mex/build/octave/linsolve/Makefile.am | 6 ++ mex/sources/Makefile.am | 3 +- mex/sources/linsolve/linsolve.cc | 120 ++++++++++++++++++++++++++ 7 files changed, 131 insertions(+), 41 deletions(-) delete mode 100644 matlab/missing/linsolve/linsolve.m create mode 100644 mex/build/octave/linsolve/Makefile.am create mode 100644 mex/sources/linsolve/linsolve.cc diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index 0644f8e7f..98ed0c0bb 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -93,11 +93,6 @@ if (exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('statistics')) ... addpath([dynareroot '/missing/nanmean']) end -% linsolve is missing in Octave -if (exist('OCTAVE_VERSION')) - addpath([dynareroot '/missing/linsolve']) -end - % Add path to MEX files if exist('OCTAVE_VERSION') addpath([dynareroot '../mex/octave/']); diff --git a/matlab/missing/linsolve/linsolve.m b/matlab/missing/linsolve/linsolve.m deleted file mode 100644 index 2ae167bd4..000000000 --- a/matlab/missing/linsolve/linsolve.m +++ /dev/null @@ -1,33 +0,0 @@ -function [x,c] = linsolve(A,B,opts) -% (very imperfect) Clone of Matlab's linsolve. - -% Copyright (C) 2010-2011 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 . - - c = []; - x = []; - if nargin == 3 - if isfield(opts,'TRANSA') - A = A'; - end - end - if nargout == 2 - c = rcond(A); - end - - x = A\B; - diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index 132ddbe15..2689afab4 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_ qzcomplex ordschur block_kalman_filter sobol local_state_space_iterations +SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv k_order_perturbation dynare_simul_ qzcomplex ordschur block_kalman_filter sobol local_state_space_iterations linsolve if HAVE_GSL if HAVE_MATIO diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac index da3c6e2ac..0a1362799 100644 --- a/mex/build/octave/configure.ac +++ b/mex/build/octave/configure.ac @@ -132,6 +132,7 @@ AC_CONFIG_FILES([Makefile ms_sbvar/Makefile block_kalman_filter/Makefile sobol/Makefile - local_state_space_iterations/Makefile]) + local_state_space_iterations/Makefile + linsolve/Makefile]) AC_OUTPUT diff --git a/mex/build/octave/linsolve/Makefile.am b/mex/build/octave/linsolve/Makefile.am new file mode 100644 index 000000000..4f3a5aaf1 --- /dev/null +++ b/mex/build/octave/linsolve/Makefile.am @@ -0,0 +1,6 @@ +EXEEXT = .oct +include ../mex.am + +noinst_PROGRAMS = linsolve + +nodist_linsolve_SOURCES = $(top_srcdir)/../../sources/linsolve/linsolve.cc diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am index 41b94ffa6..57934b164 100644 --- a/mex/sources/Makefile.am +++ b/mex/sources/Makefile.am @@ -15,7 +15,8 @@ EXTRA_DIST = \ ms-sbvar \ block_kalman_filter \ sobol \ - local_state_space_iterations + local_state_space_iterations \ + linsolve clean-local: rm -rf `find mex/sources -name *.o` diff --git a/mex/sources/linsolve/linsolve.cc b/mex/sources/linsolve/linsolve.cc new file mode 100644 index 000000000..e6329a28f --- /dev/null +++ b/mex/sources/linsolve/linsolve.cc @@ -0,0 +1,120 @@ +/* + * Oct-file for bringing MATLAB's linsolve function to Octave. + * + * The implementation is incomplete: + * - it only knows about the TRANSA, LT, UT and SYM options + * - it only works with square matrices + * - it only works on double matrices (no single precision or complex) + * + * Written by Sebastien Villemot . + */ + +/* + * Copyright (C) 2012 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 + +DEFUN_DLD(linsolve, args, nargout, "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} @var{x} = linsolve (@var{a}, @var{b})\n\ +@deftypefnx {Loadable Function} [ @var{x}, @var{r} ] = linsolve (@var{a}, @var{b})\n\ +@deftypefnx {Loadable Function} @var{x} = linsolve (@var{a}, @var{b}, @var{options})\n\ +@deftypefnx {Loadable Function} [ @var{x}, @var{r} ] = linsolve (@var{a}, @var{b}, @var{options})\n\ +\n\ +Solves the linear system @math{A*X = B} and returns @var{X}.\n\ +\n\ +Alternatively, if @var{options} is provided and has a field @code{TRANSA} equal \ +to @code{true}, then it solves the system @math{A'*X = B}.\n\ +\n\ +\n\Also, the @code{LT} field of @var{options} (resp. the @code{UT} field) can be set \ +to @code{true} to indicate that the matrix @var{a} is lower (resp. upper) \ +triangular; similarly, the @code{SYM} field can be set to @code{true} to \ +indicate that the matrix is symmetric.\n\ +\n\ +If requested, @var{r} will contain the reciprocal condition number.\n\ +@end deftypefn\n\ +") +{ + int nargin = args.length(); + octave_value_list retval; + + if (nargin > 3 || nargin < 2 || nargout > 2) + { + print_usage(); + return retval; + } + + Matrix A = args(0).matrix_value(); + Matrix B = args(1).matrix_value(); + if (error_state) + return retval; + + dim_vector dimA = A.dims(); + dim_vector dimB = B.dims(); + + if (dimA(0) != dimB(0)) + { + error("linsolve: must have same number of lines in A and B"); + return retval; + } + + if (dimA(0) != dimA(1)) + { + error("linsolve: rectangular A not yet supported"); + return retval; + } + + + MatrixType typA; + typA.mark_as_full(); + + bool transa = false; + if (nargin == 3) + { + octave_scalar_map opts = args(2).scalar_map_value(); + if (error_state) + return retval; + + octave_value tmp = opts.contents("TRANSA"); + transa = tmp.is_defined() && tmp.bool_matrix_value().elem(0); + + tmp = opts.contents("UT"); + if (tmp.is_defined() && tmp.bool_matrix_value().elem(0)) + typA.mark_as_upper_triangular(); + + tmp = opts.contents("LT"); + if (tmp.is_defined() && tmp.bool_matrix_value().elem(0)) + typA.mark_as_lower_triangular(); + + tmp = opts.contents("SYM"); + if (tmp.is_defined() && tmp.bool_matrix_value().elem(0)) + typA.mark_as_symmetric(); + } + + double rcond; + octave_idx_type info; + + retval(0) = A.solve(typA, B, info, rcond, NULL, true, transa ? blas_trans : blas_no_trans); + + if (nargout == 2) + retval(1) = rcond; + + return retval; +}