Verified Commit 6a269268 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Add new block_trust_region MEX

This MEX solves nonlinear systems of equations using a trust region algorithm.
The problem is subdivided in smaller problems by doing a block
triangularisation of the Jacobian at the guess value, using the
Dulmage-Mendelsohn algorithm.

The interface of the MEX is simply:

  [x, info] = block_trust_region(f, guess_value);

Where f is either a function handle or a string designating a function.
f must take one argument (the evaluation point), and return either one or two
arguments (the residuals and, optionally, the Jacobian).

On success, info=0; on failure, info=1.
parent 93bd817c
...@@ -101,6 +101,12 @@ doc/internals/ltxpng ...@@ -101,6 +101,12 @@ doc/internals/ltxpng
/mex/build/matlab/*/*.mod /mex/build/matlab/*/*.mod
/mex/build/octave/*/*.mod /mex/build/octave/*/*.mod
# Extra rules for trust_region MEX testfiles
/mex/sources/block_trust_region/test/*.mod
/mex/sources/block_trust_region/test/dulmage_mendelsohn_test
/mex/sources/block_trust_region/test/trust_region_test
!/mex/sources/block_trust_region/test/Makefile
# Dynare++ # Dynare++
/dynare++/integ/src/quadrature-points.dSYM/ /dynare++/integ/src/quadrature-points.dSYM/
/dynare++/src/dynare++.dSYM/ /dynare++/src/dynare++.dSYM/
......
mex_PROGRAMS = block_trust_region
nodist_block_trust_region_SOURCES = \
dulmage_mendelsohn.f08 \
matlab_fcn_closure.f08 \
trust_region.f08 \
mexFunction.f08 \
matlab_mex.F08 \
blas_lapack.F08
BUILT_SOURCES = $(nodist_block_trust_region_SOURCES)
CLEANFILES = $(nodist_block_trust_region_SOURCES)
dulmage_mendelsohn.mod: dulmage_mendelsohn.o
matlab_fcn_closure.mod: matlab_fcn_closure.o
matlab_fcn_closure.o: matlab_mex.mod
trust_region.mod: trust_region.o
trust_region.o: lapack.mod
mexFunction.o: matlab_mex.mod dulmage_mendelsohn.mod matlab_fcn_closure.mod trust_region.mod
%.f08: $(top_srcdir)/../../sources/block_trust_region/%.f08
$(LN_S) -f $< $@
ACLOCAL_AMFLAGS = -I ../../../m4 ACLOCAL_AMFLAGS = -I ../../../m4
SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem num_procs SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem num_procs block_trust_region
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_ # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if ENABLE_MEX_DYNAREPLUSPLUS if ENABLE_MEX_DYNAREPLUSPLUS
......
include ../mex.am
include ../../block_trust_region.am
...@@ -172,6 +172,7 @@ AC_CONFIG_FILES([Makefile ...@@ -172,6 +172,7 @@ AC_CONFIG_FILES([Makefile
sobol/Makefile sobol/Makefile
local_state_space_iterations/Makefile local_state_space_iterations/Makefile
perfect_foresight_problem/Makefile perfect_foresight_problem/Makefile
num_procs/Makefile]) num_procs/Makefile
block_trust_region/Makefile])
AC_OUTPUT AC_OUTPUT
ACLOCAL_AMFLAGS = -I ../../../m4 ACLOCAL_AMFLAGS = -I ../../../m4
SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem num_procs SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem num_procs block_trust_region
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_ # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if ENABLE_MEX_DYNAREPLUSPLUS if ENABLE_MEX_DYNAREPLUSPLUS
......
EXEEXT = .mex
include ../mex.am
include ../../block_trust_region.am
...@@ -129,6 +129,7 @@ AC_CONFIG_FILES([Makefile ...@@ -129,6 +129,7 @@ AC_CONFIG_FILES([Makefile
sobol/Makefile sobol/Makefile
local_state_space_iterations/Makefile local_state_space_iterations/Makefile
perfect_foresight_problem/Makefile perfect_foresight_problem/Makefile
num_procs/Makefile]) num_procs/Makefile
block_trust_region/Makefile])
AC_OUTPUT AC_OUTPUT
...@@ -15,7 +15,8 @@ EXTRA_DIST = \ ...@@ -15,7 +15,8 @@ EXTRA_DIST = \
gensylv \ gensylv \
dynare_simul_ \ dynare_simul_ \
perfect_foresight_problem \ perfect_foresight_problem \
num_procs num_procs \
block_trust_region
clean-local: clean-local:
rm -rf `find mex/sources -name *.o` rm -rf `find mex/sources -name *.o`
......
This diff is collapsed.
! Encapsulates a MATLAB function from ℝⁿ to ℝⁿ
! It must take as 1st argument the point where it is evaluated,
! and return 1 or 2 arguments: value and, optionally, the Jacobian
! The input is copied on entry, the output arguments are copied on exit.
! It may also take extra arguments, which are stored in the module (hence the
! latter is conceptually equivalent to a closure).
!
! Additionally, if x_indices and f_indices are associated, then the matlab_fcn
! procedure only exposes a restricted version of the MATLAB procedure, limited
! to the specified indices specified for x and f. In that case, x_all needs to
! be set in order to give the input values for the indices that are not passed
! to matlab_fcn.
! Copyright © 2019 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/>.
module matlab_fcn_closure
use iso_c_binding
use matlab_mex
implicit none
private
public :: func, extra_args, f_indices, x_indices, x_all, matlab_fcn
type(c_ptr), pointer :: func => null()
type(c_ptr), dimension(:), pointer :: extra_args => null()
integer, dimension(:), pointer :: f_indices => null(), x_indices => null()
real(real64), dimension(:), pointer :: x_all => null()
contains
subroutine matlab_fcn(x, fvec, fjac)
real(real64), dimension(:), intent(in) :: x
real(real64), dimension(size(x)), intent(out) :: fvec
real(real64), dimension(size(x), size(x)), intent(out), optional :: fjac
type(c_ptr), dimension(2) :: call_lhs
type(c_ptr), dimension(size(extra_args)+2) :: call_rhs
real(real64), dimension(:), pointer :: x_mat ! Needed to avoid gfortran ICE…
real(real64), dimension(:), pointer :: fvec_all
real(real64), dimension(:,:), pointer :: fjac_all
integer(c_int) :: nlhs
integer(mwSize) :: n_all
call_rhs(1) = func
if (associated(x_all)) then
n_all = size(x_all)
else
n_all = size(x)
end if
call_rhs(2) = mxCreateDoubleMatrix(n_all, 1_mwSize, mxREAL)
x_mat => mxGetPr(call_rhs(2))
if (associated(x_indices) .and. associated(x_all)) then
x_mat = x_all
x_mat(x_indices) = x
else
x_mat = x
end if
call_rhs(3:) = extra_args
if (present(fjac)) then
nlhs = 2
else
nlhs = 1
end if
! We use "feval", because it’s the only way of evaluating a function handle through mexCallMATLAB
if (mexCallMATLAB(nlhs, call_lhs, int(size(call_rhs), c_int), call_rhs, "feval") /= 0) &
call mexErrMsgTxt("Error calling function to be solved")
fvec_all => mxGetPr(call_lhs(1))
if (associated(f_indices)) then
fvec = fvec_all(f_indices)
else
fvec = fvec_all
end if
if (present(fjac)) then
fjac_all(1:n_all,1:n_all) => mxGetPr(call_lhs(2))
if (associated(x_indices) .and. associated(f_indices)) then
fjac = fjac_all(f_indices, x_indices)
else
fjac = fjac_all
end if
end if
end subroutine matlab_fcn
end module matlab_fcn_closure
! Copyright © 2019 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/>.
subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
use iso_fortran_env
use iso_c_binding
use dulmage_mendelsohn
use matlab_mex
use matlab_fcn_closure
use trust_region
implicit none
type(c_ptr), dimension(*), intent(in), target :: prhs
type(c_ptr), dimension(*), intent(out) :: plhs
integer(c_int), intent(in), value :: nlhs, nrhs
real(real64), dimension(:), allocatable, target :: x
type(dm_block), dimension(:), allocatable, target :: blocks
integer :: info, i
real(real64), parameter :: tolf = 1e-6_real64
real(real64), dimension(:), allocatable :: fvec
real(real64), dimension(:,:), allocatable :: fjac
if (nrhs < 2 .or. nlhs /= 2) then
call mexErrMsgTxt("Must have at least 2 inputs and exactly 2 outputs")
return
end if
if (.not. ((mxIsChar(prhs(1)) .and. mxGetM(prhs(1)) == 1) .or. mxIsClass(prhs(1), "function_handle"))) then
call mexErrMsgTxt("First argument should be a string or a function handle")
return
end if
if (.not. (mxIsDouble(prhs(2)) .and. (mxGetM(prhs(2)) == 1 .or. mxGetN(prhs(2)) == 1))) then
call mexErrMsgTxt("Second argument should be a real vector")
return
end if
func => prhs(1)
extra_args => prhs(3:nrhs)
associate (x_mat => mxGetPr(prhs(2)))
allocate(x(size(x_mat)))
x = x_mat
end associate
allocate(fvec(size(x)))
allocate(fjac(size(x), size(x)))
! Compute block decomposition
nullify(x_indices, f_indices, x_all)
call matlab_fcn(x, fvec, fjac)
call dm_blocks(fjac, blocks)
! Solve the system, starting from bottom-rightmost block
x_all => x
do i = size(blocks),1,-1
block
real(real64), dimension(size(blocks(i)%col_indices)) :: x_block
x_indices => blocks(i)%col_indices
f_indices => blocks(i)%row_indices
x_block = x(x_indices)
call trust_region_solve(x_block, matlab_fcn, info, tolf = tolf)
x(x_indices) = x_block
end block
end do
! Verify that we have a solution
! Note that here we use the ∞-norm, while trust region uses 2-norm; otherwise
! this check would almost always fail (because the 2-norm of the full fvec is
! larger than the 2-norm of its sub-vectors)
! If the check fails, this normally means that the block decomposition was
! incorrect (because some element of the Jacobian was numerically zero at the
! guess value while not being symbolically zero)
nullify(x_indices, f_indices, x_all)
call matlab_fcn(x, fvec)
if (maxval(abs(fvec)) > tolf) then
call trust_region_solve(x, matlab_fcn, info, tolf = tolf)
else
info = 1
end if
plhs(1) = mxCreateDoubleMatrix(int(size(x, 1), mwSize), 1_mwSize, mxREAL)
mxGetPr(plhs(1)) = x
if (info == 1) then
plhs(2) = mxCreateDoubleScalar(0._c_double)
else
plhs(2) = mxCreateDoubleScalar(1._c_double)
end if
end subroutine mexFunction
FC = gfortran
FCFLAGS = -g -Wall
all: dulmage_mendelsohn_test trust_region_test
dulmage_mendelsohn_test: dulmage_mendelsohn.o dulmage_mendelsohn_test.o
$(FC) $(FCFLAGS) $^ -o $@
dulmage_mendelsohn_test.o: dulmage_mendelsohn_test.f08 dulmage_mendelsohn.mod
$(FC) $(FCFLAGS) -c $< -o $@
dulmage_mendelsohn.mod: dulmage_mendelsohn.o
dulmage_mendelsohn.o: ../dulmage_mendelsohn.f08
$(FC) $(FCFLAGS) -c $< -o $@
trust_region_test: trust_region.o trust_region_test.o
$(FC) $(FCFLAGS) $^ -o $@ -llapack -lblas
trust_region_test.o: trust_region_test.f08 trust_region.mod
$(FC) $(FCFLAGS) -c $< -o $@
trust_region.mod: trust_region.o
trust_region.o: ../trust_region.f08 lapack.mod
$(FC) $(FCFLAGS) -c $< -o $@
blas_lapack.o: ../../blas_lapack.F08
$(FC) $(FCFLAGS) -c $< -o $@
lapack.mod: blas_lapack.o
clean:
rm -f *.o *.mod dulmage_mendelsohn_test trust_region_test
program dmperm_test
use iso_fortran_env
use dulmage_mendelsohn
implicit none
real(real64), dimension(12, 11) :: M
integer, dimension(size(M, 1)) :: row_order
integer, dimension(size(M, 2)) :: col_order
integer, dimension(:), allocatable :: row_blocks, col_blocks
type(dm_block), dimension(:), allocatable :: blocks
integer :: i, j
! Matrix given p. 305 of Pothen and Fan (1990)
M = reshape(real([ &
1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, &
0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, &
1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, &
0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, &
0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, &
0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, &
0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, &
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 &
], real64), shape(M), order = [2, 1])
! Shuffle the matrix
M = M([ 3, 9, 8, 7, 1, 2, 5, 4, 10, 12, 11, 6 ], &
[ 5, 9, 6, 1, 10, 8, 4, 11, 2, 7, 3 ])
print *, "Shuffled matrix:"
do i = 1, size(M, 1)
write (*, "(11f3.0) ") (M(i, j), j=1,size(M,2))
end do
print *
! Test dmperm
call dmperm(M, row_order, col_order, row_blocks, col_blocks)
print *, "Reordered matrix:"
do i = 1, size(M, 1)
write (*, "(11f3.0) ") (M(row_order(i), col_order(j)), j=1,size(M,2))
end do
print *
! Print fine decomposition
print *, "Column blocks:", col_blocks
print *, "Row blocks:", row_blocks
print *
! Test compute_dm_blocks
call dm_blocks(M, blocks)
do i = 1, size(blocks)
print *, "*** Block:", i
print *, "Coarse type: ", blocks(i)%coarse_type
print *, "Rows:", blocks(i)%row_indices
print *, "Columns:", blocks(i)%col_indices
print *, "Predecessors:", blocks(i)%predecessors
print *
end do
end program dmperm_test
program trust_region_test
use trust_region
use iso_fortran_env
implicit none
real(real64), dimension(2), parameter :: rosenbrock_guess = [ -10**4, 1 ]
real(real64), dimension(2), parameter :: rosenbrock_solution = [ 1, 1 ]
real(real64), dimension(2) :: rosenbrock_x
real(real64), dimension(4), parameter :: powell1_guess = [ 3, -1, 0, 1 ]
real(real64), dimension(4), parameter :: powell1_solution = 0
real(real64), dimension(4) :: powell1_x
integer :: info
rosenbrock_x = rosenbrock_guess
call trust_region_solve(rosenbrock_x, rosenbrock, info)
if (info /= 1 .or. maxval(abs(rosenbrock_x - rosenbrock_solution)) > 1e-11_real64) &
error stop "Failed to solve rosenbrock"
powell1_x = powell1_guess
call trust_region_solve(powell1_x, powell1, info, tolf = 1e-8_real64)
if (info /= 1 .or. maxval(abs(rosenbrock_x - rosenbrock_solution)) > 4e-5_real64) &
error stop "Failed to solve powell1"
contains
subroutine rosenbrock(x, fvec, fjac)
real(real64), dimension(:), intent(in) :: x
real(real64), dimension(size(x)), intent(out) :: fvec
real(real64), dimension(size(x),size(x)), intent(out), optional :: fjac
fvec(1) = 1 - x(1)
fvec(2) = 10*(x(2)-x(1)**2)
if (present(fjac)) then
fjac(1,1) = -1
fjac(1,2) = 0
fjac(2,1) = -20*x(1)
fjac(2,2) = 10
end if
end subroutine rosenbrock
subroutine powell1(x, fvec, fjac)
real(real64), dimension(:), intent(in) :: x
real(real64), dimension(size(x)), intent(out) :: fvec
real(real64), dimension(size(x),size(x)), intent(out), optional :: fjac
fvec(1) = x(1)+10*x(2)
fvec(2) = sqrt(5._real64)*(x(3)-x(4))
fvec(3) = (x(2)-2*x(3))**2
fvec(4) = sqrt(10._real64)*(x(1)-x(4))**2
if (present(fjac)) then
fjac(1,1) = 1
fjac(1,2) = 10
fjac(2,3) = sqrt(5._real64)
fjac(2,4) = -fjac(2,3)
fjac(3,2) = 2*(x(2)-2*x(3))
fjac(3,3) = -2*fjac(3,2)
fjac(4,1) = 2*sqrt(10._real64)*(x(1)-x(4))
fjac(4,4) = -fjac(4,1)
end if
end subroutine powell1
end program trust_region_test
! Solves a system of nonlinear equation with the trust region method
!
! Implementation heavily inspired from the hybrj function from MINPACK
! Copyright © 2019 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/>.
module trust_region
use iso_fortran_env
use lapack
implicit none
private
public :: trust_region_solve
contains
subroutine trust_region_solve(x, f, info, tolx, tolf, maxiter)
real(real64), dimension(:), intent(inout) :: x ! On entry: guess value; on exit: solution
interface
subroutine f(x1, fvec, fjac)
use iso_fortran_env
real(real64), dimension(:), intent(in) :: x1
real(real64), dimension(size(x)), intent(out) :: fvec
real(real64), dimension(size(x),size(x)), intent(out), optional :: fjac
end subroutine f
end interface
! Exit code:
! 1 = success (relative error between two consecutive iterates is at most tolx)
! 2 = maximum number of iterations reached
! 3 = tolx is too small, no further improvement of the approximate solution x is possible
! 4 = iteration is not making good progress, as measured by the improvement from the last maxslowiter iterations
integer, intent(out) :: info
real(real64), intent(in), optional :: tolx, tolf ! Tolerances in x and f
integer, intent(in), optional :: maxiter ! Maximum number of iterations
real(real64) :: tolx_actual, tolf_actual
integer :: maxiter_actual
integer, parameter :: maxslowiter = 15 ! Maximum number of consecutive iterations with slow progress
real(real64), parameter :: macheps = epsilon(x)
real(real64) :: delta ! Radius of the trust region
real(real64), dimension(size(x)) :: fvec ! Current function value
real(real64) :: fn ! Norm of the current function value
real(real64), dimension(size(x)) :: jcn, dg ! Jacobian column-norms and rescaling factors
real(real64), dimension(size(x), size(x)) :: fjac ! Jacobian matrix
real(real64), dimension(size(x)) :: gn ! Gauss-Newton direction
logical :: recompute_gn ! Whether to recompute Gauss-Newton direction at next dogleg
integer :: niter ! Current iteration
integer :: ncsucc ! Number of consecutive successful iterations
integer :: ncslow ! Number of consecutive iterations with slow progress
! Initialize variables associated to optional arguments
if (present(tolx)) then
tolx_actual = tolx
else
tolx_actual = 1e-6_real64
end if
if (present(tolf)) then
tolf_actual = tolf
else
tolf_actual = 1e-6_real64
end if
if (present(maxiter)) then
maxiter_actual = maxiter
else
maxiter_actual = 50
end if
niter = 1
ncsucc = 0