diff --git a/license.txt b/license.txt
index 61a2f6a5d3cc4c7230ccb888e2d6c22f06cb6b99..a479aeeb7154f00226c9604a57ee1a5d293d398c 100644
--- a/license.txt
+++ b/license.txt
@@ -118,11 +118,6 @@ Copyright: 2011 Lawrence J. Christiano, Mathias Trabandt and Karl Walentin
2013-2017 Dynare Team
License: GPL-3+
-Files: matlab/trust_region.m
-Copyright: 2008-2012 VZLU Prague, a.s.
- 2014-2021 Dynare Team
-License: GPL-3+
-
Files: matlab/one_sided_hp_filter.m
Copyright: 2010-2015 Alexander Meyer-Gohde
2015-2017 Dynare Team
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 0244855004ea6de900af46dc39facc4e3f1f4a11..f68d14cce3d26c6dfb7597c478417591c97180b5 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -14,6 +14,7 @@ function [x, errorflag, fvec, fjac, exitflag] = dynare_solve(f, x, options, vara
% - errorflag [logical] scalar, true iff the model can not be solved.
% - fvec [double] n×1 vector, function value at x (f(x), used for debugging when errorflag is true).
% - fjac [double] n×n matrix, Jacobian value at x (J(x), used for debugging when errorflag is true).
+% - exitflag [integer] scalar,
% Copyright © 2001-2022 Dynare Team
%
@@ -94,7 +95,7 @@ if jacobian_flag
if ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) ...
|| any(~isreal(fvec)) || any(~isreal(fjac(:)))
if max(abs(fvec)) < tolf %return if initial value solves problem
- info = 0;
+ exitflag = -1;
return;
end
disp_verbose('Randomize initial guess...',options.verbosity)
@@ -231,10 +232,7 @@ if options.solve_algo == 0
elseif options.solve_algo==1
[x, errorflag, exitflag] = solve1(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, [], options.debug, arguments{:});
elseif options.solve_algo==9
- [x, errorflag] = trust_region(f,x, 1:nn, 1:nn, jacobian_flag, options.gstep, ...
- tolf, tolx, maxit, ...
- options.trust_region_initial_step_bound_factor, ...
- options.debug, arguments{:});
+ [x, errorflag, exitflag] = trust_region(f, x, 1:nn, 1:nn, jacobian_flag, options.gstep, tolf, tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, arguments{:});
elseif ismember(options.solve_algo, [2, 12, 4])
if ismember(options.solve_algo, [2, 12])
solver = @solve1;
@@ -310,11 +308,11 @@ elseif ismember(options.solve_algo, [2, 12, 4])
dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u with trust_region routine.', i);
end
end
- [x, errorflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
- options.gstep, ...
- tolf, options.solve_tolx, maxit, ...
- options.trust_region_initial_step_bound_factor, ...
- options.debug, arguments{:});
+ [x, errorflag, exitflag] = solver(f, x, j1(j), j2(j), jacobian_flag, ...
+ options.gstep, ...
+ tolf, options.solve_tolx, maxit, ...
+ options.trust_region_initial_step_bound_factor, ...
+ options.debug, arguments{:});
fre = true;
if errorflag
return
@@ -323,10 +321,10 @@ elseif ismember(options.solve_algo, [2, 12, 4])
fvec = feval(f, x, arguments{:});
if max(abs(fvec))>tolf
disp_verbose('Call solver on the full nonlinear problem.',options.verbosity)
- [x, errorflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
- options.gstep, tolf, options.solve_tolx, maxit, ...
- options.trust_region_initial_step_bound_factor, ...
- options.debug, arguments{:});
+ [x, errorflag, exitflag] = solver(f, x, 1:nn, 1:nn, jacobian_flag, ...
+ options.gstep, tolf, options.solve_tolx, maxit, ...
+ options.trust_region_initial_step_bound_factor, ...
+ options.debug, arguments{:});
end
elseif options.solve_algo==3
if jacobian_flag
diff --git a/matlab/trust_region.m b/matlab/trust_region.m
index ab5d6a908d68fe87c3f82405c6452ece90a092fd..8fcaaaa664752688b4f19249b75e8d6de317d498 100644
--- a/matlab/trust_region.m
+++ b/matlab/trust_region.m
@@ -1,32 +1,42 @@
-function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,factor,debug,varargin)
+function [x, errorflag, info] = trust_region(objfun, x, j1, j2, jacobianflag, gstep, tolf, tolx, maxiter, factor, debug, varargin)
+
% Solves systems of non linear equations of several variables, using a
% trust-region method.
%
% INPUTS
-% fcn: name of the function to be solved
-% x0: guess values
-% j1: equations index for which the model is solved
-% j2: unknown variables index
-% jacobian_flag=true: jacobian given by the 'func' function
-% jacobian_flag=false: jacobian obtained numerically
-% gstep increment multiplier in numercial derivative
-% computation
-% tolf tolerance for residuals
-% tolx tolerance for solution variation
-% maxiter maximum number of iterations
-% factor real scalar, determines the initial step bound
-% debug debug flag
-% varargin: list of arguments following bad_cond_flag
+% - objfun [function handle, char] name of the routine evaluating the system of nonlinear equations (and possibly jacobian).
+% - x [double] n×1 vector, initial guess for the solution.
+% - j1 [integer] vector, equation indices defining a subproblem to be solved.
+% - j2 [integer] vector, unknown variable indices to be solved for in the subproblem.
+% - jacobianflag [logical] scalar, if true the jacobian matrix is expected to be returned as a second output argument when calling objfun, otherwise
+% the jacobian is computed numerically.
+% - gstep [double] scalar, increment multiplier in numerical derivative computation (only used if jacobianflag value is false).
+% - tolf [double] scalar, tolerance for residuals.
+% - tolx [double] scalar, tolerance for solution variation.
+% - maxiter [integer] scalar, maximum number of iterations;
+% - factor [double] scalar, determines the initial step bound.
+% - debug [logical] scalar, dummy argument.
+% - varargin: [cell] list of additional arguments to be passed to objfun.
%
% OUTPUTS
-% x: results
-% check=1: the model can not be solved
-% info: detailed exitcode
-% SPECIAL REQUIREMENTS
-% none
+% - x [double] n⨱1 vector, solution of the nonlinear system of equations.
+% - errorflag [logical] scalar, false iff nonlinear solver is successful.
+% - info [integer] scalar, information about the failure.
+%
+% REMARKS
+% [1] j1 and j2 muyst have the same number of elements.
+% [2] debug is here for compatibility purpose (see solve1), it does not affect the output.
+% [3] Possible values for info are:
+%
+% -1 if the initial guess is a solution of the nonlinear system of equations.
+% 0 if the nonlinear solver failed because the problem is ill behaved at the initial guess.
+% 1 if the nonlinear solver found a solution.
+% 2 if the maximum number of iterations has been reached.
+% 3 if spurious convergence (trust region radius is too small).
+% 4 if iteration is not making good progress, as measured by the improvement from the last 15 iterations.
+% 5 if no further improvement in the approximate solution x is possible (xtol is too small).
-% Copyright (C) 2008-2012 VZLU Prague, a.s.
-% Copyright (C) 2014-2021 Dynare Team
+% Copyright © 2014-2022 Dynare Team
%
% This file is part of Dynare.
%
@@ -42,186 +52,275 @@ function [x,check,info] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tol
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-%
-% Initial author: Jaroslav Hajek <highegg@gmail.com>, for GNU Octave
-if (ischar (fcn))
- fcn = str2func (fcn);
+% Convert to function handle if necessary
+if ischar(objfun)
+ objfun = str2func(objfun);
end
-n = length(j1);
+%
+% Set constants
+%
+
+radiusfactor = .5;
+actualreductionthreshold = .001;
+relativereductionthreshold = .1;
-% These defaults are rather stringent. I think that normally, user
-% prefers accuracy to performance.
+% number of equations
+n = length(j1);
-macheps = eps (class (x0));
+%
+% Initialization
+%
-niter = 1;
+errorflag = true;
-x = x0;
+iter = 1;
info = 0;
+delta = 0.0;
-% Initial evaluation.
-% Handle arbitrary shapes of x and f and remember them.
-fvec = fcn (x, varargin{:});
-fvec = fvec(j1);
-fn = norm (fvec);
-recompute_jacobian = true;
+ncsucc = 0; ncslow = 0;
-% Outer loop.
-while (niter < maxiter && ~info)
+%
+% Attempt to evaluate the residuals and jacobian matrix on the initial guess
+%
- % Calculate Jacobian (possibly via FD).
- if recompute_jacobian
- if jacobian_flag
- [~, fjac] = fcn (x, varargin{:});
- fjac = fjac(j1,j2);
- else
- dh = max(abs(x(j2)),gstep(1)*ones(n,1))*eps^(1/3);
-
- for j = 1:n
- xdh = x ;
- xdh(j2(j)) = xdh(j2(j))+dh(j) ;
- t = fcn(xdh,varargin{:});
- fjac(:,j) = (t(j1) - fvec)./dh(j) ;
- end
+try
+ if jacobianflag
+ [fval, fjac] = objfun(x, varargin{:});
+ fval = fval(j1);
+ fjac = fjac(j1,j2);
+ else
+ fval = objfun(x, varargin{:});
+ fval = fval(j1);
+ dh = max(abs(x(j2)), gstep(1)*ones(n,1))*eps^(1/3);
+ fjac = zeros(n);
+ for j = 1:n
+ xdh = x;
+ xdh(j2(j)) = xdh(j2(j))+dh(j);
+ Fval = objfun(xdh, varargin{:});
+ fjac(:,j) = (Fval(j1)-fval)./dh(j);
end
- recompute_jacobian = false;
end
+ fnorm = norm(fval);
+catch
+ % System of equation cannot be evaluated at the initial guess.
+ return
+end
- % Get column norms, use them as scaling factors.
- jcn = sqrt(sum(fjac.*fjac))';
- if (niter == 1)
- dg = jcn;
- dg(dg == 0) = 1;
- else
- % Rescale adaptively.
- dg = max (dg, jcn);
- end
+if any(isnan(fval)) || any(isinf(fval)) || any(~isreal(fval)) || any(isnan(fjac(:))) || any(isinf(fjac(:))) || any(~isreal(fjac(:)))
+ % System of equations is ill-behaved at the initial guess.
+ return
+end
- if (niter == 1)
- xn = norm (dg .* x(j2));
- % FIXME: something better?
- delta = max (xn, 1)*factor;
- end
+if fnorm<tolf
+ % Initial guess is a solution.
+ errorflag = false;
+ info = -1;
+ return
+end
- % Get trust-region model (dogleg) minimizer.
- s = - dogleg (fjac, fvec, dg, delta);
- w = fvec + fjac * s;
+%
+% Main loop
+%
- sn = norm (dg .* s);
- if (niter == 1)
- delta = min (delta, sn);
+while iter<=maxiter && ~info
+ % Compute the columns norm ofr the Jacobian matrix.
+ fjacnorm = transpose(sqrt(sum((fjac.*fjac))));
+ if iter==1
+ % On the first iteration, calculate the norm of the scaled vector of unknowns x
+ % and initialize the step bound delta. Scaling is done according to the norms of
+ % the columns of the initial jacobian.
+ fjacnorm__ = fjacnorm;
+ fjacnorm__(fjacnorm<eps(1.0)) = 1.0;
+ xnorm = norm(fjacnorm__.*x(j2));
+ if xnorm>0
+ delta = xnorm*factor;
+ else
+ delta = factor;
+ end
+ else
+ fjacnorm__ = max(.1*fjacnorm__, fjacnorm);
+ xnorm = norm(fjacnorm__.*x(j2));
end
-
- x2 = x;
- x2(j2) = x2(j2) + s;
- fvec1 = fcn (x2, varargin{:});
- fvec1 = fvec1(j1);
- fn1 = norm (fvec1);
-
- if (fn1 < fn)
- % Scaled actual reduction.
- actred = 1 - (fn1/fn)^2;
+ % Determine the direction p (with trust region model defined in dogleg routine).
+ p = dogleg(fjac, fval, fjacnorm__, delta);
+ % Compute the norm of p.
+ pnorm = norm(fjacnorm__.*p);
+ xx = x;
+ x0 = x;
+ % Move along the direction p. Set a candidate value for x and predicted improvement for f.
+ xx(j2) = x(j2) - p;
+ ww = fval - fjac*p;
+ % Evaluate the function at xx and calculate its norm.
+ try
+ fval1 = objfun(xx, varargin{:});
+ fval1 = fval1(j1);
+ catch
+ % If evaluation of the residuals returns an error, then restart but with a smaller radius of the trust region.
+ delta = delta*radiusfactor;
+ iter = iter+1;
+ continue
+ end
+ if any(isnan(fval1)) || any(isinf(fval1)) || any(~isreal(fval1))
+ % If evaluation of the residuals returns a NaN, an infinite number or a complex number, then restart but with a smaller radius of the trust region.
+ delta = delta*radiusfactor;
+ iter = iter+1;
+ continue
+ end
+ fnorm1 = norm(fval1);
+ if fnorm1<tolf
+ x = xx;
+ errorflag = false;
+ info = 1;
+ continue
+ end
+ % Compute the scaled actual reduction.
+ if fnorm1<fnorm
+ actualreduction = 1.0-(fnorm1/fnorm)^2;
else
- actred = -1;
+ actualreduction = -1.0;
end
-
- % Scaled predicted reduction, and ratio.
- t = norm (w);
- if (t < fn)
- prered = 1 - (t/fn)^2;
- ratio = actred / prered;
+ % Compute the scaled predicted reduction and the ratio of the actual to the
+ % predicted reduction.
+ tt = norm(ww);
+ if tt<fnorm
+ predictedreduction = 1.0 - (tt/fnorm)^2;
+ else
+ predictedreduction = 0.0;
+ end
+ if predictedreduction>0
+ ratio = actualreduction/predictedreduction;
else
- prered = 0;
ratio = 0;
end
-
- % Update delta.
- if (ratio < 0.1)
- delta = 0.5*delta;
- if (delta <= 1e1*macheps*xn)
- % Trust region became uselessly small.
- if (fn1 <= tolf)
- info = 1;
- else
- info = -3;
- end
- break
+ % Update the radius of the trust region if need be.
+ if iter==1
+ % On first iteration adjust the initial step bound.
+ delta = min(delta, pnorm);
+ end
+ if ratio<relativereductionthreshold
+ % Reduction is much smaller than predicted… Reduce the radius of the trust region.
+ ncsucc = 0;
+ delta = delta*radiusfactor;
+ else
+ ncsucc = ncsucc + 1;
+ if abs(ratio-1.0)<relativereductionthreshold
+ delta = pnorm/radiusfactor;
+ elseif ratio>=radiusfactor || ncsucc>1
+ delta = max(delta, pnorm/radiusfactor);
end
- elseif (abs (1-ratio) <= 0.1)
- delta = 1.4142*sn;
- elseif (ratio >= 0.5)
- delta = max (delta, 1.4142*sn);
end
-
- if (ratio >= 1e-4)
- % Successful iteration.
- x(j2) = x(j2) + s;
- xn = norm (dg .* x(j2));
- fvec = fvec1;
- fn = fn1;
- recompute_jacobian = true;
+ if ratio>1e-4
+ % Successful iteration. Update x, xnorm, fval, fnorm and fjac.
+ x = xx;
+ fval = fval1;
+ xnorm = norm(fjacnorm__.*x(j2));
+ fnorm = fnorm1;
end
-
- niter = niter + 1;
-
-
- % Tests for termination condition
- if (fn <= tolf)
- info = 1;
+ % Determine the progress of the iteration.
+ ncslow = ncslow+1;
+ if actualreduction>=actualreductionthreshold
+ ncslow = 0;
+ end
+ iter = iter+1;
+ if iter==maxiter
+ info = 2;
+ x(:) = inf;
+ continue
+ end
+ if delta<tolx*xnorm
+ info = 3;
+ x(:) = inf;
+ errorflag = true;
+ continue
+ end
+ % Tests for termination and stringent tolerances.
+ if max(.1*delta, pnorm)<=10*eps(xnorm)*xnorm
+ % xtol is too small. no further improvement in
+ % the approximate solution x is possible.
+ info = 5;
+ x(:) = inf;
+ errorflag = true;
+ continue
+ end
+ if ncslow==15
+ info = 4;
+ x(:) = inf;
+ errorflag = true;
+ continue
+ end
+ % Compute the jacobian for the next iteration.
+ if jacobianflag
+ try
+ [~, fjac] = objfun(x, varargin{:});
+ fjac = fjac(j1,j2);
+ catch
+ % If evaluation of the Jacobian matrix returns an error, then restart but with a smaller radius of the trust region.
+ x = x0;
+ delta = delta*radiusfactor;
+ continue
+ end
+ else
+ dh = max(abs(x(j2)), gstep(1)*ones(n,1))*eps^(1/3);
+ for j = 1:n
+ xdh = x;
+ xdh(j2(j)) = xdh(j2(j))+dh(j);
+ Fval = objfun(xdh, varargin{:});
+ fjac(:,j) = (Fval(j1)-fval)./dh(j);
+ end
+ end
+ if any(isnan(fjac(:))) || any(isinf(fjac(:))) || any(~isreal(fjac(:)))
+ % If evaluation of the Jacobian matrix returns NaNs, an infinite numbers or a complex numbers, then restart but with a smaller radius of the trust region.
+ x = x0;
+ delta = delta*radiusfactor;
end
end
-if info==1
- check = 0;
-else
- check = 1;
-end
-end
-
-% Solve the double dogleg trust-region least-squares problem:
-% Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
-% x being a convex combination of the gauss-newton and scaled gradient.
-
-% TODO: error checks
-% TODO: handle singularity, or leave it up to mldivide?
function x = dogleg (r, b, d, delta)
-% Get Gauss-Newton direction.
+% Compute the Gauss-Newton direction.
if isoctave || matlab_ver_less_than('9.3')
- % The decomposition() function does not exist in Octave and MATLAB < R2017b
+ % The decomposition() function does not exist in Octave and MATLAB < R2017b
x = r \ b;
else
x = decomposition(r, 'CheckCondition', false) \ b;
end
-xn = norm (d .* x);
-if (xn > delta)
- % GN is too big, get scaled gradient.
- s = (r' * b) ./ d;
- sn = norm (s);
- if (sn > 0)
- % Normalize and rescale.
- s = (s / sn) ./ d;
+% Compute norm of scaled x
+qnorm = norm(d.*x);
+if qnorm<=delta
+ % Gauss-Newton direction is acceptable. There is nothing to do here.
+else
+ % Gauss-Newton direction is not acceptable…
+ % Compute the scale gradient direction and its norm
+ s = (r'*b)./d;
+ gnorm = norm(s);
+ if gnorm>0
+ % Normalize and rescale → gradient direction.
+ s = (s/gnorm)./d;
% Get the line minimizer in s direction.
- tn = norm (r*s);
- snm = (sn / tn) / tn;
- if (snm < delta)
- % Get the dogleg path minimizer.
- bn = norm (b);
- dxn = delta/xn; snmd = snm/delta;
- t = (bn/sn) * (bn/xn) * snmd;
- t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
- alpha = dxn*(1-snmd^2) / t;
+ temp0 = norm(r*s);
+ sgnorm = gnorm/(temp0*temp0);
+ if sgnorm<delta
+ % The scaled gradient direction is not acceptable…
+ % Compute the point along the dogleg at which the
+ % quadratic is minimized.
+ bnorm = norm(b);
+ temp1 = delta/qnorm;
+ temp2 = sgnorm/delta;
+ temp0 = bnorm*bnorm*temp2/(gnorm*qnorm);
+ temp0 = temp0 - temp1*temp2*temp2 + sqrt((temp0-temp1)^2+(1.0-temp1*temp1)*(1.0-temp2*temp2));
+ alpha = temp1*(1.0-temp2*temp2)/temp0;
else
- alpha = 0;
+ % The scaled gradient direction is acceptable.
+ alpha = 0.0;
end
else
- alpha = delta / xn;
- snm = 0;
+ % If the norm of the scaled gradient direction is zero.
+ alpha = delta/qnorm;
+ sgnorm = 0.0;
end
- % Form the appropriate convex combination.
- x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
-end
-end
-
+ % Form the appropriate convex combination of the Gauss-Newton direction and the
+ % scaled gradient direction.
+ x = alpha*x + (1.0-alpha)*min(sgnorm, delta)*s;
+end
\ No newline at end of file
diff --git a/tests/nonlinearsolvers.m b/tests/nonlinearsolvers.m
index 7d820076454928031980264e5db5319753ca0943..5231dcd7b02465d8cbd8c2a29ea9b8f5898a9d90 100644
--- a/tests/nonlinearsolvers.m
+++ b/tests/nonlinearsolvers.m
@@ -33,176 +33,120 @@ factor = 10;
auxstruct = struct();
-t0 = clock;
-
-try
- NumberOfTests = NumberOfTests+1;
- x = rosenbrock();
- [x, errorflag] = block_trust_region(@rosenbrock, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
+% List of function handles
+objfun = { @rosenbrock,
+ @powell1,
+ @powell2,
+ @wood,
+ @helicalvalley,
+ @watson,
+ @chebyquad,
+ @brown,
+ @discreteboundaryvalue,
+ @discreteintegralequation,
+ @trigonometric,
+ @variablydimensioned,
+ @broydentridiagonal,
+ @broydenbanded };
+
+% FIXME block_trust_region (mex) or trust_region (matlab) do not work for all n (not sure we can fix that).
+% FIXME block_trust_region (mex) and trust_region (matlab) do not behave the same (spurious convergence for powell2 and trigonometric with block_trust_region).
+
+%
+% Test mex routine
+%
-try
- NumberOfTests = NumberOfTests+1;
- x = powell1();
- [x, errorflag] = block_trust_region(@powell1, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- x = powell2();
- [x, errorflag] = block_trust_region(@powell2, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- x = wood();
- [x, errorflag] = block_trust_region(@wood, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- % FIXME block_trust_region is diverging if x(1)<0. Note that trust_region is not finding the
- % solution for the same initial conditions.
- x = helicalvalley();
- x(1) = 5;
- [x, errorflag] = block_trust_region(@helicalvalley, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- n = 10;
- x = watson(nan(n,1));
- [x, errorflag] = block_trust_region(@watson, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- % FIXME block_trust_region does not work for all n.
- n = 9;
- x = chebyquad(nan(n,1));
- [x, errorflag] = block_trust_region(@chebyquad, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- n = 10;
- x = brown(nan(n,1));
- [x, errorflag] = block_trust_region(@brown, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
+t0 = clock;
-try
+for i=1:length(objfun)
NumberOfTests = NumberOfTests+1;
- n = 10;
- x = discreteboundaryvalue(nan(n,1));
- [x, errorflag] = block_trust_region(@discreteboundaryvalue, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
+ switch func2str(objfun{i})
+ case 'helicalvalley'
+ % FIXME block_trust_region is diverging if x(1)<0.
+ x = helicalvalley();
+ x(1) = 5;
+ case 'chebyquad'
+ % Fails with a system of 10 equations.
+ x = objfun{i}(nan(9,1));
+ case {'watson', 'brown', 'discreteintegralequation', 'discreteboundaryvalue', 'chebyquad', 'trigonometric', 'variablydimensioned', 'broydenbanded', 'broydentridiagonal'}
+ x = objfun{i}(nan(4,1));
+ otherwise
+ x = objfun{i}();
end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- n = 10;
- x = discreteintegralequation(nan(n,1));
- [x, errorflag] = block_trust_region(@discreteintegralequation, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
+ try
+ [x, errorflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, false, auxstruct);
+ if isequal(func2str(objfun{i}), 'powell2')
+ if ~errorflag
+ testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver is expected to fail on %s function but did not return an error.', func2str(objfun{i}))
+ end
+ end
+ else
+ if errorflag || norm(objfun{i}(x))>tolf
+ testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver (mex) failed on %s function (norm(f(x))=%s).', func2str(objfun{i}), num2str(norm(objfun{i}(x))))
+ end
+ end
+ end
+ catch
testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver (mex) failed on %s function.', func2str(objfun{i}))
+ end
end
-catch
- testFailed = testFailed+1;
end
-try
- NumberOfTests = NumberOfTests+1;
- n = 10;
- x = trigonometric(nan(n,1));
- [x, errorflag] = block_trust_region(@trigonometric, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
+t1 = clock; etime(t1, t0)
-try
- NumberOfTests = NumberOfTests+1;
- n = 10;
- x = variablydimensioned(nan(n,1));
- [x, errorflag] = block_trust_region(@variablydimensioned, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
- end
-catch
- testFailed = testFailed+1;
-end
+%
+% Test matlab routine
+%
-try
+for i=1:length(objfun)
NumberOfTests = NumberOfTests+1;
- n = 10;
- x = broydentridiagonal(nan(n,1));
- [x, errorflag] = block_trust_region(@broydentridiagonal, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
- testFailed = testFailed+1;
+ switch func2str(objfun{i})
+ case 'chebyquad'
+ % Fails with a system of 10 equations.
+ x = objfun{i}(nan(9,1));
+ case {'watson', 'brown', 'discreteintegralequation', 'discreteboundaryvalue', 'trigonometric', 'variablydimensioned', 'broydenbanded', 'broydentridiagonal'}
+ x = objfun{i}(nan(10,1));
+ otherwise
+ x = objfun{i}();
end
-catch
- testFailed = testFailed+1;
-end
-
-try
- NumberOfTests = NumberOfTests+1;
- n = 10;
- x = broydenbanded(nan(n,1));
- [x, errorflag] = block_trust_region(@broydenbanded, x, tolf, tolx, maxit, factor, false, auxstruct);
- if errorflag
+ try
+ [x, errorflag, info] = trust_region(objfun{i}, x, 1:length(x), 1:length(x), true, [], tolf, tolx, maxit, factor);
+ if isequal(func2str(objfun{i}), 'powell2')
+ if ~errorflag
+ testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver is expected to fail on %s function but did not return an error.', func2str(objfun{i}))
+ end
+ end
+ if info~=3
+ testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver is expected to fail on %s function with info==3 but did not the correct value of info.', func2str(objfun{i}))
+ end
+ end
+ else
+ if errorflag
+ testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver failed on %s function (info=%s).', func2str(objfun{i}), int2str(info))
+ end
+ end
+ end
+ catch
testFailed = testFailed+1;
+ if debug
+ dprintf('Nonlinear solver failed on %s function.', func2str(objfun{i}))
+ end
end
-catch
- testFailed = testFailed+1;
end
-t1 = clock;
+t2 = clock; etime(t2, t1)
if ~debug
cd(getenv('TOP_TEST_DIR'));
@@ -216,14 +160,14 @@ else
fid = fopen('nonlinearsolvers.m.trs', 'w+');
end
if testFailed
- fprintf(fid,':test-result: FAIL\n');
+ fprintf(fid,':test-result: FAIL\n');
else
- fprintf(fid,':test-result: PASS\n');
+ fprintf(fid,':test-result: PASS\n');
end
fprintf(fid,':number-tests: %i\n', NumberOfTests);
fprintf(fid,':number-failed-tests: %i\n', testFailed);
fprintf(fid,':list-of-passed-tests: nonlinearsolvers.m\n');
-fprintf(fid,':elapsed-time: %f\n', etime(t1, t0));
+fprintf(fid,':elapsed-time: %f\n', etime(t2, t0));
fclose(fid);
if ~debug