From 3025a14ed9074c684a96b428be939b60a3c2d15a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Mon, 25 Mar 2013 12:03:09 +0100 Subject: [PATCH] Adapt for removal of luinc in MATLAB R2013a --- matlab/dynare_config.m | 5 ++++ matlab/missing/ilu/ilu.m | 43 ++++++++++++++++++++++++++++ matlab/solve_one_boundary.m | 12 ++++---- matlab/solve_two_boundaries.m | 10 +++---- mex/sources/bytecode/SparseMatrix.cc | 10 +++++-- 5 files changed, 66 insertions(+), 14 deletions(-) create mode 100644 matlab/missing/ilu/ilu.m diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index 7bb701b56..6d77c52fb 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -89,6 +89,11 @@ if ~exist('OCTAVE_VERSION') && matlab_ver_less_than('7.4') addpath([dynareroot '/missing/bsxfun']) end +% ilu is missing in old versions of MATLAB and in Octave +if exist('OCTAVE_VERSION') || matlab_ver_less_than('7.4') + addpath([dynareroot '/missing/ilu']) +end + % nanmean is in Octave Forge Statistics package and in MATLAB Statistics % toolbox if (exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('statistics')) ... diff --git a/matlab/missing/ilu/ilu.m b/matlab/missing/ilu/ilu.m new file mode 100644 index 000000000..1e4dfcf1e --- /dev/null +++ b/matlab/missing/ilu/ilu.m @@ -0,0 +1,43 @@ +function [L, U, P] = ilu(A, setup) +% Partially implement function ilu using the (now deprecated) luinc + +% Copyright (C) 2013 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 . + + if nargout ~= 2 + error('Only two output arguments supported') + end + if nargin ~= 2 + error('Only two input arguments supported') + end + + if isfield(setup, 'type') + if setup.type ~= 'ilutp' + error(['Unsupported type: ' setup.type ]) + end + 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 diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 8ae660889..3e8825575 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -55,7 +55,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ... % none. % -% Copyright (C) 1996-2012 Dynare Team +% Copyright (C) 1996-2013 Dynare Team % % This file is part of Dynare. % @@ -77,7 +77,7 @@ Blck_size=size(y_index_eq,2); g2 = []; g3 = []; correcting_factor=0.01; -luinc_tol=1e-10; +ilu_setup.droptol=1e-10; max_resa=1e100; reduced = 0; if(forward_backward) @@ -316,7 +316,7 @@ for it_=start:incr:finish disp('steady: GMRES '); end while(flag1>0) - [L1, U1]=luinc(g1,luinc_tol); + [L1, U1]=ilu(g1,ilu_setup); [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1); if (flag1>0 || reduced) if(flag1==1) @@ -326,7 +326,7 @@ for it_=start:incr:finish elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]); end; - luinc_tol = luinc_tol/10; + ilu_setup.droptol = ilu_setup.droptol/10; reduced = 0; else ya = ya + lambda*dx; @@ -343,7 +343,7 @@ for it_=start:incr:finish disp('steady: BiCGStab'); end while(flag1>0) - [L1, U1]=luinc(g1,luinc_tol); + [L1, U1]=ilu(g1,ilu_setup); phat = ya - U1 \ (L1 \ r); if(is_dynamic) y(it_,y_index_eq) = phat; @@ -370,7 +370,7 @@ for it_=start:incr:finish elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]); end; - luinc_tol = luinc_tol/10; + ilu_setup.droptol = ilu_setup.droptol/10; reduced = 0; else ya = ya + lambda*dx; diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 27db7f4a5..942d0aa80 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -43,7 +43,7 @@ function y = solve_two_boundaries(fname, y, x, params, steady_state, y_index, nz % none. % -% Copyright (C) 1996-2011 Dynare Team +% Copyright (C) 1996-2013 Dynare Team % % This file is part of Dynare. % @@ -67,7 +67,7 @@ g2 = []; g3 = []; Blck_size=size(y_index,2); correcting_factor=0.01; -luinc_tol=1e-10; +ilu_setup.droptol=1e-10; max_resa=1e100; Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods); g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); @@ -198,7 +198,7 @@ while ~(cvg==1 || iter>maxit_), flag1=1; while(flag1>0) if preconditioner == 2 - [L1, U1]=luinc(g1a,luinc_tol); + [L1, U1]=ilu(g1a,ilu_setup); elseif preconditioner == 3 Size = Blck_size; gss1 = g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size); @@ -233,7 +233,7 @@ while ~(cvg==1 || iter>maxit_), elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); end; - luinc_tol = luinc_tol/10; + ilu_setup.droptol = ilu_setup.droptol/10; reduced = 0; else dx = za - ya; @@ -271,7 +271,7 @@ while ~(cvg==1 || iter>maxit_), elseif(flag1==3) disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); end; - luinc_tol = luinc_tol/10; + ilu_setup.droptol = ilu_setup.droptol/10; reduced = 0; else dx = za - ya; diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index baf1d9e89..6fb6f02b7 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2007-2012 Dynare Team + * Copyright (C) 2007-2013 Dynare Team * * This file is part of Dynare. * @@ -4589,11 +4589,15 @@ dynSparseMatrix::Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double throw FatalExceptionHandling(tmp.str()); #endif size_t n = mxGetM(A_m); + const char *field_names[] = {"droptol"}; + mwSize dims[1] = { 1 }; + mxArray *Setup = mxCreateStructArray(1, dims, 1, field_names); + mxSetFieldByNumber(Setup, 0, 0, mxCreateDoubleScalar(lu_inc_tol)); mxArray *lhs0[2]; mxArray *rhs0[2]; rhs0[0] = A_m; - rhs0[1] = mxCreateDoubleScalar(lu_inc_tol); - mexCallMATLAB(2, lhs0, 2, rhs0, "luinc"); + rhs0[1] = Setup; + mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"); mxArray *L1 = lhs0[0]; mxArray *U1 = lhs0[1]; /*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/