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

Remove fallbacks for ilu, linsolve and ordschur under Octave

These functions are all included in Octave 4.2.
parent a3301777
......@@ -103,11 +103,6 @@ if isoctave && ~compare_versions(version(), supported_octave_version(),'>=')
skipline()
end
% ilu is missing in Octave < 4.0
if isoctave && octave_ver_less_than('4.0')
p{end+1} = '/missing/ilu';
end
% corrcoef with two outputs is missing in Octave (ticket #796)
if isoctave && ~user_has_octave_forge_package('nan')
p{end+1} = '/missing/corrcoef';
......
function [L, U, P] = ilu(A, setup)
% Partially implement function ilu using the (now deprecated) luinc
% Copyright (C) 2013-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/>.
if nargout ~= 2
error('Only two output arguments supported')
end
if nargin ~= 2
error('Only two input arguments supported')
end
if isfield(setup, 'milu')
if setup.milu == 'off'
setup.milu = 0;
else
error(['Unsupported milu: ' setup.milu ])
end
end
[L, U] = luinc(A, setup);
end
......@@ -4,14 +4,6 @@ ACLOCAL_AMFLAGS = -I ../../../m4
if DO_SOMETHING
SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv qzcomplex block_kalman_filter sobol local_state_space_iterations
if COMPILE_LINSOLVE
SUBDIRS += linsolve
endif
if COMPILE_ORDSCHUR
SUBDIRS += ordschur
endif
if HAVE_MATIO
SUBDIRS += k_order_perturbation dynare_simul_
endif
......
......@@ -33,13 +33,8 @@ if test "x$MKOCTFILE" != "x"; then
LDFLAGS="`$MKOCTFILE -p LFLAGS` `$MKOCTFILE -p LDFLAGS`"
OCTAVE_VERSION=`$MKOCTFILE -v 2>&1 | sed 's/mkoctfile, version //'`
AX_COMPARE_VERSION([$OCTAVE_VERSION], [lt], [3.6], [AC_MSG_ERROR([Your Octave is too old, please upgrade to version 3.6 at least.])])
AX_COMPARE_VERSION([$OCTAVE_VERSION], [ge], [3.8], [OCTAVE38=yes])
AX_COMPARE_VERSION([$OCTAVE_VERSION], [ge], [4.0], [OCTAVE40=yes])
fi
AM_CONDITIONAL([COMPILE_LINSOLVE], [test "$OCTAVE38" != "yes"])
AM_CONDITIONAL([COMPILE_ORDSCHUR], [test "$OCTAVE40" != "yes"])
CFLAGS="$CFLAGS -Wall -Wno-parentheses"
FFLAGS="$FFLAGS -Wall"
CXXFLAGS="$CXXFLAGS -Wall -Wno-parentheses"
......@@ -109,18 +104,6 @@ else
BUILD_MS_SBVAR_MEX_OCTAVE="no (missing GSL or MatIO library)"
fi
if test -n "$MKOCTFILE" -a "$OCTAVE38" != "yes"; then
BUILD_LINSOLVE_OCTAVE="yes"
else
BUILD_LINSOLVE_OCTAVE="no (Octave >= 3.8)"
fi
if test -n "$MKOCTFILE" -a "$OCTAVE40" != "yes"; then
BUILD_ORDSCHUR_OCTAVE="yes"
else
BUILD_ORDSCHUR_OCTAVE="no (Octave >= 4.0)"
fi
AC_ARG_ENABLE([openmp], AS_HELP_STRING([--enable-openmp], [use OpenMP for parallelization of some MEX files]), [
if test "x$enable_openmp" = "xyes"; then
CPPFLAGS="$CPPFLAGS -DUSE_OMP"
......@@ -138,8 +121,6 @@ Binaries (with "make"):
MS-SBVAR MEX files for Octave: $BUILD_MS_SBVAR_MEX_OCTAVE
Kalman Steady State MEX file for Octave: $BUILD_KALMAN_STEADY_STATE_OCTAVE
K-order and dynare_simul MEX for Octave: $BUILD_KORDER_DYNSIMUL_MEX_OCTAVE
Linsolve for Octave: $BUILD_LINSOLVE_OCTAVE
Ordschur for Octave: $BUILD_ORDSCHUR_OCTAVE
])
......@@ -152,12 +133,10 @@ AC_CONFIG_FILES([Makefile
k_order_perturbation/Makefile
dynare_simul_/Makefile
qzcomplex/Makefile
ordschur/Makefile
kalman_steady_state/Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
local_state_space_iterations/Makefile
linsolve/Makefile])
local_state_space_iterations/Makefile])
AC_OUTPUT
EXEEXT = .oct
include ../mex.am
mex_PROGRAMS = linsolve
nodist_linsolve_SOURCES = $(top_srcdir)/../../sources/linsolve/linsolve.cc
EXEEXT = .oct
include ../mex.am
mex_PROGRAMS = ordschur
nodist_ordschur_SOURCES = $(top_srcdir)/../../sources/ordschur/ordschur.cc
......@@ -8,13 +8,11 @@ EXTRA_DIST = \
bytecode \
qzcomplex \
k_order_perturbation \
ordschur \
kalman_steady_state \
ms-sbvar \
block_kalman_filter \
sobol \
local_state_space_iterations \
linsolve
local_state_space_iterations
clean-local:
rm -rf `find mex/sources -name *.o`
......
/*
* 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 Sébastien Villemot <sebastien@dynare.org>.
*/
/*
* Copyright (C) 2012-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/>.
*/
#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\
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)
{
if (dimB(1) > 0)
retval(1) = rcond;
else // If B has zero columns, A.solve() apparently does not compute A's rcond
retval(1) = A.rcond();
}
return retval;
}
/*
* 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 Sébastien Villemot <sebastien@dynare.org>.
*/
/*
* Copyright (C) 2010-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/f77-fcn.h>
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("ordschur: dtrsen failed");
return retval;
}
retval(0) = octave_value(U);
retval(1) = octave_value(T);
return retval;
}
/*
%!test
%! A = [1 2 3 -2; 4 5 6 -5 ; 7 8 9 -5; 10 11 12 4 ];
%! [U, T] = schur(A);
%! [US, TS] = ordschur(U, T, [ 0 0 1 1 ]);
%! assert(US*TS*US', A, sqrt(eps))
*/
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