diff --git a/license.txt b/license.txt
index 10f10415c..94e618186 100644
--- a/license.txt
+++ b/license.txt
@@ -88,6 +88,25 @@ License: GPL-3+
You should have received a copy of the GNU General Public License
along with Dynare. If not, see .
+Files: matlab/missing/bicgstab/bicgstab.m
+Copyright: 2008, Radek Salac
+License: GPL-3+
+ This file is part of Octave.
+ .
+ Octave 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.
+ .
+ Octave 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 Octave; see the file COPYING. If not, see
+ .
+
Files: doc/manual.xml, doc/*.tex, doc/*.svg, doc/*.dia
Copyright: 1996-2010, Dynare Team
License: GFDL-1.3+
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 0027f7db0..8fc38b9f6 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -61,9 +61,10 @@ if exist('OCTAVE_VERSION') || matlab_ver_less_than('7.0.1')
addpath([dynareroot '/missing/ordeig'])
end
-% rcond() was introduced in Octave 3.2.0
+% rcond() and bicgstable() were introduced in Octave 3.2.0
if exist('OCTAVE_VERSION') && octave_ver_less_than('3.2.0')
addpath([dynareroot '/missing/rcond'])
+ addpath([dynareroot '/missing/bicgstab'])
end
% orschur() is missing in Octave; we don't have a real replacement;
diff --git a/matlab/missing/bicgstab/bicgstab.m b/matlab/missing/bicgstab/bicgstab.m
new file mode 100644
index 000000000..07083a1e4
--- /dev/null
+++ b/matlab/missing/bicgstab/bicgstab.m
@@ -0,0 +1,172 @@
+## Copyright (C) 2008 Radek Salac
+##
+## This file is part of Octave.
+##
+## Octave 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.
+##
+## Octave 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 Octave; see the file COPYING. If not, see
+## .
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} bicgstab (@var{A}, @var{b})
+## @deftypefnx {Function File} {} bicgstab (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## This procedure attempts to solve a system of linear equations A*x = b for x.
+## The @var{A} must be square, symmetric and positive definite real matrix N*N.
+## The @var{b} must be a one column vector with a length of N.
+## The @var{tol} specifies the tolerance of the method, the default value is 1e-6.
+## The @var{maxit} specifies the maximum number of iterations, the default value is min(20,N).
+## The @var{M1} specifies a preconditioner, can also be a function handler which returns M\X.
+## The @var{M2} combined with @var{M1} defines preconditioner as preconditioner=M1*M2.
+## The @var{x0} is the initial guess, the default value is zeros(N,1).
+##
+## The value @var{x} is a computed result of this procedure.
+## The value @var{flag} can be 0 when we reach tolerance in @var{maxit} iterations, 1 when
+## we don't reach tolerance in @var{maxit} iterations and 3 when the procedure stagnates.
+## The value @var{relres} is a relative residual - norm(b-A*x)/norm(b).
+## The value @var{iter} is an iteration number in which x was computed.
+## The value @var{resvec} is a vector of @var{relres} for each iteration.
+##
+## @end deftypefn
+
+function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, M1, M2, x0)
+
+ if (nargin < 2 || nargin > 7 || nargout > 5)
+ print_usage ();
+ elseif (!isnumeric (A) || rows (A) != columns (A))
+ error ("bicgstab: the first argument must be a n-by-n matrix");
+ elseif (!isvector (b))
+ error ("bicgstab: b must be a vector");
+ elseif (!any (b))
+ error ("bicgstab: b shuldn't be a vector of zeros");
+ elseif (rows (A) != rows (b))
+ error ("bicgstab: the first and second argument must have the same number of rows");
+ elseif (nargin > 2 && !isscalar (tol))
+ error ("bicgstab: tol must be a scalar");
+ elseif (nargin > 3 && !isscalar (maxit))
+ error ("bicgstab: maxit must be a scalar");
+ elseif (nargin > 4 && ismatrix (M1) && (rows (M1) != rows (A) || columns (M1) != columns (A)))
+ error ("bicgstab: M1 must have the same number of rows and columns as A");
+ elseif (nargin > 5 && (!ismatrix (M2) || rows (M2) != rows (A) || columns (M2) != columns (A)))
+ error ("bicgstab: M2 must have the same number of rows and columns as A");
+ elseif (nargin > 6 && !isvector (x0))
+ error ("bicgstab: x0 must be a vector");
+ elseif (nargin > 6 && rows (x0) != rows (b))
+ error ("bicgstab: x0 must have the same number of rows as b");
+ endif
+
+ ## Default tolerance.
+ if (nargin < 3)
+ tol = 1e-6;
+ endif
+
+ ## Default maximum number of iteration.
+ if (nargin < 4)
+ maxit = min (rows (b), 20);
+ endif
+
+ ## Left preconditioner.
+ if (nargin == 5)
+ if (isnumeric (M1))
+ precon = @(x) M1 \ x;
+ endif
+ elseif (nargin > 5)
+ if (issparse (M1) && issparse (M2))
+ precon = @(x) M2 \ (M1 \ x);
+ else
+ M = M1*M2;
+ precon = @(x) M \ x;
+ endif
+ else
+ precon = @(x) x;
+ endif
+
+ ## specifies initial estimate x0
+ if (nargin < 7)
+ x = zeros (rows (b), 1);
+ else
+ x = x0;
+ endif
+
+ norm_b = norm (b);
+
+ res = b - A*x;
+ rr = res;
+
+ ## Vector of the residual norms for each iteration.
+ resvec = [norm(res)/norm_b];
+
+ ## Default behaviour we don't reach tolerance tol within maxit iterations.
+ flag = 1;
+
+ for iter = 1:maxit
+ rho_1 = res' * rr;
+
+ if (iter == 1)
+ p = res;
+ else
+ beta = (rho_1 / rho_2) * (alpha / omega);
+ p = res + beta * (p - omega * v);
+ endif
+
+ phat = precon (p);
+
+ v = A * phat;
+ alpha = rho_1 / (rr' * v);
+ s = res - alpha * v;
+
+ shat = precon (s);
+
+ t = A * shat;
+ omega = (t' * s) / (t' * t);
+ x = x + alpha * phat + omega * shat;
+ res = s - omega * t;
+ rho_2 = rho_1;
+
+ relres = norm (res) / norm_b;
+ resvec = [resvec; relres];
+
+ if (relres <= tol)
+ ## We reach tolerance tol within maxit iterations.
+ flag = 0;
+ break;
+ elseif (resvec (end) == resvec (end - 1))
+ ## The method stagnates.
+ flag = 3;
+ break;
+ endif
+ endfor
+
+ if (nargout < 2)
+ if (flag == 0)
+ printf (["bicgstab converged at iteration %i ",
+ "to a solution with relative residual %e\n"],iter,relres);
+ elseif (flag == 3)
+ printf (["bicgstab stopped at iteration %i ",
+ "without converging to the desired tolerance %e\n",
+ "because the method stagnated.\n",
+ "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres);
+ else
+ printf (["bicgstab stopped at iteration %i ",
+ "without converging to the desired toleranc %e\n",
+ "because the maximum number of iterations was reached.\n",
+ "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres);
+ endif
+ endif
+
+endfunction
+
+%!demo
+%! % Solve system of A*x=b
+%! A = [5 -1 3;-1 2 -2;3 -2 3]
+%! b = [7;-1;4]
+%! [x, flag, relres, iter, resvec] = bicgstab(A, b)
+