Provide a better implementation of linsolve for Octave

Closes: #273

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/']);

matlab/missing/linsolve/linsolve.m

-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 <http://www.gnu.org/licenses/>.
-
-    c = [];
-    x = [];
-    if nargin == 3
-        if isfield(opts,'TRANSA')
-            A = A';
-        end
-    end
-    if nargout == 2
-        c = rcond(A);
-    end
-    
-    x = A\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

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

mex/build/octave/linsolve/Makefile.am

+EXEEXT = .oct
+include ../mex.am
+
+noinst_PROGRAMS = linsolve
+
+nodist_linsolve_SOURCES = $(top_srcdir)/../../sources/linsolve/linsolve.cc

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`

mex/sources/linsolve/linsolve.cc

+/*
+ * 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 <sebastien.villemot@ens.fr>.
+ */
+
+/*
+ * 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 <http://www.gnu.org/licenses/>.
+ */
+
+
+#include <octave/oct.h>
+#include <octave/ov-struct.h>
+
+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;
+}