diff --git a/scripts/sparse/bicg.m b/scripts/sparse/bicg.m
--- a/scripts/sparse/bicg.m
+++ b/scripts/sparse/bicg.m
@@ -1,6 +1,5 @@
## Copyright (C) 2006 Sylvain Pelissier
-## Copyright (C) 2012-2015 Carlo de Falco
-## Copyright (C) 2016 Cristiano Dorigo, Octave Arena
+## Copyright (C) 2012-2016 Carlo de Falco
##
## This file is part of Octave.
##
@@ -19,61 +18,35 @@
## .
## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
-## @deftypefnx {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
-## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, @dots{})
+## @deftypefn {} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## @deftypefnx {} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
+## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, @dots{})
## Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
##
-## The input parameters are:
-##
## @itemize @minus
-##
-## @item @var{A} it is a square matrix. @var{A} can be passed as a matrix or
-## as a function handle or inline function
-## @code{f} such that @code{f(x, "notransp") = A*x} and
-## @code{f(x, "transp") = A'*x}.
-##
-## @item @var{b} is the right hand side vector. It must be a column vector
-## with same number of rows of @var{A}.
-##
## @item @var{rtol} is the relative tolerance, if not given or set to [] the
## default value 1e-6 is used.
##
## @item @var{maxit} the maximum number of outer iterations, if not given or
## set to [] the default value @code{min (20, numel (b))} is used.
##
-## @item @var{M1}, @var{M2} are the preconditioners. The
-## preconditioner @var{P} is given as @code{@var{P} = @var{M1} * @var{M2}}.
-## Both @var{M1} and @var{M2} can be passed as a matrix or as a
-## function handle or inline
-## function @code{g} such that @code{g(@var{x}, "notransp") =
-## @var{M1} \ @var{x}} or
-## @code{g(@var{x}, "notransp") = @var{M2} \ @var{x}}
-## and @code{g(@var{x}, "transp") = @var{M1}' \ @var{x}} or
-## @code{g(@var{x}, "transp") = @var{M2}' \ @var{x}}.
-## If @var{M1} is empty or not passed, then preconditioning is not applied.
-## The preconditioned system is theoretically equivalent to apply the
-## @code{bicg} method to the linear systems
-## @code{inv (@var{M1}) * A * inv (@var{M2}) * @var{y} = inv
-## (@var{M1}) * @var{b}} and
-## @code{inv (@var{M2'}) * A' * inv (@var{M1'}) * @var{z} =
-## inv (@var{M2'}) * @var{b}}.
-##
## @item @var{x0} the initial guess, if not given or set to [] the default
## value @code{zeros (size (b))} is used.
## @end itemize
##
-## The arguments which follow @var{x0} are treated as parameters, and passed in
-## a proper way to any of the functions (@var{A} or @var{P}) which are passed
-## to @code{bicg}.
+## @var{A} can be passed as a matrix or as a function handle or inline function
+## @code{f} such that @code{f(x, "notransp") = A*x} and
+## @code{f(x, "transp") = A'*x}.
##
-## The output parameters are:
+## The preconditioner @var{P} is given as @code{P = M1 * M2}. Both @var{M1}
+## and @var{M2} can be passed as a matrix or as a function handle or inline
+## function @code{g} such that @code{g(x, "notransp") = M1 \ x} or
+## @code{g(x, "notransp") = M2 \ x} and @code{g(x, "transp") = M1' \ x} or
+## @code{g(x, "transp") = M2' \ x}.
##
-## @itemize
+## If called with more than one output parameter
##
-## @item @var{x} is the approximation computed. If the method doesn't
-## converge then it is the iterated with the minimum residual.
-##
+## @itemize @minus
## @item @var{flag} indicates the exit status:
##
## @itemize @minus
@@ -81,257 +54,174 @@
##
## @item 1: the maximum number of iterations was reached before convergence
##
-## @item 2: the preconditioner matrix is singular
-##
## @item 3: the algorithm reached stagnation
-##
-## @item 4: the algorithm can't continue due to a division by zero
## @end itemize
##
-## @item @var{relres} is the relative residual obtained with as @code{(@var{A}*@var{x}-@var{b}) / @code{norm(@var{b})}}.
+## (the value 2 is unused but skipped for compatibility).
##
-## @item @var{iter} is the iteration which @var{x} is computed.
+## @item @var{relres} is the final value of the relative residual.
##
-## @item @var{resvec} is a vector containing the residual at each iteration.
-## Doing @code{length(@var{resvec}) - 1} is possible to see the total number
-## of iterations performed.
+## @item @var{iter} is the number of iterations performed.
+##
+## @item @var{resvec} is a vector containing the relative residual at each
+## iteration.
## @end itemize
##
-## Let us consider a trivial problem with a tridiagonal matrix
-##
-## @example
-## @group
-## n = 20;
-## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-## toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-## sparse (1, 2, 1, 1, n) * n / 2);
-## b = A * ones (n, 1);
-## restart = 5;
-## [M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
-## M = M1 * M2;
-## Afun = @@(x, string) strcmp (string, "notransp") * (A * x) + ...
-## strcmp (string, "transp") * (A' * x);
-## Mfun = @@(x, string) strcmp (string, "notransp") * (M \ x) + ...
-## strcmp (string, "transp") * (M' \ x);
-## M1fun = @@(x, string) strcmp (string, "notransp") * (M1 \ x) + ...
-## strcmp (string, "transp") * (M1' \ x);
-## M2fun = @@(x, string) strcmp (string, "notransp") * (M2 \ x) + ...
-## strcmp (string, "transp") * (M2' \ x);;
-## @end group
-## @end example
-##
-## @sc{Example 1:} simplest usage of @code{bicg}
-##
-## @example
-## x = bicg (A, b, [], n)
-## @end example
-##
-## @sc{Example 2:} @code{bicg} with a function which computes
-## @code{@var{A} * @var{x}} and @code{@var{A'} * @var{x}}
-##
-## @example
-## x = bicg (Afun, b, [], n)
-## @end example
-##
-## @sc{Example 3:} @code{bicg} with a preconditioner matrix @var{M}
-##
-## @example
-## x = bicg (A, b, [], 1e-06, n, M)
-## @end example
-##
-## @sc{Example 4:} @code{bicg} with a function as preconditioner
-##
-## @example
-## x = bicg (Afun, b, 1e-6, n, Mfun)
-## @end example
-##
-## @sc{Example 5:} @code{bicg} with preconditioner matrices @var{M1}
-## and @var{M2}
-##
-## @example
-## x = bicg (A, b, [], 1e-6, n, M1, M2)
-## @end example
-##
-## @sc{Example 6:} @code{bicg} with functions as preconditioners
-##
-## @example
-## x = bicg (Afun, b, 1e-6, n, M1fun, M2fun)
-## @end example
-##
-## @sc{Example 7:} @code {bicg} with as input a function requiring an argument
-##
-## @example
-## @group
-## function y = Ap (A, x, string, z) # compute A^z * x or (A^z)' * x
-## y = x;
-## if (strcmp (string, "notransp"))
-## for i = 1:z
-## y = A * y;
-## endfor
-## elseif (strcmp (string, "transp"))
-## for i = 1:z
-## y = A' * y;
-## endfor
-## endif
-## endfunction
-## Apfun = @(x, string, p) Ap (A, x, string, p);
-## x = bicg (Apfun, b, [], [], [], [], [], 2);
-## @end group
-## @end example
-##
-## References:
-##
-## @enumerate
-##
-## @item @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear
-## Systems}, Second edition, 2003, SIAM
-##
-## @end enumerate
-##
-## @seealso{bicgstab, cgs, gmres, pcg, qmr, tfqmr}
+## @seealso{bicgstab, cgs, gmres, pcg, qmr}
##
## @end deftypefn
-function [x_min, flag, relres, iter_min, resvec] = ...
- bicg (A, b, tol = [], maxit = [], M1 = [], M2 = [], x0 = [], varargin)
+## Author: Sylvain Pelissier
+## Author: Carlo de Falco
- [A_input, M1_input, M2_input, is_prec, ...
- is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, M1, M2);
+function [x, flag, res1, k, resvec] = bicg (A, b, rtol, maxit, M1, M2, x0)
- [tol, maxit, x0] = __default__input__ ({1e-06, min(rows(b), 20), ...
- zeros(rows (b),1)}, tol, maxit, x0);
+ if (nargin >= 2 && isvector (full (b)))
- norm_b = norm (b, 2);
+ if (ischar (A))
+ fun = str2func (A);
+ Ax = @(x) feval (fun, x, "notransp");
+ Atx = @(x) feval (fun, x, "transp");
+ elseif (isnumeric (A) && issquare (A))
+ Ax = @(x) A * x;
+ Atx = @(x) A' * x;
+ elseif (isa (A, "function_handle"))
+ Ax = @(x) feval (A, x, "notransp");
+ Atx = @(x) feval (A, x, "transp");
+ else
+ error ("bicg: A must be a square matrix or function");
+ endif
- if (norm_b == 0) # the only (only iff det(A) == 0) solution is x = 0
- if (nargout < 2)
- printf("The right hand side vector is all zero so bicg \n")
- printf ("returned an all zero solution without iterating.\n")
+ if (nargin < 3 || isempty (rtol))
+ rtol = 1e-6;
endif
- x_min = zeros (numel (b), 1);
- flag = 0;
- relres = 0;
- iter_min = 0;
- resvec = 0;
- return
+
+ if (nargin < 4 || isempty (maxit))
+ maxit = min (rows (b), 20);
+ else
+ maxit = fix (maxit);
+ endif
+
+ if (nargin < 5 || isempty (M1))
+ M1m1x = @(x, ignore) x;
+ M1tm1x = M1m1x;
+ elseif (ischar (M1))
+ fun = str2func (M1);
+ M1m1x = @(x) feval (fun, x, "notransp");
+ M1tm1x = @(x) feval (fun, x, "transp");
+ elseif (isnumeric (M1) && ismatrix (M1))
+ M1m1x = @(x) M1 \ x;
+ M1tm1x = @(x) M1' \ x;
+ elseif (isa (M1, "function_handle"))
+ M1m1x = @(x) feval (M1, x, "notransp");
+ M1tm1x = @(x) feval (M1, x, "transp");
+ else
+ error ("bicg: preconditioner M1 must be a function or matrix");
+ endif
+
+ if (nargin < 6 || isempty (M2))
+ M2m1x = @(x, ignore) x;
+ M2tm1x = M2m1x;
+ elseif (ischar (M2))
+ fun = str2func (M2);
+ M2m1x = @(x) feval (fun, x, "notransp");
+ M2tm1x = @(x) feval (fun, x, "transp");
+ elseif (isnumeric (M2) && ismatrix (M2))
+ M2m1x = @(x) M2 \ x;
+ M2tm1x = @(x) M2' \ x;
+ elseif (isa (M2, "function_handle"))
+ M2m1x = @(x) feval (M2, x, "notransp");
+ M2tm1x = @(x) feval (M2, x, "transp");
+ else
+ error ("bicg: preconditioner M2 must be a function or matrix");
+ endif
+
+ Pm1x = @(x) M2m1x (M1m1x (x));
+ Ptm1x = @(x) M1tm1x (M2tm1x (x));
+
+ if (nargin < 7 || isempty (x0))
+ x0 = zeros (size (b));
+ endif
+
+ y = x = x0;
+ c = b;
+
+ r0 = b - Ax (x);
+ s0 = c - Atx (y);
+
+ d = Pm1x (r0);
+ f = Ptm1x (s0);
+
+ bnorm = norm (b);
+ res0 = Inf;
+
+ if (any (r0 != 0))
+
+ for k = 1:maxit
+
+ a = (s0' * Pm1x (r0)) ./ (f' * Ax (d));
+
+ x += a * d;
+ y += conj (a) * f;
+
+ r1 = r0 - a * Ax (d);
+ s1 = s0 - conj (a) * Atx (f);
+
+ beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0));
+
+ d = Pm1x (r1) + beta * d;
+ f = Ptm1x (s1) + conj (beta) * f;
+
+ r0 = r1;
+ s0 = s1;
+
+ res1 = norm (b - Ax (x)) / bnorm;
+ if (res1 < rtol)
+ flag = 0;
+ if (nargout < 2)
+ printf ("bicg converged at iteration %i ", k);
+ printf ("to a solution with relative residual %e\n", res1);
+ endif
+ break;
+ endif
+
+ if (res0 <= res1)
+ flag = 3;
+ printf ("bicg stopped at iteration %i ", k);
+ printf ("without converging to the desired tolerance %e\n", rtol);
+ printf ("because the method stagnated.\n");
+ printf ("The iterate returned (number %i) ", k-1);
+ printf ("has relative residual %e\n", res0);
+ break
+ endif
+ res0 = res1;
+ if (nargout > 4)
+ resvec(k) = res0;
+ endif
+ endfor
+
+ if (k == maxit)
+ flag = 1;
+ printf ("bicg stopped at iteration %i ", maxit);
+ printf ("without converging to the desired tolerance %e\n", rtol);
+ printf ("because the maximum number of iterations was reached. ");
+ printf ("The iterate returned (number %i) has ", maxit);
+ printf ("relative residual %e\n", res1);
+ endif
+
+ else
+ flag = 0;
+ if (nargout < 2)
+ printf ("bicg converged after 0 interations\n");
+ endif
+ endif
+
+ else
+ print_usage ();
endif
- if (is_matrix_A && is_matrix_M1 && is_matrix_M2)
- A = A_input;
- M1 = M1_input;
- M2 = M2_input;
- [x_min, flag, iter, iter_min, resvec] = ...
- bicg_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
- else
- if (is_matrix_A)
- A = @(x, trans) A_sub (A_input, x, trans);
- endif
- if (is_prec)
- if (is_matrix_M1)
- M1 = @(x, trans) M_sub (M1_input, x, trans);
- endif
- if (!isempty (M2_input) && is_matrix_M2)
- M2 = @(x, trans) M_sub (M2_input, x, trans);
- endif
- endif
- [x_min, flag, iter, iter_min, resvec] = ...
- bicg_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, varargin{:});
- endif
+endfunction;
- if (flag == 2)
- relres = 1;
- else
- relres = resvec (iter_min + 1) / norm_b;
- endif
-
- if ((flag == 1) && (relres <= tol))
- flag = 0;
- endif
-
- if (nargout < 2)
- switch (flag)
- case {0}
- printf ("bicg converged at iteration %i ", iter_min);
- printf ("to a solution with relative residual %e\n", relres);
- case {1}
- printf ("bicg stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the maximum number of iterations was reached. ");
- printf ("The iterate returned (number %i) has ", iter_min);
- printf ("relative residual %e\n", relres);
- case {2}
- printf ("bicg stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the preconditioner matrix is singular.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {3}
- printf ("bicg stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method stagnated.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {4}
- printf ("bicg stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method can't continue.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- endswitch
- endif
-
-endfunction
-
-function y = A_sub (A_input, x, trans)
- if (strcmp (trans, "transp"))
- y = A_input' * x;
- else
- y = A_input * x;
- endif
-endfunction
-
-function y = M_sub (M_input, x, trans)
- if (strcmp (trans, "transp"))
- y = M_input' \ x;
- else
- y = M_input \ x;
- endif
-endfunction
-
-
-%!test
-%! ## Check that all the subscripts works
-%! A = magic (5);
-%! b = rand (5, 1);
-%! M1 = diag (sqrt (diag(A)));
-%! M2 = M1;
-%! Afun = @(z, string) strcmp(string,"notransp")*(A*z) + ...
-%! strcmp(string,"transp")*(A'*z);
-%! M1_fun = @(z, string) strcmp(string,"notransp")*(M1 \ z) + ...
-%! strcmp(string,"transp")*(M1' \ z);
-%! M2_fun = @(z, string) strcmp(string,"notransp")*(M2 \ z) + ...
-%! strcmp(string,"transp")*(M2' \ z);
-%! [x, flag] = bicg (A,b);
-%! assert(flag, 0);
-%! [x, flag] = bicg (A, b, [], [], M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicg (A, b, [], [], M1_fun, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = bicg (A, b,[],[], M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicg (A, b,[],[], M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = bicg (Afun, b);
-%! assert(flag, 0);
-%! [x, flag] = bicg (Afun, b,[],[], M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicg (Afun, b,[],[], M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicg (Afun, b,[],[], M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = bicg (Afun, b,[],[], M1_fun, M2_fun);
-%! assert(flag, 0);
%!test
%! n = 100;
@@ -375,64 +265,3 @@
%! [x, flag, relres, iter, resvec] = bicg (A, b, rtol, [], diag (diag (A)));
%! assert (x, ones (size (b)), 1e-7);
-%!test
-%! ## Check that if the preconditioner is singular, the method doesn't work
-%! A = rand(5);
-%! b = ones(5,1);
-%! M = ones(5);
-%! [x,flag] = bicg (A,b,[],[],M);
-%! assert(flag, 2)
-
-%!test
-%! ## If A singular, the algorithm doesn't work due to division by zero
-%! A = ones(5);
-%! b = [1:5]';
-%! [x,flag] = bicg (A,b);
-%! assert(flag, 4)
-
-%!test
-%! ## test for a complex linear system
-%! A = rand(5) + 1i*rand(5);
-%! b = sum(A,2);
-%! [x,flag] = bicg(A,b);
-%! assert(flag,0)
-
-%!demo # simplest use
-%! n = 20;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.1 * eye (n));
-%! M = M1 * M2;
-%! x = bicg (A, b, [], n);
-%! function y = Ap (A, x, string, z) # compute A^z * x or (A^z)' * x
-%! y = x;
-%! if (strcmp (string, "notransp"))
-%! for i = 1:z
-%! y = A * y;
-%! endfor
-%! elseif (strcmp (string, "transp"))
-%! for i = 1:z
-%! y = A' * y;
-%! endfor
-%! endif
-%! endfunction
-%! Afun = @(x, string) Ap (A, x, string, 1);
-%! x = bicg (Afun, b, [], n);
-%! x = bicg (A, b, 1e-6, n, M);
-%! x = bicg (A, b, 1e-6, n, M1, M2);
-%! function y = Mfun(M, x, string)
-%! if (strcmp (string, "notransp"))
-%! y = M \ x;
-%! else
-%! y = M' \ x;
-%! endif
-%! endfunction
-%! M1fun = @(x, string) Mfun (M, x, string);
-%! x = bicg (Afun, b, 1e-6, n, M1fun);
-%! M1fun = @(x, string) Mfun (M1, x, string);
-%! M2fun = @(x, string) Mfun (M2, x, string);
-%! x = bicg (Afun, b, 1e-6, n, M1fun, M2fun);
-%! Afun = @(x, string, p) Ap (A, x, string, p);
-%! x = bicg (Afun, b, [], 2*n, [], [], [], 2); # solution of A^2 * x = b
diff --git a/scripts/sparse/bicgstab.m b/scripts/sparse/bicgstab.m
--- a/scripts/sparse/bicgstab.m
+++ b/scripts/sparse/bicgstab.m
@@ -1,6 +1,5 @@
-## Copyright (C) 2008-2015 Radek Salac
+## Copyright (C) 2008-2016 Radek Salac
## Copyright (C) 2012 Carlo de Falco
-## Copyright (C) 2016 Cristiano Dorigo, Octave Arena
##
## This file is part of Octave.
##
@@ -19,56 +18,33 @@
## .
## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
-## @deftypefnx {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
-## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, @dots{})
+## @deftypefn {} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## @deftypefnx {} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
+## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, @dots{})
## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative
## method.
##
-## The input parameters are:
-##
## @itemize @minus
-##
-## @item @var{A} is the square matrix of the linear system. @var{A}
-## can be passed
-## as a matrix or as a function handle or inline function @code{f} such that
-## @code{f(@var{x}) = @var{A} * @var{x}}.
-##
-## @item @var{b} is the right hand side vector. It must be a column vector
-## with the same number of rows as @var{A}.
-##
## @item @var{rtol} is the relative tolerance, if not given or set to [] the
## default value 1e-6 is used.
##
## @item @var{maxit} the maximum number of outer iterations, if not given or
## set to [] the default value @code{min (20, numel (b))} is used.
##
-## @item @var{M1}, @var{M2} are the preconditioners. The preconditioner
-## @var{P} is given as @code{@var{P} = @var{M1} * @var{M2}}.
-## Both @var{M1} and @var{M2} can be passed as a matrix or as a function
-## handle or inline function @code{g} such that
-## @code{g(@var{x}) = @var{M1} \ @var{x}} or
-## @code {g(@var{x}) = @var{M2} \ @var{x}}.
-## The techinque used is the right preconditioning, i.e. it is
-## solved @code{@var{A} * inv (@var{P}) * @var{y} = @var{b}} and then
-## @code{@var{x} = inv (@var{P}) * @var{y}}.
-##
## @item @var{x0} the initial guess, if not given or set to [] the default
-## value @code{zeros (size (@var{b}))} is used.
-##
+## value @code{zeros (size (b))} is used.
## @end itemize
##
-## The arguments which follow @var{x0} are treated as parameters, and passed in
-## a proper way to any of the functions (@var{A} or @var{P}) which are passed
-## to @code{bicstab}.
+## @var{A} can be passed as a matrix or as a function handle or inline
+## function @code{f} such that @code{f(x) = A*x}.
##
-## The output parameters are:
+## The preconditioner @var{P} is given as @code{P = M1 * M2}. Both @var{M1}
+## and @var{M2} can be passed as a matrix or as a function handle or inline
+## function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}.
+##
+## If called with more than one output parameter
##
## @itemize @minus
-##
-## @item @var{x} is the approximation computed. If the method doesn't
-## converge then it is the iterated with the minimum residual.
-##
## @item @var{flag} indicates the exit status:
##
## @itemize @minus
@@ -76,251 +52,161 @@
##
## @item 1: the maximum number of iterations was reached before convergence
##
-## @item 2: the preconditioner matrix is singular
-##
## @item 3: the algorithm reached stagnation
-##
-## @item 4: the algorithm can't continue due to a division by zero
## @end itemize
##
-## @item @var{relres} is the relative residual obtained with as
-## @code{(@var{A}*@var{x}-@var{b}) / @code{norm(@var{b})}}.
+## (the value 2 is unused but skipped for compatibility).
##
-## @item @var{iter} is the (possibily half) iteration which @var{x} is
-## computed. If it is an half iteration then it is @code{@var{iter} + 0.5}
+## @item @var{relres} is the final value of the relative residual.
##
-## @item @var{resvec} is a vector containing the residual of each half and
-## total iteration (There are also the half iterations since @var{x} is
-## computed in two steps at each iteration).
-## Doing @code{(length(@var{resvec}) - 1) / 2} is possible to see the
-## total number of (total) iterations performed.
+## @item @var{iter} is the number of iterations performed.
##
+## @item @var{resvec} is a vector containing the relative residual at each
+## iteration.
## @end itemize
##
-## Let us consider a trivial problem with a tridiagonal matrix
-##
-## @example
-## @group
-## n = 20;
-## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-## toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-## sparse (1, 2, 1, 1, n) * n / 2);
-## b = A * ones (n, 1);
-## restart = 5;
-## [M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
-## M = M1 * M2;
-## Afun = @@(x) A * x;
-## Mfun = @@(x) M \ x;
-## M1fun = @@(x) M1 \ x;
-## M2fun = @@(x) M2 \ x;
-## @end group
-## @end example
-##
-## @sc{Example 1:} simplest usage of @code{bicgstab}
-##
-## @example
-## x = bicgstab (A, b, [], n)
-## @end example
-##
-## @sc{Example 2:} @code{bicgstab} with a function which computes
-## @code{@var{A} * @var{x}}
-##
-## @example
-## x = bicgstab (Afun, b, [], n)
-## @end example
-##
-## @sc{Example 3:} @code{bicgstab} with a preconditioner matrix @var{M}
-##
-## @example
-## x = bicgstab (A, b, [], 1e-06, n, M)
-## @end example
-##
-## @sc{Example 4:} @code{bicgstab} with a function as preconditioner
-##
-## @example
-## x = bicgstab (Afun, b, 1e-6, n, Mfun)
-## @end example
-##
-## @sc{Example 5:} @code{bicgstab} with preconditioner matrices @var{M1}
-## and @var{M2}
-##
-## @example
-## x = bicgstab (A, b, [], 1e-6, n, M1, M2)
-## @end example
-##
-## @sc{Example 6:} @code{bicgstab} with functions as preconditioners
-##
-## @example
-## x = bicgstab (Afun, b, 1e-6, n, M1fun, M2fun)
-## @end example
-##
-## @sc{Example 7:} @code {bicgstab} with as input a function requiring
-## an argument
-##
-## @example
-## @group
-## function y = Ap (A, x, z) # compute A^z * x
-## y = x;
-## for i = 1:z
-## y = A * y;
-## endfor
-## endfunction
-## Apfun = @(x, string, p) Ap (A, x, string, p);
-## x = bicgstab (Apfun, b, [], [], [], [], [], 2);
-## @end group
-## @end example
-##
-## @sc{Example 8:} explicit example to show that @code{bicgstab} uses a
-## right preconditioner
-##
-## @example
-## @group
-## [M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
-## M = M1 * M2;
-##
-## ## reference solution computed by bicgstab after one iteration
-## [x_ref, fl] = bicgstab (A, b, [], 1, M)
-##
-## ## rigth preconditioning
-## [y, fl] = bicgstab (A / M, b, [], 1)
-## x = M \ y # compare x and x_ref
-##
-## @end group
-## @end example
-##
-## References:
-##
-## @enumerate
-##
-## @item @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear
-## Systems}, Second edition, 2003, SIAM
-##
-## @end enumerate
-##
-## @seealso{bicg, cgs, gmres, pcg, qmr, tfqmr}
+## @seealso{bicg, cgs, gmres, pcg, qmr}
##
## @end deftypefn
-function [x_min, flag, relres, iter_min, resvec] = ...
- bicgstab (A, b, tol = [], maxit = [], M1 = [], M2 = [], ...
- x0 = [], varargin)
+function [x, flag, relres, iter, resvec] = bicgstab (A, b, rtol, maxit,
+ M1, M2, x0)
- ## Check consistency and type of A, M1, M2
- [A, M1, M2, is_prec, is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, M1, M2);
-
- # Check if input tol are empty (set them to default if necessary)
- [tol, maxit, x0] = __default__input__ ({1e-06, min(rows(b), 20), ...
- zeros(rows(b), 1)}, tol, maxit, x0);
-
- norm_b = norm (b, 2);
- if (norm_b == 0)
- if (nargout < 2)
- printf("The right hand side vector is all zero so bicgstab \n")
- printf ("returned an all zero solution without iterating.\n")
- endif
- x_min = zeros (numel (b), 1);
- iter_ min = 0;
- flag = 0;
- resvec = 0;
- relres = 0;
- return
+ if (nargin < 2 || nargin > 7 || ! isvector (full (b)))
+ print_usage ();
endif
- ## different script for different types of A, M1, M2
-
- if (is_matrix_A && is_matrix_M1 && is_matrix_M2)
- [x_min, flag, iter, iter_min, resvec] = ...
- bicgstab_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
+ if (ischar (A))
+ A = str2func (A);
+ elseif (isnumeric(A) && issquare (A))
+ Ax = @(x) A * x;
+ elseif (isa (A, "function_handle"))
+ Ax = @(x) feval (A, x);
else
- if (is_matrix_A)
- A = @(x) A * x;
- endif
- if (is_prec)
- if (is_matrix_M1)
- M1 = @(x) M1 \ x;
- endif
- if (!isempty (M2) && is_matrix_M2)
- M2 = @(x) M2 \ x;
- endif
- endif
- [x_min, flag, iter, iter_min, resvec] = ...
- bicgstab_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, ...
- varargin{:});
- endif
- relres = resvec (iter_min + 1) / norm_b; ## I set the relative residual
- iter /= 2;
- iter_min /= 2;
-
- if (flag == 1) && (relres <= tol)
- flag = 0;
+ error ("bicgstab: A must be a square matrix or function");
endif
- ## output strings to print when the outputs requested are less than 2
+ if (nargin < 3 || isempty (rtol))
+ rtol = 1e-6;
+ endif
+
+ if (nargin < 4 || isempty (maxit))
+ maxit = min (rows (b), 20);
+ endif
+
+ if (nargin < 5 || isempty (M1))
+ M1m1x = @(x) x;
+ elseif (ischar (M1))
+ M1m1x = str2func (M1);
+ elseif (isnumeric(M1) && ismatrix (M1))
+ M1m1x = @(x) M1 \ x;
+ elseif (isa (M1, "function_handle"))
+ M1m1x = @(x) feval (M1, x);
+ else
+ error ("bicgstab: preconditioner M1 must be a function or matrix");
+ endif
+
+ if (nargin < 6 || isempty (M2))
+ M2m1x = @(x) x;
+ elseif (ischar (M2))
+ M2m1x = str2func (M2);
+ elseif (isnumeric(M2) && ismatrix (M2))
+ M2m1x = @(x) M2 \ x;
+ elseif (isa (M2, "function_handle"))
+ M2m1x = @(x) feval (M2, x);
+ else
+ error ("bicgstab: preconditioner M2 must be a function or matrix");
+ endif
+
+ precon = @(x) M2m1x (M1m1x (x));
+
+ if (nargin < 7 || isempty (x0))
+ x0 = zeros (size (b));
+ endif
+
+ ## specifies initial estimate x0
+ if (nargin < 7)
+ x = zeros (rows (b), 1);
+ else
+ x = x0;
+ endif
+
+ norm_b = norm (b);
+
+ res = b - Ax (x);
+ rr = res;
+
+ ## Vector of the residual norms for each iteration.
+ resvec = norm (res) / norm_b;
+
+ ## Default behavior we don't reach tolerance rtol within maxit iterations.
+ flag = 1;
+
+ for iter = 1:maxit
+ rho_1 = rr' * res;
+
+ if (iter == 1)
+ p = res;
+ else
+ beta = (rho_1 / rho_2) * (alpha / omega);
+ p = res + beta * (p - omega * v);
+ endif
+
+ phat = precon (p);
+
+ v = Ax (phat);
+ alpha = rho_1 / (rr' * v);
+ s = res - alpha * v;
+
+ shat = precon (s);
+
+ t = Ax (shat);
+ omega = (s' * t) / (t' * t);
+ x += alpha * phat + omega * shat;
+ res = s - omega * t;
+ rho_2 = rho_1;
+
+ relres = norm (res) / norm_b;
+ resvec = [resvec; relres];
+
+ if (relres <= rtol)
+ ## We reach tolerance rtol within maxit iterations.
+ flag = 0;
+ break;
+ elseif (resvec(end) == resvec(end - 1))
+ ## The method stagnates.
+ flag = 3;
+ break;
+ endif
+ endfor
+
if (nargout < 2)
- switch (flag)
- case {0}
- printf ("bicgstab converged at iteration %i ", iter_min);
- printf ("to a solution with relative residual %e\n", relres);
- case {1}
- printf ("bicgstab stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the maximum number of iterations was reached.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {2}
- printf ("bicgstab stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the preconditioner matrix is singular.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {3}
- printf ("bicgstab stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method stagnated.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {4}
- printf ("bicgstab stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method can't continue.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- endswitch
+ if (flag == 0)
+ printf ("bicgstab converged at iteration %i ", iter);
+ printf ("to a solution with relative residual %e\n", relres);
+ elseif (flag == 3)
+ printf ("bicgstab stopped at iteration %i ", iter);
+ printf ("without converging to the desired tolerance %e\n", rtol);
+ printf ("because the method stagnated.\n");
+ printf ("The iterate returned (number %i) ", iter);
+ printf ("has relative residual %e\n", relres);
+ else
+ printf ("bicgstab stopped at iteration %i ", iter);
+ printf ("without converging to the desired toleranc %e\n", rtol);
+ printf ("because the maximum number of iterations was reached.\n");
+ printf ("The iterate returned (number %i) ", iter);
+ printf ("has relative residual %e\n", relres);
+ endif
endif
endfunction
-%!test
-%! ## Check that all the subscripts works
-%! A = magic (5);
-%! b = rand (5, 1);
-%! M1 = diag (sqrt (diag(A)));
-%! M2 = M1;
-%! maxit = 10;
-%! Afun = @(z) A*z;
-%! M1_fun = @(z) M1 \ z;
-%! M2_fun = @(z) M2 \ z;
-%! [x, flag] = bicgstab (A,b);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (A, b, [], maxit, M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (A, b, [], maxit, M1_fun, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (A, b, [], maxit, M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (A, b, [], maxit, M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (Afun, b);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (Afun, b, [], maxit, M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (Afun, b, [], maxit, M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (Afun, b, [], maxit, M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = bicgstab (Afun, b, [], maxit, M1_fun, M2_fun);
-%! assert(flag, 0);
+
+%!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)
%!shared A, b, n, M1, M2
%!
@@ -357,66 +243,8 @@
%! assert (x, ones (size (b)), 1e-7);
%!test
-%! ## bicgstab solves complex linear systems
%! A = [1 + 1i, 1 + 1i; 2 - 1i, 2 + 1i];
%! b = A * [1; 1];
%! [x, flag, relres, iter, resvec] = bicgstab (A, b);
%! assert (x, [1; 1], 1e-6);
-%!test
-%! ## test with a non symmetric matrix
-%! A = diag(1:50);
-%! A (1,50) = 10000;
-%! b = ones (50,1);
-%! [x, flag, relres, iter, resvec] = bicgstab (A, b, [], 100);
-%! assert (flag, 0)
-%! assert (x, A\b, 1e-05)
-%! ## test that bicgstab detects a singular preconditioner
-%! M = ones(50);
-%! M(1,1) = 0;
-%! [x, flag] = bicgstab (A, b, [], 100, M);
-%! assert(flag, 2)
-
-%!demo # simplest use
-%! n = 20;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.1 * eye (n));
-%! M = M1 * M2;
-%! x = bicgstab (A, b, [], n);
-%! Afun = @(x) A * x;
-%! x = bicgstab (Afun, b, [], n);
-%! x = bicgstab (A, b, 1e-6, n, M);
-%! x = bicgstab (A, b, 1e-6, n, M1, M2);
-%! Mfun = @(z) M \ z;
-%! x = bicgstab (Afun, b, 1e-6, n, Mfun);
-%! M1fun = @(z) M1 \ z;
-%! M2fun = @(z) M2 \ z;
-%! x = bicgstab (Afun, b, 1e-6, n, M1fun, M2fun);
-%! function y = Ap (A, x, z) # compute A^z * x or (A^z)' * x
-%! y = x;
-%! for i = 1:z
-%! y = A * y;
-%! endfor
-%! endfunction
-%! Afun = @(x, p) Ap (A, x, p);
-%! x = bicgstab (Afun, b, [], 2*n, [], [], [], 2); # solution of A^2 * x = b
-
-%!demo
-%! n = 10;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
-%! M = M1 * M2;
-%!
-%! ## reference solution computed by bicgstab after one iteration
-%! [x_ref, fl] = bicgstab (A, b, [], 1, M);
-%! x_ref
-%!
-%! ## right preconditioning
-%! [y, fl] = bicgstab (A / M, b, [], 1);
-%! x = M \ y # compare x and x_ref
diff --git a/scripts/sparse/cgs.m b/scripts/sparse/cgs.m
--- a/scripts/sparse/cgs.m
+++ b/scripts/sparse/cgs.m
@@ -1,6 +1,5 @@
-## Copyright (C) 2008-2015 Radek Salac
+## Copyright (C) 2008-2016 Radek Salac
## Copyright (C) 2012 Carlo de Falco
-## Copyright (C) 2016 Cristiano Dorigo, Octave Arena
##
## This file is part of Octave.
##
@@ -19,50 +18,33 @@
## .
## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} cgs (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
-## @deftypefnx {Function File} {@var{x} =} cgs (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
-## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} cgs (@var{A}, @var{b}, @dots{})
+## @deftypefn {} {@var{x} =} cgs (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## @deftypefnx {} {@var{x} =} cgs (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
+## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} cgs (@var{A}, @var{b}, @dots{})
## Solve @code{A x = b}, where @var{A} is a square matrix, using the
## Conjugate Gradients Squared method.
##
-## The input arguments are:
-##
## @itemize @minus
-##
-## @item @var{A} it is a square matrix. It can be passed as a matrix or as a function handle or inline function @code{f} such that @code{f(x) = A*x}.
-##
-## @item @var{b} is the right hand side vector. It must be a column vector
-## with same number of rows of @var{A}.
-##
## @item @var{rtol} is the relative tolerance, if not given or set to [] the
## default value 1e-6 is used.
##
## @item @var{maxit} the maximum number of outer iterations, if not given or
## set to [] the default value @code{min (20, numel (b))} is used.
##
-## @item @var{M1}, @var{M2} are the preconditioners. The preconditioner
-## matrix is given as @code{P = M1 * M2}. Both @var{M1}
-## and @var{M2} can be passed as a matrix or as a function handle or inline
-## function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}.
-## If M1 is empty or not passed then no preconditioners are applied.
-## The techinque used is the right preconditioning, i.e. it is solved
-## @code{@var{A}*inv(@var{P})*y = b} and then @code{@var{x} = inv(@var{P})*y}.
-##
## @item @var{x0} the initial guess, if not given or set to [] the default
## value @code{zeros (size (b))} is used.
## @end itemize
##
-## The arguments which follow @var{x0} are treated as parameters, and passed in
-## a proper way to any of the functions (@var{A} or @var{P}) which are passed
-## to @code{cgs}.
+## @var{A} can be passed as a matrix or as a function handle or inline
+## function @code{f} such that @code{f(x) = A*x}.
##
-## The output parameters are:
+## The preconditioner @var{P} is given as @code{P = M1 * M2}. Both @var{M1}
+## and @var{M2} can be passed as a matrix or as a function handle or inline
+## function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}.
+##
+## If called with more than one output parameter
##
## @itemize @minus
-##
-## @item @var{x} is the approximation computed. If the method doesn't
-## converge then it is the iterated with the minimum residual.
-##
## @item @var{flag} indicates the exit status:
##
## @itemize @minus
@@ -70,290 +52,147 @@
##
## @item 1: the maximum number of iterations was reached before convergence
##
-## @item 2: the preconditioner matrix is singular
-##
## @item 3: the algorithm reached stagnation
-##
-## @item 4: the algorithm can't continue due to a division by zero
## @end itemize
##
-## @item @var{relres} is the relative residual obtained with as
-## @code{(@var{A}*@var{x}-@var{b}) / @code{norm(@var{b})}}.
+## (the value 2 is unused but skipped for compatibility).
##
-## @item @var{iter} is the iteration which @var{x} is computed.
+## @item @var{relres} is the final value of the relative residual.
##
-## @item @var{resvec} is a vector containing the residual at each iteration.
-## Doing @code{length(@var{resvec}) - 1} is possible to see the total number
-## of iterations performed.
+## @item @var{iter} is the number of iterations performed.
+##
+## @item @var{resvec} is a vector containing the relative residual at
+## each iteration.
## @end itemize
##
-## Let us consider a trivial problem with a tridiagonal matrix
-##
-## @example
-## @group
-## n = 20;
-## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-## toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-## sparse (1, 2, 1, 1, n) * n / 2);
-## b = A * ones (n, 1);
-## restart = 5;
-## [M1, M2] = ilu (A); # in this tridiagonal case it corresponds to chol (A)'
-## M = M1 * M2;
-## Afun = @@(x) A * x;
-## Mfun = @@(x) M \ x;
-## M1fun = @@(x) M1 \ x;
-## M2fun = @@(x) M2 \ x;
-## @end group
-## @end example
-##
-## @sc{Example 1:} simplest usage of @code{cgs}
-##
-## @example
-## x = cgs (A, b, [], n)
-## @end example
-##
-## @sc{Example 2:} @code{cgs} with a function which computes
-## @code{@var{A} * @var{x}}
-##
-## @example
-## x = cgs (Afun, b, [], n)
-## @end example
-##
-## @sc{Example 3:} @code{cgs} with a preconditioner matrix @var{M}
-##
-## @example
-## x = cgs (A, b, [], 1e-06, n, M)
-## @end example
-##
-## @sc{Example 4:} @code{cgs} with a function as preconditioner
-##
-## @example
-## x = cgs (Afun, b, 1e-6, n, Mfun)
-## @end example
-##
-## @sc{Example 5:} @code{cgs} with preconditioner matrices @var{M1}
-## and @var{M2}
-##
-## @example
-## x = cgs (A, b, [], 1e-6, n, M1, M2)
-## @end example
-##
-## @sc{Example 6:} @code{cgs} with functions as preconditioners
-##
-## @example
-## x = cgs (Afun, b, 1e-6, n, M1fun, M2fun)
-## @end example
-##
-## @sc{Example 7:} @code {cgs} with as input a function requiring an argument
-##
-## @example
-## @group
-## function y = Ap (A, x, z) # compute A^z * x
-## y = x;
-## for i = 1:z
-## y = A * y;
-## endfor
-## endfunction
-## Apfun = @(x, string, p) Ap (A, x, string, p);
-## x = cgs (Apfun, b, [], [], [], [], [], 2);
-## @end group
-## @end example
-##
-## @sc{Example 8:} explicit example to show that @code{cgs} uses a
-## right preconditioner
-##
-## @example
-## @group
-## [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
-## M = M1 * M2;
-##
-## ## reference solution computed by cgs after one iteration
-## [x_ref, fl] = cgs (A, b, [], 1, M)
-##
-## ## rigth preconditioning
-## [y, fl] = cgs (A / M, b, [], 1)
-## x = M \ y # compare x and x_ref
-##
-## @end group
-## @end example
-##
-## References:
-##
-## @enumerate
-##
-## @item @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear Systems},
-## Second edition, 2003, SIAM
-##
-## @end enumerate
-##
-## @seealso{pcg, bicgstab, bicg, gmres, qmr, tfqmr}
+## @seealso{pcg, bicgstab, bicg, gmres, qmr}
## @end deftypefn
-function [x_min, flag, relres, iter_min, resvec] = ...
- cgs (A, b, tol = [], maxit = [], M1 = [] , M2 = [], x0 = [], varargin)
+function [x, flag, relres, iter, resvec] = cgs (A, b, rtol, maxit, M1, M2, x0)
- [A, M1, M2, is_prec, is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, M1, M2);
+ if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
- [tol, maxit, x0] = __default__input__ ({1e-06, min( rows(b), 20), ...
- zeros(size(b))}, tol, maxit, x0);
+ if (ischar (A))
+ A = str2func (A);
+ elseif (isnumeric (A) && issquare (A))
+ Ax = @(x) A * x;
+ elseif (isa (A, "function_handle"))
+ Ax = @(x) feval (A, x);
+ else
+ error ("cgs: A must be a square matrix or function");
+ endif
- norm_b = norm (b, 2);
- if (norm_b == 0)
- if (nargout < 2)
- printf("The right hand side vector is all zero so cgs \n")
- printf ("returned an all zero solution without iterating.\n")
+ if (nargin < 3 || isempty (rtol))
+ rtol = 1e-6;
endif
- x_min = zeros (numel (b), 1);
- iter_ min = 0;
- flag = 0;
- resvec = 0;
- relres = 0;
- return
- endif
- resvec = zeros (maxit, 1); # Preallocation of resvec
+ if (nargin < 4 || isempty (maxit))
+ maxit = min (rows (b), 20);
+ endif
- if (is_matrix_A && is_matrix_M1 && is_matrix_M2)
- [x_min, flag, iter, iter_min, resvec] = ...
- cgs_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
- else
- if (is_matrix_A)
- A = @(x) A * x;
+ if (nargin < 5 || isempty (M1))
+ M1m1x = @(x) x;
+ elseif (ischar (M1))
+ M1m1x = str2func (M1);
+ elseif (isnumeric (M1) && ismatrix (M1))
+ M1m1x = @(x) M1 \ x;
+ elseif (isa (M1, "function_handle"))
+ M1m1x = @(x) feval (M1, x);
+ else
+ error ("cgs: preconditioner M1 must be a function or matrix");
endif
- if (is_prec)
- if (is_matrix_M1)
- M1 = @(x) M1 \ x;
+
+ if (nargin < 6 || isempty (M2))
+ M2m1x = @(x) x;
+ elseif (ischar (M2))
+ M2m1x = str2func (M2);
+ elseif (isnumeric (M2) && ismatrix (M2))
+ M2m1x = @(x) M2 \ x;
+ elseif (isa (M2, "function_handle"))
+ M2m1x = @(x) feval (M2, x);
+ else
+ error ("cgs: preconditioner M2 must be a function or matrix");
+ endif
+
+ precon = @(x) M2m1x (M1m1x (x));
+
+ if (nargin < 7 || isempty (x0))
+ x0 = zeros (size (b));
+ endif
+
+
+ x = x0;
+
+ res = b - Ax (x);
+ norm_b = norm (b);
+ ## Vector of the residual norms for each iteration.
+ resvec = norm (res) / norm_b;
+ ro = 0;
+ ## Default behavior we don't reach tolerance rtol within maxit iterations.
+ flag = 1;
+ for iter = 1:maxit
+
+ z = precon (res);
+
+ ## Cache.
+ ro_old = ro;
+ ro = res' * z;
+ if (iter == 1)
+ p = z;
+ else
+ beta = ro / ro_old;
+ p = z + beta * p;
endif
- if (!isempty (M2) && is_matrix_M2)
- M2 = @(x) M2 \ x;
+ ## Cache.
+ q = Ax (p);
+ alpha = ro / (p' * q);
+ x += alpha * p;
+
+ res -= alpha * q;
+ relres = norm (res) / norm_b;
+ resvec = [resvec; relres];
+
+ if (relres <= rtol)
+ ## We reach tolerance rtol within maxit iterations.
+ flag = 0;
+ break
+ elseif (resvec(end) == resvec(end - 1))
+ ## The method stagnates.
+ flag = 3;
+ break
+ endif
+ endfor
+
+ if (nargout < 1)
+ if (flag == 0)
+ printf ("cgs converged at iteration %i to a solution with relative residual %e\n",
+ iter, relres);
+ elseif (flag == 3)
+ printf (["cgs 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, rtol, iter, relres);
+ else
+ printf (["cgs stopped at iteration %i without converging to the desired tolerance %e\n",
+ "because the maximum number of iterations was reached.\n",
+ "The iterate returned (number %i) has relative residual %e\n"],
+ iter, rtol, iter, relres);
endif
endif
- [x_min, flag, iter, iter_min, resvec] = ...
- cgs_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, ...
- varargin{:});
- endif
- relres = resvec (iter_min + 1) / norm_b;
- if (relres <= tol) && (flag = 1)
- flag = 0;
- endif
-
- if (nargout < 2)
- switch (flag)
- case {0}
- printf ("cgs converged at iteration %i ", iter_min);
- printf ("to a solution with relative residual %e\n", relres);
- case {1}
- printf ("cgs stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the maximum number of iterations was reached.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {2}
- printf ("cgs stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the preconditioner matrix is singular.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {3}
- printf ("cgs stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method stagnated.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {4}
- printf ("cgs stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method can't continue.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- endswitch
+ else
+ print_usage ();
endif
endfunction
+
%!demo
%! % Solve system of A*x=b
%! A = [5 -1 3;-1 2 -2;3 -2 3];
%! b = [7;-1;4];
%! [a,b,c,d,e] = cgs (A,b)
-%!demo # simplest use
-%! n = 20;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.1 * eye (n));
-%! M = M1 * M2;
-%! x = cgs (A, b, [], n);
-%! Afun = @(x) A * x;
-%! x = cgs (Afun, b, [], n);
-%! x = cgs (A, b, 1e-6, n, M);
-%! x = cgs (A, b, 1e-6, n, M1, M2);
-%! Mfun = @(z) M \ z;
-%! x = cgs (Afun, b, 1e-6, n, Mfun);
-%! M1fun = @(z) M1 \ z;
-%! M2fun = @(z) M2 \ z;
-%! x = cgs (Afun, b, 1e-6, n, M1fun, M2fun);
-%! function y = Ap (A, x, z) # compute A^z * x or (A^z)' * x
-%! y = x;
-%! for i = 1:z
-%! y = A * y;
-%! endfor
-%! endfunction
-%! Afun = @(x, p) Ap (A, x, p);
-%! x = cgs (Afun, b, [], 2*n, [], [], [], 2); # solution of A^2 * x = b
-
-%!demo
-%! n = 10;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
-%! M = M1 * M2;
-%!
-%! ## reference solution computed by cgs after one iteration
-%! [x_ref, fl] = cgs (A, b, [], 1, M);
-%! x_ref
-%!
-%! ## right preconditioning
-%! [y, fl] = cgs (A / M, b, [], 1);
-%! x = M \ y # compare x and x_ref
-
-%!test
-%! ## Check that all the subscripts works
-%! A = magic (5);
-%! b = rand (5,1);
-%! M1 = diag (sqrt (diag (A)));
-%! M2 = M1;
-%! maxit = 10;
-%! Afun = @(z) A * z;
-%! M1_fun = @(z) M1 \ z;
-%! M2_fun = @(z) M2 \ z;
-%! [x, flag] = cgs (A,b);
-%! assert(flag, 0);
-%! [x, flag] = cgs (A, b, [], maxit, M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = cgs (A, b, [], maxit, M1_fun, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = cgs (A, b, [], maxit, M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = cgs (A, b, [], maxit, M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = cgs (Afun, b);
-%! assert(flag, 0);
-%! [x, flag] = cgs (Afun, b, [], maxit, M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = cgs (Afun, b, [], maxit, M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = cgs (Afun, b, [], maxit, M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = cgs (Afun, b, [], maxit, M1_fun, M2_fun);
-%! assert(flag, 0);
-
%!shared A, b, n, M
%!
%!test
@@ -382,10 +221,3 @@
%! [x, flag, relres, iter, resvec] = cgs (A, b, rtol, [], diag (diag (A)));
%! assert (x, ones (size (b)), 1e-7);
-%!test
-%! n = 5;
-%! A = rand(n);
-%! b = sum(A,2);
-%! M = ones(n);
-%! [x,flag] = cgs (A, b, [], [], M);
-%! assert (flag, 2)
diff --git a/scripts/sparse/gmres.m b/scripts/sparse/gmres.m
--- a/scripts/sparse/gmres.m
+++ b/scripts/sparse/gmres.m
@@ -1,5 +1,4 @@
-## Copyright (C) 2009-2015 Carlo de Falco
-## Copyright (C) 2016 Cristiano Dorigo
+## Copyright (C) 2009-2016 Carlo de Falco
##
## This file is part of Octave.
##
@@ -18,68 +17,36 @@
## .
## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{restart}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
-## @deftypefnx {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{restart}, @var{rtol}, @var{maxit}, @var{P}, [], @var{x0})
-## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@dots{})
+## @deftypefn {} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## @deftypefnx {} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{P})
+## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@dots{})
## Solve @code{A x = b} using the Preconditioned GMRES iterative method with
-## restart, a.k.a. PGMRES(restart).
+## restart, a.k.a. PGMRES(m).
##
-## The input arguments are:
## @itemize @minus
-##
-## @item @var{A} is the matrix of the linear system and it must be square.
-## @var{A} can be passed as a matrix, function handle, or inline
-## function @code{f} such that @code{f(x) = A * x}.
-##
-## @item @var{b} is the right hand side vector. It must be a column vector
-## with the same numbers of rows as @var{A}.
-##
-## @item @var{restart} is the number of iterations before that the
-## method restarts. If it is [] or N = numel (b), then the restart
-## is not applied.
-##
-## @item @var{RTOL} is the required relative tolerance for the
-## preconditioned residual error,
-## @code{inv (@var{M}) * (@var{b} - @var{a} * @var{x})}. The iteration stops if
-## @code{norm (inv (@var{M}) * (@var{b} - @var{a} * @var{x})) <=
-## @var{RTOL} * norm (inv (@var{M}) * @var{B})}. If @var{TOL} is omitted or
-## empty, then a tolerance of 1e-6 is used.
+## @item @var{rtol} is the relative tolerance,
+## if not given or set to [] the default value 1e-6 is used.
##
## @item @var{maxit} is the maximum number of outer iterations, if not given or
-## set to [], then the default value @code{min (10, @var{N} / @var{restart})}
-## is used.
-## Note that, if @var{restart} is empty, then @var{maxit} is the maximum number
-## of iterations. If @var{restart} and @var{maxit} are not empty, then
-## the maximum number of iterations is @code{@var{restart} * @var{maxit}}.
-## If both @var{restart} and @var{maxit} are empty, then the maximum
-## number of iterations is set to @code{min (10, @var{N})}.
-##
-## @item @var{M1}, @var{M2} are the preconditioners. The preconditioner
-## @var{P} is given as @code{P = M1 * M2}. Both @var{M1} and @var{M2} can
-## be passed as a matrix, function handle, or inline function @code{g} such
-## that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}. If @var{M1} is [] or not
-## given, then the preconditioner is not applied.
-## The technique used is the left-preconditioning, i.e., it is solved
-## @code{inv(@var{P}) * @var{A} * @var{x} = inv(@var{P}) * @var{b}} instead of
-## @code{@var{A} * @var{x} = @var{b}}.
+## set to [] the default value @code{min (10, numel (b) / restart)} is used.
##
## @item @var{x0} is the initial guess,
-## if not given or set to [], then the default value
-## @code{zeros (size (@var{b}))} is used.
+## if not given or set to [] the default value @code{zeros (size (b))} is used.
##
+## @item @var{m} is the restart parameter,
+## if not given or set to [] the default value @code{numel (b)} is used.
## @end itemize
##
-## The arguments which follow @var{x0} are treated as parameters, and passed in
-## a proper way to any of the functions (@var{A} or @var{P} or
-## @var{M1} or @var{M2}) which are passed to @code{gmres}.
+## Argument @var{A} can be passed as a matrix, function handle, or inline
+## function @code{f} such that @code{f(x) = A*x}.
##
-## The outputs are:
+## The preconditioner @var{P} is given as @code{P = M1 * M2}. Both @var{M1}
+## and @var{M2} can be passed as a matrix, function handle, or inline function
+## @code{g} such that @code{g(x) = M1\x} or @code{g(x) = M2\x}.
+##
+## Besides the vector @var{x}, additional outputs are:
##
## @itemize @minus
-##
-## @item @var{x} the computed approximation. If the method does not
-## converge, then it is the iterated with minimum residual.
-##
## @item @var{flag} indicates the exit status:
##
## @table @asis
@@ -87,463 +54,183 @@
##
## @item 1 : maximum number of iterations exceeded
##
-## @item 2 : the preconditioner matrix is singular
+## @item 2 : unused, but skipped for compatibility
##
-## @item 3 : algorithm reached stagnation (the relative difference between two
-## consecutive iterations is less than eps)
+## @item 3 : algorithm reached stagnation (no change between iterations)
## @end table
##
-## @item @var{relres} is the value of the relative preconditioned
-## residual of the approximation @var{x}.
+## @item @var{relres} is the final value of the relative residual.
##
## @item @var{iter} is a vector containing the number of outer iterations and
-## inner iterations performed to compute @var{x}. That is:
+## total iterations performed.
##
-## @itemize
-## @item @var{iter(1)}: number of outer iterations, i.e. how many
-## times the method restarted. (if @var{restart} is empty or @var{N},
-## then it is 1, if not 1 <= @var{iter(1)} <= @var{maxit}).
-##
-## @item @var{iter(2)}: the number of iterations performed before the
-## restart, i.e., the method restarts when
-## @code{@var{iter(2)} = @var{restart}}. If @var{restart} is empty or
-## @var{N}, then 1 <= @var{iter(2)} <= @var{maxit}.
+## @item @var{resvec} is a vector containing the relative residual at each
+## iteration.
## @end itemize
##
-## To be more clear, the approximation @var{x} is computed at the iteration
-## @code{(@var{iter(1)} - 1) * @var{restart} + @var{iter(2)}}.
-## Since the output @var{x} corresponds to the minimal preconditioned
-## residual solution, the total number of iterations that
-## the method performed is given by @code{length (resvec) - 1}.
-##
-## @item @var{resvec} is a vector containing the preconditioned
-## relative residual at each iteration, including the 0-th iteration
-## @code{norm (@var{A} * @var{x0} - @var{b})}.
-## @end itemize
-##
-## Let us consider a trivial problem with a tridiagonal matrix
-##
-## @example
-## @group
-## n = 20;
-## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-## toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-## sparse (1, 2, 1, 1, n) * n / 2);
-## b = A * ones (n, 1);
-## restart = 5;
-## [M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
-## M = M1 * M2;
-## Afun = @@(x) A * x;
-## Mfun = @@(x) M \ x;
-## M1fun = @@(x) M1 \ x;
-## M2fun = @@(x) M2 \ x;
-## @end group
-## @end example
-##
-## @sc{Example 1:} simplest usage of @code{gmres}
-##
-## @example
-## x = gmres (A, b, [], [], n)
-## @end example
-##
-## @sc{Example 2:} @code{gmres} with a function which computes
-## @code{@var{A} * @var{x}}
-##
-## @example
-## x = gmres (Afun, b, [], [], n)
-## @end example
-##
-## @sc{Example 3:} usage of @code{gmres} with the restart
-##
-## @example
-## x = gmres (A, b, restart);
-## @end example
-##
-## @sc{Example 4:} @code{gmres} with a preconditioner matrix @var{M}
-## with and without restart
-## @example
-## @group
-## x = gmres (A, b, [], 1e-06, n, M)
-## x = gmres (A, b, restart, 1e-06, n, M)
-## @end group
-## @end example
-##
-## @sc{Example 5:} @code{gmres} with a function as preconditioner
-##
-## @example
-## x = gmres (Afun, b, [], 1e-6, n, Mfun)
-## @end example
-##
-## @sc{Example 6:} @code{gmres} with preconditioner matrices @var{M1}
-## and @var{M2}
-##
-## @example
-## x = gmres (A, b, [], 1e-6, n, M1, M2)
-## @end example
-##
-## @sc{Example 7:} @code{gmres} with functions as preconditioners
-##
-## @example
-## x = gmres (Afun, b, 1e-6, n, M1fun, M2fun)
-## @end example
-##
-## @sc{Example 8:} @code {gmres} with as input a function requiring an argument
-##
-## @example
-## @group
-## function y = Ap (A, x, p) # compute A^p * x
-## y = x;
-## for i = 1:p
-## y = A * y;
-## endfor
-## endfunction
-## Apfun = @@(x, p) Ap (A, x, p);
-## x = gmres (Apfun, b, [], [], [], [], [], [], 2);
-## @end group
-## @end example
-##
-## @sc{Example 9:} explicit example to show that @code{gmres} uses a
-## left preconditioner
-##
-## @example
-## @group
-## [M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
-## M = M1 * M2;
-##
-## ## reference solution computed by gmres after two iterations
-## [x_ref, fl] = gmres (A, b, [], [], 1, M)
-##
-## ## left preconditioning
-## [x, fl] = gmres (M \ A, M \ b, [], [], 1)
-## x # compare x and x_ref
-##
-## @end group
-## @end example
-##
-## References:
-##
-## @enumerate
-##
-## @item @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear
-## Systems}, Second edition, 2003, SIAM
-##
-## @end enumerate
## @seealso{bicg, bicgstab, cgs, pcg, pcr, qmr}
## @end deftypefn
+function [x, flag, relres, it, resvec] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
-function [x_min, flag, relres, it, resvec] = ...
- gmres (A, b, restart = [], rtol = [], maxit = [], M1 = [],
- M2 = [], x0 = [], varargin)
-
- [A, M1, M2, is_prec, is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, M1, M2);
-
- ## Check if the inputs are empty, and in case set them
- [rtol, x0] = __default__input__ ({1e-06, zeros(size (b))}, rtol, x0);
-
- empty_restart = isempty (restart);
- empty_maxit = isempty (maxit);
- size_b = rows (b);
-
- if (rtol >= 1)
- warning("Input tol is bigger than 1. \n Try to use a smaller tolerance.");
- elseif (rtol <= eps / 2)
- warning("Input tol may not be achievable by gmres. \n Try to use a bigger tolerance.");
+ if (nargin < 2 || nargin > 8)
+ print_usage ();
endif
- ## This big "if block" is to set maxit and restart in the proper way
-
- if ((empty_restart) && (empty_maxit))
- restart = size_b;
- maxit = 1;
- max_iter_number = min (size_b, 10);
- elseif (restart <= 0) || (maxit <= 0)
- error ("gmres: MAXIT and RESTART must be positive integers")
- elseif (restart < size_b) && (empty_maxit)
- maxit = min (size_b / restart, 10);
- max_iter_number = maxit * restart;
- elseif (restart == size_b) && (empty_maxit)
- maxit = 1;
- max_iter_number = min (size_b, 10);
- elseif (restart > size_b) && (empty_maxit)
- warning ("RESTART is %d but it should be bounded by SIZE(A,2).\n Setting restart to %d. \n", restart, size_b)
- restart = size_b;
- maxit = 1;
- max_iter_number = restart;
- elseif (empty_restart) && (maxit <= size_b)
- restart = size_b;
- max_iter_number = maxit;
- elseif (empty_restart) && (maxit > size_b)
- warning ("MAXIT is %d but it should be bounded by SIZE(A,2). \n Setting MAXIT to %d", maxit, size_b);
- restart = size_b;
- maxit = size_b;
- max_iter_number = size_b;
- elseif (restart > size_b) && (!empty_maxit)
- warning ("RESTART is %d but it should be bounded by SIZE(A,2).\n Setting restart to %d. \n", restart, size_b)
- restart = size_b;
- max_iter_number = restart * maxit;
- elseif (restart == size_b) && (maxit <= size_b)
- max_iter_number = maxit;
+ if (ischar (A))
+ Ax = str2func (A);
+ elseif (isnumeric (A) && issquare (A))
+ Ax = @(x) A*x;
+ elseif (isa (A, "function_handle"))
+ Ax = A;
else
- max_iter_number = restart*maxit;
+ error ("gmres: A must be a function or square matrix");
endif
- prec_b_norm = norm (b, 2);
- if (prec_b_norm == 0)
- if (nargout < 2)
- printf("The right hand side vector is all zero so gmres\nreturned an all zero solution without iterating.\n")
- endif
- x_min = b;
- flag = 0;
- relres = 0;
- resvec = 0;
- it = [0, 0];
- return
+ if (nargin < 3 || isempty (restart))
+ restart = rows (b);
endif
- if (is_matrix_A && is_matrix_M1 && is_matrix_M2)
- [x_min, flag, it, iter, iter_min, resvec] = ...
- gmres_matrix (A, b, prec_b_norm, x0, restart, rtol, maxit, ...
- max_iter_number, is_prec, M1, M2);
- else
- if (is_matrix_A)
- A = @(x) A * x;
- endif
- if (is_prec)
- if (is_matrix_M1)
- M1 = @(x) M1 \ x;
- endif
- if (!isempty (M2) && is_matrix_M2)
- M2 = @(x) M2 \ x;
- endif
- endif
- [x_min, flag, it, iter, iter_min, resvec] = ...
- gmres_function (A, b, prec_b_norm, x0, restart, rtol, maxit, ...
- max_iter_number, is_prec, M1, M2, varargin{:});
+ if (nargin < 4 || isempty (rtol))
+ rtol = 1e-6;
endif
- if (flag == 2)
- resvec = norm (b);
- relres = 1;
- else
- resvec = resvec (1:iter);
- relres = resvec (iter) / prec_b_norm;
+ if (nargin < 5 || isempty (maxit))
+ maxit = min (rows (b)/restart, 10);
endif
- if ((relres <= rtol) && (flag == 1))
- flag = 0; # Converged to solution within tolerance
+ if (nargin < 6 || isempty (M1))
+ M1m1x = @(x) x;
+ elseif (ischar (M1))
+ M1m1x = str2func (M1);
+ elseif (isnumeric (M1) && ismatrix (M1))
+ M1m1x = @(x) M1 \ x;
+ elseif (isa (M1, "function_handle"))
+ M1m1x = M1;
+ else
+ error ("gmres: preconditioner M1 must be a function or matrix");
endif
- if ((nargout < 2) && (restart != size_b)) # restart applied
- switch (flag)
- case {0} # gmres converged
- printf ("gmres(%d) converged at outer iteration %d (inner iteration %d) ",restart, it (1), it (2));
- printf ("to a solution with relative residual %d \n", relres);
- case {1} # max number of iteration reached
- printf ("gmres(%d) stopped at outer iteration %d (inner iteration %d) ", restart, outer_it, restart_it-1);
- printf ("without converging to the desired tolerance %d ", rtol);
- printf ("because the maximum number of iterations was reached \n");
- printf ("The iterated returned (number %d(%d)) ", it(1), it(2));
- printf ("has relative residual %d \n", relres);
- case {2} # preconditioner singular
- printf ("gmres(%d) stopped at outer iteration %d (inner iteration %d) ",restart, outer_it, restart_it-1);
- printf ("without converging to the desired tolerance %d ", rtol);
- printf ("because the preconditioner matrix is singular \n");
- printf ("The iterated returned (number %d(%d)) ", it(1), it(2));
- printf ("has relative residual %d \n", relres);
- case {3} # stagnation
- printf ("gmres(%d) stopped at outer iteration %d (inner iteration %d) ", restart, outer_it, restart_it - 1);
- printf ("without converging to the desired tolerance %d", rtol);
- printf ("because it stagnates. \n");
- printf ("The iterated returned (number %d(%d)) ", it(1), it(2));
- printf ("has relative residual %d \n", relres);
- endswitch
- elseif ((nargout < 2) && (restart == size_b)) # no restart
- switch (flag)
- case {0} # gmres converged
- printf ("gmres converged at iteration %d ", it(2));
- printf ("to a solution with relative residual %d \n", relres);
- case {1} # max number of iteration reached
- printf ("gmres stopped at iteration %d ", restart_it - 1);
- printf ("without converging to the desired tolerance %d ", rtol);
- printf ("because the maximum number of iterations was reached \n");
- printf ("The iterated returned (number %d) ", it(2));
- printf ("has relative residual %d \n", relres);
- case {2} # preconditioner ill-conditioned
- printf ("gmres stopped at iteration %d ", restart_it - 1);
- printf ("without converging to the desired tolerance %d ", rtol);
- printf ("because the preconditioner matrix is singular \n")
- printf ("The iterated returned (number %d) ", it (2));
- printf ("has relative residual %d \n", relres);
- case {3} # stagnation
- printf ("gmres stopped at iteration %d ", restart_it - 1);
- printf ("without converging at the desired tolerance %d ", rtol);
- printf ("because it stagnates\n");
- printf ("The iterated returned (number %d) ", it(2));
- printf ("has relative residual %d \n", relres);
- endswitch
+ if (nargin < 7 || isempty (M2))
+ M2m1x = @(x) x;
+ elseif (ischar (M2))
+ M2m1x = str2func (M2);
+ elseif (isnumeric (M2) && ismatrix (M2))
+ M2m1x = @(x) M2 \ x;
+ elseif (isa (M2, "function_handle"))
+ M2m1x = M2;
+ else
+ error ("gmres: preconditioner M2 must be a function or matrix");
endif
+
+ Pm1x = @(x) M2m1x (M1m1x (x));
+
+ if (nargin < 8 || isempty (x0))
+ x0 = zeros (size (b));
+ endif
+
+ x_old = x0;
+ x = x_old;
+ prec_res = Pm1x (b - Ax (x_old));
+ presn = norm (prec_res, 2);
+
+ B = zeros (restart + 1, 1);
+ V = zeros (rows (x), restart);
+ H = zeros (restart + 1, restart);
+
+ ## begin loop
+ iter = 1;
+ restart_it = restart + 1;
+ resvec = zeros (maxit, 1);
+ resvec(1) = presn;
+ prec_b_norm = norm (Pm1x (b), 2);
+ flag = 1; # Default flag is maximum # of iterations exceeded
+
+ while (iter <= maxit * restart && presn > rtol * prec_b_norm)
+
+ ## restart
+ if (restart_it > restart)
+ restart_it = 1;
+ x_old = x;
+ prec_res = Pm1x (b - Ax (x_old));
+ presn = norm (prec_res, 2);
+ B(1) = presn;
+ H(:) = 0;
+ V(:, 1) = prec_res / presn;
+ endif
+
+ ## basic iteration
+ tmp = Pm1x (Ax (V(:, restart_it)));
+ [V(:,restart_it+1), H(1:restart_it+1, restart_it)] = ...
+ mgorth (tmp, V(:,1:restart_it));
+
+ Y = (H(1:restart_it+1, 1:restart_it) \ B(1:restart_it+1));
+
+ little_res = B(1:restart_it+1) - ...
+ H(1:restart_it+1, 1:restart_it) * Y(1:restart_it);
+
+ presn = norm (little_res, 2);
+
+ x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
+
+ resvec(iter+1) = presn;
+ if (norm (x - x_old, inf) <= eps)
+ flag = 3; # Stagnation: no change between iterations
+ break;
+ endif
+
+ restart_it += 1;
+ iter += 1;
+ endwhile
+
+ if (nargout > 1)
+ ## Calculate extra outputs as requested
+ relres = presn / prec_b_norm;
+ if (relres <= rtol)
+ flag = 0; # Converged to solution within tolerance
+ endif
+
+ it = [floor(iter/restart), restart_it-1];
+ endif
+
endfunction
+
%!demo
%! dim = 20;
%! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
%! b = ones (dim, 1);
-%! [x, flag, relres, iter, resvec] = ...
-%! gmres (A, b, 10, 1e-10, dim, @(x) x ./ diag (A), [], b)
+%! [x, flag, relres, iter, resvec] = gmres (A, b, 10, 1e-10, dim, @(x) x ./ diag (A), [], b)
-%!demo # simplest use
-%! n = 20;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! restart = 5;
-%! [M1, M2] = ilu (A + 0.1 * eye (n));
-%! M = M1 * M2;
-%! x = gmres (A, b, [], [], n);
-%! x = gmres (A, b, restart, [], n); # gmres with restart
-%! Afun = @(x) A * x;
-%! x = gmres (Afun, b, [], [], n);
-%! x = gmres (A, b,[], 1e-6, n, M); # gmres without restart
-%! x = gmres (A, b, [], 1e-6, n, M1, M2);
-%! Mfun = @(x) M \ x;
-%! x = gmres (Afun, b, [], 1e-6, n, Mfun);
-%! M1fun = @(x) M1 \ x;
-%! M2fun = @(x) M2 \ x;
-%! x = gmres (Afun, b, [], 1e-6, n, M1fun, M2fun);
-%! function y = Ap (A, x, p) # compute A^p * x
-%! y = x;
-%! for i = 1:p
-%! y = A * y;
-%! endfor
-%! endfunction
-%! Afun = @(x, p) Ap (A, x, p);
-%! x = gmres (Afun, b, [], [], n, [], [], [], 2); # solution of A^2 * x = b
-
-%!demo
-%! n = 10;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
-%! M = M1 * M2;
-%!
-%! ## reference solution computed by gmres after one iteration
-%! [x_ref, fl] = gmres (A, b, [], [], 1, M);
-%! x_ref
-%!
-%! ## left preconditioning
-%! [x, fl] = gmres ( M \ A, M \ b, [], [], 1);
-%! x # compare x and x_ref
-
+%!shared A, b, dim
+%! dim = 100;
%!test
-%! ## Check that all the subscripts works (matrix, funtion handle and mixed)
-%! A = magic (5);
-%! b = rand (5, 1);
-%! M1 = diag (sqrt (diag (A)));
-%! M2 = M1;
-%! Afun = @(z) A * z;
-%! M1_fun = @(z) M1 \ z;
-%! M2_fun = @(z) M2 \ z;
-%! [x, flag] = gmres (A, b);
-%! assert(flag, 0);
-%! [x, flag] = gmres (A, b, [],[],[], M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = gmres (A, b, [],[],[], M1_fun, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = gmres (A, b, [],[],[], M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = gmres (A, b, [],[],[], M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = gmres (Afun, b);
-%! assert(flag, 0);
-%! [x, flag] = gmres (Afun, b, [],[],[], M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = gmres (Afun, b, [],[],[], M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = gmres (Afun, b, [],[],[], M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = gmres (Afun, b, [],[],[], M1_fun, M2_fun);
-%! assert(flag, 0);
-
-%!test
-%! dim = 100;
%! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
%! b = ones (dim, 1);
-%! [x, flag] = gmres (A, b, 10, 1e-10, dim, @(x) x ./ diag (A), [], b);
+%! x = gmres (A, b, 10, 1e-10, dim, @(x) x ./ diag (A), [], b);
%! assert (x, A\b, 1e-9*norm (x, Inf));
-%! [x, flag] = gmres (A, b, dim, 1e-10, 1e4, @(x) diag (diag (A)) \ x, [], b);
+%!
+%!test
+%! x = gmres (A, b, dim, 1e-10, 1e4, @(x) diag (diag (A)) \ x, [], b);
%! assert(x, A\b, 1e-7*norm (x, Inf));
-
+%!
%!test
-%! dim = 100;
-%! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); ...
-%! [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim);
+%! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim);
%! A = A'*A;
%! b = rand (dim, 1);
%! [x, resvec] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag (A), [], []);
%! assert (x, A\b, 1e-9*norm (x, Inf));
-%! [x, flag] = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag (diag (A)) \ x, [], []);
+%! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag (diag (A)) \ x, [], []);
%! assert (x, A\b, 1e-9*norm (x, Inf));
-%! [x, flag] = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x ./ diag (A), [], []);
+%!test
+%! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x ./ diag (A), [], []);
%! assert (x, A\b, 1e-7*norm (x, Inf));
-%!test
-%! ## gmres solves complex linear systems
-%! N = 20;
-%! A = 2*rand(N)-1 + 1i * (2 * rand(N) - 1);
-%! x_exact = ones (N, 1);
-%! b = A * x_exact;
-%! [x, flag] = gmres(A, b, [], [], N);
-%! assert (flag, 0);
-%! assert (x, x_exact, -1e-06)
-%!test
-%! ## Maximum number of iteration reached
-%! A = rand(100);
-%! b = rand(100,1);
-%! [x, flag] = gmres (A, b);
-%! assert(flag, 1);
+%!error gmres (1)
+%!error gmres (1,2,3,4,5,6,7,8,9)
+%!error gmres ({1},2)
+%!error gmres ({1},2)
+%!error gmres (1,2,3,4,5,{6})
+%!error gmres (1,2,3,4,5,6,{7})
-%!test
-%! ## gmres recognizes that the preconditioner matrix is singular
-%! AA = rand (3);
-%! bb = ones (3, 1);
-%! I = eye (3);
-%! M = [1 0 0; 0 1 0; 0 0 0]; # the last row is zero
-%! [x, flag] = gmres(@(y) AA*y, bb, [], [], [], @(y) M\y, @(y) I\y);
-%! assert (flag, 2)
-
-%!test
-%! A = rand (4);
-%! A = A' * A;
-%! [x, flag] = gmres (A, zeros (4, 1), [], [], [], [], [], ones (4, 1));
-%! assert (x, zeros (4, 1))
-
-%!test
-%! A = rand (4);
-%! b = zeros (4, 1);
-%! [x, flag, relres, iter] = gmres (A, b);
-%! assert (relres, 0)
-
-%!test
-%! A = rand (4);
-%! b = A * ones(4,1);
-%! [x, flag, relres, iter] = gmres (A, b,[],[],[],[],[],...
-%! ones(4,1)+1e-8*(rand(4,1)*2-1));
-%! assert (iter, [0, 0])
-
-%!test
-%! A = rand(20);
-%! b = A * ones(20,1);
-%! [x,flag,relres,iter,resvec] = gmres(A,b,[],[],1);
-%! assert (iter, [1, 1])
-
-%!test
-%! A = rand(20);
-%! b = A * ones(20,1);
-%! [x,flag,relres,iter,resvec] = gmres (A,b,5);
-%! assert (iter, [4, 5])
diff --git a/scripts/sparse/ilu.m b/scripts/sparse/ilu.m
--- a/scripts/sparse/ilu.m
+++ b/scripts/sparse/ilu.m
@@ -167,11 +167,11 @@
endif
if (! (issparse (A) && issquare (A)))
- error ("ilu: A must be a sparse square matrix");
+ error ("ichol: A must be a sparse square matrix");
endif
if (! isstruct (opts))
- error ("ilu: OPTS must be a structure.");
+ error ("ichol: OPTS must be a structure.");
endif
## If A is empty then return empty L, U and P for Matlab compatibility
diff --git a/scripts/sparse/module.mk b/scripts/sparse/module.mk
--- a/scripts/sparse/module.mk
+++ b/scripts/sparse/module.mk
@@ -3,21 +3,7 @@
scripts/sparse/private
scripts_sparse_PRIVATE_FCN_FILES = \
- scripts/sparse/private/bicg_function.m \
- scripts/sparse/private/bicg_matrix.m \
- scripts/sparse/private/bicgstab_function.m \
- scripts/sparse/private/bicgstab_matrix.m \
- scripts/sparse/private/cgs_function.m \
- scripts/sparse/private/cgs_matrix.m \
- scripts/sparse/private/__sprand__.m \
- scripts/sparse/private/__checkandstring__.m \
- scripts/sparse/private/__default__input__.m \
- scripts/sparse/private/gmres_function.m \
- scripts/sparse/private/gmres_matrix.m \
- scripts/sparse/private/pcg_function.m \
- scripts/sparse/private/pcg_matrix.m \
- scripts/sparse/private/tfqmr_function.m \
- scripts/sparse/private/tfqmr_matrix.m
+ scripts/sparse/private/__sprand__.m
scripts_sparse_FCN_FILES = \
scripts/sparse/bicg.m \
@@ -46,7 +32,6 @@
scripts/sparse/spstats.m \
scripts/sparse/spy.m \
scripts/sparse/svds.m \
- scripts/sparse/tfqmr.m \
scripts/sparse/treelayout.m \
scripts/sparse/treeplot.m
diff --git a/scripts/sparse/pcg.m b/scripts/sparse/pcg.m
--- a/scripts/sparse/pcg.m
+++ b/scripts/sparse/pcg.m
@@ -1,5 +1,4 @@
-## Copyright (C) 2004-2015 Piotr Krzyzanowski
-## Copyright (C) 2016 Cristiano Dorigo
+## Copyright (C) 2004-2016 Piotr Krzyzanowski
##
## This file is part of Octave.
##
@@ -18,22 +17,22 @@
## .
## -*- texinfo -*-
-## @deftypefn {Function File} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{})
-## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{})
+## @deftypefn {} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{})
+## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{})
##
## Solve the linear system of equations @w{@code{@var{A} * @var{x} = @var{b}}}
## by means of the Preconditioned Conjugate Gradient iterative method.
##
-## The input arguments are:
+## The input arguments are
##
## @itemize
## @item
## @var{A} can be either a square (preferably sparse) matrix or a function
## handle, inline function or string containing the name of a function which
-## computes @w{@code{@var{A} * @var{x}}}. @var{A} has to be
-## Hermitian and Positive Definite (HPD); if @code{pcg} finds @var{A} not to be
-## positive definite, a warning is printed and the @var{flag} output
-## is set.
+## computes @w{@code{@var{A} * @var{x}}}. In principle, @var{A} should be
+## symmetric and positive definite; if @code{pcg} finds @var{A} not to be
+## positive definite, a warning is printed and the @var{flag} output will be
+## set.
##
## @item
## @var{b} is the right-hand side vector.
@@ -42,36 +41,25 @@
## @var{tol} is the required relative tolerance for the residual error,
## @w{@code{@var{b} - @var{A} * @var{x}}}. The iteration stops if
## @w{@code{norm (@var{b} - @var{A} * @var{x})} @leq{}
-## @w{@code{@var{tol} * norm (@var{b})}}}.
-## If @var{tol} is omitted or empty, then a tolerance of 1e-6 is used.
+## @w{@var{tol} * norm (@var{b})}}.
+## If @var{tol} is omitted or empty then a tolerance of 1e-6 is used.
##
## @item
-## @var{maxit} is the maximum allowed number of iterations; if @var{maxit}
+## @var{maxit} is the maximum allowable number of iterations; if @var{maxit}
## is omitted or empty then a value of 20 is used.
##
## @item
-## @var{m} is a HPD preconditioning matrix. For any decomposition
-## @code{@var{m} = @var{p1} * @var {p2}} such that
-## @w{@code{inv (@var{p1}) * @var{A} * inv (@var{p2})}} is HPD, the
-## conjugate gradient method is formally applied to the linear system
-## @w{@code{inv (@var{p1}) * @var{A} * inv (@var{p2}) * @var{y} = inv
-## (@var{p1}) * @var{b}}},
-## with @code{@var{x} = inv (@var{p2}) * @var{y}} (split preconditioning).
-## In practice, at each iteration of the conjugate gradient method a
-## linear system with matrix @var{m} is solved with @code{mldivide}.
-## If a particular factorization
-## @code{@var{m} = @var{m1} * @var{m2}} is available (for instance, an
-## incomplete Cholesky factorization of @var{a}), the two matrices
-## @var{m1} and @var{m2} can be passed and the relative linear systems
-## are solved with the @code{mldivide} operator.
+## @var{m} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that
+## the iteration is (theoretically) equivalent to solving by @code{pcg}
+## @w{@code{@var{P} * @var{x} = @var{m} \ @var{b}}}, with
+## @w{@code{@var{P} = @var{m} \ @var{A}}}.
## Note that a proper choice of the preconditioner may dramatically improve
## the overall performance of the method. Instead of matrices @var{m1} and
## @var{m2}, the user may pass two functions which return the results of
-## applying the inverse of @var{m1} and @var{m2} to a vector.
-## If @var{m1} is omitted or empty @code{[]}, then no preconditioning
-## is applied. If no factorization of @var{m} is available, @var{m2}
-## can be omitted or left [], and the input variable @var{m1} can be
-## used to pass the preconditioner @var{m}.
+## applying the inverse of @var{m1} and @var{m2} to a vector (usually this is
+## the preferred way of using the preconditioner). If @var{m1} is omitted or
+## empty @code{[]} then no preconditioning is applied. If @var{m2} is
+## omitted, @var{m} = @var{m1} will be used as a preconditioner.
##
## @item
## @var{x0} is the initial guess. If @var{x0} is omitted or empty then the
@@ -79,54 +67,39 @@
## @end itemize
##
## The arguments which follow @var{x0} are treated as parameters, and passed in
-## a proper way to any of the functions (@var{A} or @var{m1} or
-## @var{m2}) which are passed to @code{pcg}.
-## See the examples below for further details.
-##
-## The output arguments are:
+## a proper way to any of the functions (@var{A} or @var{m}) which are passed
+## to @code{pcg}. See the examples below for further details. The output
+## arguments are
##
## @itemize
## @item
## @var{x} is the computed approximation to the solution of
-## @w{@code{@var{A} * @var{x} = @var{b}}}. If the algorithm did not converge,
-## then @var{x} is the iterated which has the minimum residual.
+## @w{@code{@var{A} * @var{x} = @var{b}}}.
##
## @item
-## @var{flag} reports on the convergence:
-## @itemize
-## @item 0: The algorithm converged at the prescribed tolerance.
-## @item 1: The algorithm did not converge and it reached the maximum
-## number of iterations.
-## @item 2: The preconditioner matrix is singular.
-## @item 3: The algorithm stagnated, i.e. the absolute value of the
-## difference between
-## the actual iteration @var{x} and the previous is less than
-## @code{@var{eps} * norm (@var{x},2)}.
-## @item 4: The algorithm detects that the input (preconditioned) matrix is not
-## HPD.
-## @end itemize
+## @var{flag} reports on the convergence. A value of 0 means the solution
+## converged and the tolerance criterion given by @var{tol} is satisfied.
+## A value of 1 means that the @var{maxit} limit for the iteration count was
+## reached. A value of 3 indicates that the (preconditioned) matrix was found
+## not to be positive definite.
##
## @item
## @var{relres} is the ratio of the final residual to its initial value,
## measured in the Euclidean norm.
##
## @item
-## @var{iter} indicates the iteration of @var{x} which it was
-## computed. Since the output @var{x} corresponds to the minimal
-## residual solution, the total number of iterations that
-## the method performed is given by @code{length(resvec) - 1}.
+## @var{iter} is the actual number of iterations performed.
##
## @item
## @var{resvec} describes the convergence history of the method.
-## @code{@var{resvec} (@var{i}, 1)} is the Euclidean norm of the residual, and
-## @code{@var{resvec} (@var{i}, 2)} is the preconditioned residual
-## norm, after the
+## @code{@var{resvec}(i,1)} is the Euclidean norm of the residual, and
+## @code{@var{resvec}(i,2)} is the preconditioned residual norm, after the
## (@var{i}-1)-th iteration, @code{@var{i} = 1, 2, @dots{}, @var{iter}+1}.
## The preconditioned residual norm is defined as
-## @code{@var{r}' * (@var{m} \ @var{r})} where
+## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where
## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the
## description of @var{m}. If @var{eigest} is not required, only
-## @code{@var{resvec} (:, 1)} is returned.
+## @code{@var{resvec}(:,1)} is returned.
##
## @item
## @var{eigest} returns the estimate for the smallest @code{@var{eigest}(1)}
@@ -138,23 +111,20 @@
## @code{@var{eigest}(2) / @var{eigest}(1)} is a lower bound for
## @code{cond (@var{P}, 2)}, which nevertheless in the limit should
## theoretically be equal to the actual value of the condition number.
+## The method which computes @var{eigest} works only for symmetric positive
+## definite @var{A} and @var{m}, and the user is responsible for verifying this
+## assumption.
## @end itemize
##
-##
-## Let us consider a trivial problem with a tridiagonal matrix
+## Let us consider a trivial problem with a diagonal matrix (we exploit the
+## sparsity of A)
##
## @example
## @group
## n = 10;
-## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
-## b = A * ones (n, 1);
-## M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
-## M2 = M1';
-## M = M1 * M2;
-## Afun = @@(x) A * x;
-## Mfun = @@(x) M \ x;
-## M1fun = @@(x) M1 \ x;
-## M2fun = @@(x) M2 \ x;
+## A = diag (sparse (1:n));
+## b = rand (n, 1);
+## [l, u, p] = ilu (A, struct ("droptol", 1.e-3));
## @end group
## @end example
##
@@ -168,65 +138,60 @@
## @code{@var{A} * @var{x}}
##
## @example
-## x = pcg (Afun, b)
+## @group
+## function y = apply_a (x)
+## y = [1:N]' .* x;
+## endfunction
+##
+## x = pcg ("apply_a", b)
+## @end group
## @end example
##
-## @sc{Example 3:} @code{pcg} with a preconditioner matrix @var{M}
+## @sc{Example 3:} @code{pcg} with a preconditioner: @var{l} * @var{u}
##
## @example
-## x = pcg (A, b, 1e-06, 100, M)
+## x = pcg (A, b, 1.e-6, 500, l*u)
## @end example
##
-## @sc{Example 4:} @code{pcg} with a function as preconditioner
+## @sc{Example 4:} @code{pcg} with a preconditioner: @var{l} * @var{u}.
+## Faster than @sc{Example 3} since lower and upper triangular matrices are
+## easier to invert
##
## @example
-## x = pcg (Afun, b, 1e-6, 100, Mfun)
+## x = pcg (A, b, 1.e-6, 500, l, u)
## @end example
##
-## @sc{Example 5:} @code{pcg} with preconditioner matrices @var{M1}
-## and @var{M2}
-##
-## @example
-## x = pcg (A, b, 1e-6, 100, M1, M2)
-## @end example
-##
-## @sc{Example 6:} @code{pcg} with functions as preconditioners
-##
-## @example
-## x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun)
-## @end example
-##
-## @sc{Example 7:} @code{pcg} with as input a function requiring an argument
+## @sc{Example 5:} Preconditioned iteration, with full diagnostics. The
+## preconditioner (quite strange, because even the original matrix @var{A} is
+## trivial) is defined as a function
##
## @example
## @group
-## function y = Ap (A, x, p) # compute A^p * x
-## y = x;
-## for i = 1:p
-## y = A * y;
-## endfor
-## endfunction
-## Apfun = @@(x, p) Ap (A, x, p);
-## x = pcg (Apfun, b, [], [], [], [], [], 2);
+## function y = apply_m (x)
+## k = floor (length (x) - 2);
+## y = x;
+## y(1:k) = x(1:k) ./ [1:k]';
+## endfunction
+##
+## [x, flag, relres, iter, resvec, eigest] = ...
+## pcg (A, b, [], [], "apply_m");
+## semilogy (1:iter+1, resvec);
## @end group
## @end example
##
-## @sc{Example 8:} explicit example to show that @code{pcg} uses a
-## split preconditioner
+## @sc{Example 6:} Finally, a preconditioner which depends on a parameter
+## @var{k}.
##
## @example
## @group
-## M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
-## M2 = M1';
-## M = M1 * M2;
+## function y = apply_M (x, varargin)
+## K = varargin@{1@};
+## y = x;
+## y(1:K) = x(1:K) ./ [1:K]';
+## endfunction
##
-## ## reference solution computed by pcg after two iterations
-## [x_ref, fl] = pcg (A, b, [], 2, M)
-##
-## ## split preconditioning
-## [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
-## x = M2 \ y # compare x and x_ref
-##
+## [x, flag, relres, iter, resvec, eigest] = ...
+## pcg (A, b, [], [], "apply_m", [], [], 3)
## @end group
## @end example
##
@@ -244,157 +209,279 @@
## @url{http://www-users.cs.umn.edu/~saad/books.html}
## @end enumerate
##
-## @seealso{sparse, pcr, gmres, bicg, bicgstab, cgs}
+## @seealso{sparse, pcr}
## @end deftypefn
## Author: Piotr Krzyzanowski
## Modified by: Vittoria Rezzonico
## - Add the ability to provide the pre-conditioner as two separate matrices
+function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, tol, maxit, m1, m2, x0, varargin)
-function [x_min, flag, relres, iter_min, resvec, eigest] =...
- pcg (A, b, tol = [], maxit = [], m1 = [], m2 = [], x0 = [], varargin)
+ ## M = M1*M2
- ## Insert the default input (if necessary)
- [tol, maxit, x0] = __default__input__ ({1e-6, min(rows (b), 20),...
- zeros(size (b))}, tol, maxit, x0);
-
- if (tol >= 1)
- warning ("Input tol is bigger than 1. \n Try to use a smaller tolerance.");
- elseif (tol <= eps / 2)
- warning ("Input tol may not be achievable by pcg. \n Try to use a bigger tolerance");
+ if (nargin < 7 || isempty (x0))
+ x = zeros (size (b));
+ else
+ x = x0;
endif
- ## Check if the input data A,b,m1,m2 are consistent (i.e. if they are
- ## matrix or function handle)
+ if (nargin < 5 || isempty (m1))
+ exist_m1 = 0;
+ else
+ exist_m1 = 1;
+ endif
- [A, m1, m2, is_prec, is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, m1, m2);
+ if (nargin < 6 || isempty (m2))
+ exist_m2 = 0;
+ else
+ exist_m2 = 1;
+ endif
+
+ if (nargin < 4 || isempty (maxit))
+ maxit = min (rows (b), 20);
+ endif
maxit += 2;
- n_arg_out = nargout;
- ## Set Initial data
- b_norm = norm (b);
- if (b_norm == 0)
- if (n_arg_out < 2)
- printf("The right hand side vector is all zero so pcg \n");
- printf ("returned an all zero solution without iterating.\n");
- endif
- x_min = b;
- flag = 0;
- relres = 0;
- resvec = 0;
- iter_min = 0;
- eigest = [NaN, NaN];
- return
+ if (nargin < 3 || isempty (tol))
+ tol = 1e-6;
endif
- if (is_matrix_A && is_matrix_M1 && is_matrix_M2)
- [x_min, flag, iter, iter_min, resvec, eigest] = ...
- pcg_matrix (A, b, b_norm, x0, tol, maxit, is_prec, m1, m2, n_arg_out);
- else
- if (is_matrix_A)
- A = @(x) A * x;
- endif
- if (is_prec)
- if (is_matrix_M1)
- m1 = @(x) m1 \ x;
- endif
- if (!isempty (m2) && is_matrix_M2)
- m2 = @(x) m2 \ x;
- endif
- endif
- [x_min, flag, iter, iter_min, resvec, eigest] = ...
- pcg_function (A, b, b_norm, x0, tol, maxit, is_prec, m1, m2, ...
- n_arg_out, varargin{:});
+ preconditioned_residual_out = false;
+ if (nargout > 5)
+ T = zeros (maxit, maxit);
+ preconditioned_residual_out = true;
endif
- ## Set the last variables
+ ## Assume A is positive definite.
+ matrix_positive_definite = true;
- if (flag == 2)
- relres = 1;
- elseif (resvec (1, 1) == 0)
- relres = 0;
+ p = zeros (size (b));
+ oldtau = 1;
+ if (isnumeric (A))
+ ## A is a matrix.
+ r = b - A*x;
else
- relres = resvec(iter_min+1, 1) ./ resvec(1, 1);
+ ## A should be a function.
+ r = b - feval (A, x, varargin{:});
endif
- iter -= 2; # compatibility
+ b_norm = norm (b);
+ resvec(1,1) = norm (r);
+ alpha = 1;
+ iter = 2;
- ## Set the flag in the proper way if flag not 3, 4 or 2
- if (flag == 2)
- flag = 2;
- elseif (flag == 1) && (relres <= tol)
- flag = 0;
+ while (resvec(iter-1,1) > tol * b_norm && iter < maxit)
+ if (exist_m1)
+ if (isnumeric (m1))
+ y = m1 \ r;
+ else
+ y = feval (m1, r, varargin{:});
+ endif
+ else
+ y = r;
+ endif
+ if (exist_m2)
+ if (isnumeric (m2))
+ z = m2 \ y;
+ else
+ z = feval (m2, y, varargin{:});
+ endif
+ else
+ z = y;
+ endif
+ tau = z' * r;
+ resvec(iter-1,2) = sqrt (tau);
+ beta = tau / oldtau;
+ oldtau = tau;
+ p = z + beta * p;
+ if (isnumeric (A))
+ ## A is a matrix.
+ w = A * p;
+ else
+ ## A should be a function.
+ w = feval (A, p, varargin{:});
+ endif
+ ## Needed only for eigest.
+ oldalpha = alpha;
+ alpha = tau / (p'*w);
+ if (alpha <= 0.0)
+ ## Negative matrix.
+ matrix_positive_definite = false;
+ endif
+ x += alpha * p;
+ r -= alpha * w;
+ if (nargout > 5 && iter > 2)
+ T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
+ [1 sqrt(beta); sqrt(beta) beta]./oldalpha;
+ ## EVS = eig (T(2:iter-1,2:iter-1));
+ ## fprintf (stderr,"PCG condest: %g (iteration: %d)\n", max (EVS)/min (EVS),iter);
+ endif
+ resvec(iter,1) = norm (r);
+ iter += 1;
+ endwhile
+
+ if (nargout > 5)
+ if (matrix_positive_definite)
+ if (iter > 3)
+ T = T(2:iter-2,2:iter-2);
+ l = eig (T);
+ eigest = [min(l), max(l)];
+ ## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1));
+ else
+ eigest = [NaN, NaN];
+ warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
+ endif
+ else
+ eigest = [NaN, NaN];
+ endif
+
+ ## Apply the preconditioner once more and finish with the precond
+ ## residual.
+ if (exist_m1)
+ if (isnumeric (m1))
+ y = m1 \ r;
+ else
+ y = feval (m1, r, varargin{:});
+ endif
+ else
+ y = r;
+ endif
+ if (exist_m2)
+ if (isnumeric (m2))
+ z = m2 \ y;
+ else
+ z = feval (m2, y, varargin{:});
+ endif
+ else
+ z = y;
+ endif
+
+ resvec(iter-1,2) = sqrt (r' * z);
+ else
+ resvec = resvec(:,1);
endif
- if (n_arg_out < 2)
- switch (flag)
- case {0}
- printf ("pcg converged at iteration %d ", iter_min);
- printf ("with relative residual %d\n", relres);
- case {1}
- printf ("pcg stopped at iteration %d ", iter+1);
- printf ("without converging to the desired tolerance %d ", tol);
- printf ("because the maximum number of iteration was reached, \n");
- printf ("The iterated returned (number %d) ",iter_min);
- printf ("has relative residual %d \n", relres);
- case {2}
- printf ("pcg stopped at iteration %d ", iter+1)
- printf ("without converging to the desired tolerance %d ", tol);
- printf ("because the preconditioned matrix is singular.\n");
- printf ("The iterated returned (number %d) ", iter_min);
- printf ("has relative residual %d \n", relres);
- case {3}
- printf ("pcg stopped at iteration %d ", iter+1);
- printf ("without converging to the desired tolerance %d ", tol);
- printf ("because of stagnation. \n");
- printf ("The iterated returned (number %d) ", iter_min);
- printf ("has relative residual %d.\n", relres);
- case {4}
- printf ("pcg stopped at iteration %d ", iter + 1);
- printf ("without converging to the desired tolerance %d ",tol);
- printf ("because the (preconditioned) matrix is not positive definite. \n");
- printf ("The iterate returned (number %d) ", iter_min);
- printf ("has relative residual %d \n", relres);
- endswitch
+ flag = 0;
+ relres = resvec(iter-1,1) ./ resvec(1,1);
+ iter -= 2;
+ if (iter >= maxit - 2)
+ flag = 1;
+ if (nargout < 2)
+ warning ("pcg: maximum number of iterations (%d) reached\n", iter);
+ warning ("pcg: the initial residual norm was reduced %g times.\n",
+ 1.0 / relres);
+ endif
+ elseif (nargout < 2)
+ fprintf (stderr, "pcg: converged in %d iterations. ", iter);
+ fprintf (stderr, "pcg: the initial residual norm was reduced %g times.\n",
+ 1.0/relres);
endif
+
+ if (! matrix_positive_definite)
+ flag = 3;
+ if (nargout < 2)
+ warning ("pcg: matrix not positive definite?\n");
+ endif
+ endif
+
endfunction
-%!test
-%! ## Check that all the subscripts works
-%! A = toeplitz (sparse ([2, 1 ,0, 0, 0]));
-%! b = A * ones (5, 1);
-%! M1 = diag (sqrt (diag (A)));
-%! M2 = M1; # M1 * M2 is the Jacobi preconditioner
-%! Afun = @(z) A*z;
-%! M1_fun = @(z) M1 \ z;
-%! M2_fun = @(z) M2 \ z;
-%! [x, flag, ~, iter] = pcg (A,b);
-%! assert(flag, 0);
-%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2);
-%! assert(flag, 0);
-%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = pcg (A, b, [], [], M1_fun, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = pcg (A, b,[],[], M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = pcg (A, b,[],[], M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = pcg (Afun, b);
-%! assert(flag, 0);
-%! [x, flag] = pcg (Afun, b,[],[], M1 * M2);
-%! assert(flag, 0);
-%! [x, flag] = pcg (Afun, b,[],[], M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = pcg (Afun, b,[],[], M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2_fun);
-%! assert(flag, 0);
+
+%!demo
+%! ## Simplest usage of pcg (see also 'help pcg')
+%!
+%! N = 10;
+%! A = diag ([1:N]); b = rand (N, 1);
+%! y = A \ b; # y is the true solution
+%! x = pcg (A, b);
+%! printf ("The solution relative error is %g\n", norm (x - y) / norm (y));
+%!
+%! ## You shouldn't be afraid if pcg issues some warning messages in this
+%! ## example: watch out in the second example, why it takes N iterations
+%! ## of pcg to converge to (a very accurate, by the way) solution
+
+%!demo
+%! ## Full output from pcg, except for the eigenvalue estimates
+%! ## We use this output to plot the convergence history
+%!
+%! N = 10;
+%! A = diag ([1:N]); b = rand (N, 1);
+%! X = A \ b; # X is the true solution
+%! [x, flag, relres, iter, resvec] = pcg (A, b);
+%! printf ("The solution relative error is %g\n", norm (x - X) / norm (X));
+%! title ("Convergence history");
+%! semilogy ([0:iter], resvec / resvec(1), "o-g");
+%! xlabel ("Iteration"); ylabel ("log(||b-Ax||/||b||)");
+%! legend ("relative residual");
+
+%!demo
+%! ## Full output from pcg, including the eigenvalue estimates
+%! ## Hilbert matrix is extremely ill-conditioned, so pcg WILL have problems
+%!
+%! N = 10;
+%! A = hilb (N); b = rand (N, 1);
+%! X = A \ b; # X is the true solution
+%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200);
+%! printf ("The solution relative error is %g\n", norm (x - X) / norm (X));
+%! printf ("Condition number estimate is %g\n", eigest(2) / eigest(1));
+%! printf ("Actual condition number is %g\n", cond (A));
+%! title ("Convergence history");
+%! semilogy ([0:iter], resvec, ["o-g";"+-r"]);
+%! xlabel ("Iteration"); ylabel ("log(||b-Ax||)");
+%! legend ("absolute residual", "absolute preconditioned residual");
+
+%!demo
+%! ## Full output from pcg, including the eigenvalue estimates
+%! ## We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2)
+%! ## and that's the reason we need some preconditioner; here we take
+%! ## a very simple and not powerful Jacobi preconditioner,
+%! ## which is the diagonal of A.
+%!
+%! N = 100;
+%! A = zeros (N, N);
+%! for i = 1 : N - 1 # form 1-D Laplacian matrix
+%! A(i:i+1, i:i+1) = [2 -1; -1 2];
+%! endfor
+%! b = rand (N, 1);
+%! X = A \ b; # X is the true solution
+%! maxit = 80;
+%! printf ("System condition number is %g\n", cond (A));
+%! ## No preconditioner: the convergence is very slow!
+%!
+%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit);
+%! printf ("System condition number estimate is %g\n", eigest(2) / eigest(1));
+%! title ("Convergence history");
+%! semilogy ([0:iter], resvec(:,1), "o-g");
+%! xlabel ("Iteration"); ylabel ("log(||b-Ax||)");
+%! legend ("NO preconditioning: absolute residual");
+%!
+%! pause (1);
+%! ## Test Jacobi preconditioner: it will not help much!!!
+%!
+%! M = diag (diag (A)); # Jacobi preconditioner
+%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
+%! printf ("JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1));
+%! hold on;
+%! semilogy ([0:iter], resvec(:,1), "o-r");
+%! legend ("NO preconditioning: absolute residual", ...
+%! "JACOBI preconditioner: absolute residual");
+%!
+%! pause (1);
+%! ## Test nonoverlapping block Jacobi preconditioner: it will help much!
+%!
+%! M = zeros (N, N); k = 4;
+%! for i = 1 : k : N # form 1-D Laplacian matrix
+%! M(i:i+k-1, i:i+k-1) = A(i:i+k-1, i:i+k-1);
+%! endfor
+%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
+%! printf ("BLOCK JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1));
+%! semilogy ([0:iter], resvec(:,1), "o-b");
+%! legend ("NO preconditioning: absolute residual", ...
+%! "JACOBI preconditioner: absolute residual", ...
+%! "BLOCK JACOBI preconditioner: absolute residual");
+%! hold off;
%!test
%! ## solve small diagonal system
@@ -407,15 +494,16 @@
%! assert (flag, 0);
%!test
-%! ## A not positive definite
+%! ## solve small indefinite diagonal system
+%! ## Despite A being indefinite, the iteration continues and converges.
%! ## The indefiniteness of A is detected.
%!
%! N = 10;
-%! A = -diag([1:N]); b = rand (N, 1);
+%! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1);
%! X = A \ b; # X is the true solution
%! [x, flag] = pcg (A, b, [], N+1);
-%! assert (flag, 4)
-
+%! assert (norm (x - X) / norm (X), 0, 1e-10);
+%! assert (flag, 3);
%!test
%! ## solve tridiagonal system, do not converge in default 20 iterations
@@ -429,9 +517,8 @@
%! X = A \ b; # X is the true solution
%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12);
%! assert (flag);
-%! assert (relres >= 1.0);
-%! #assert (iter, 20); # should perform max allowable
-%! # default number of iterations
+%! assert (relres > 1.0);
+%! assert (iter, 20); # should perform max allowable default number of iterations
%!warning
%! ## solve tridiagonal system with "perfect" preconditioner which converges
@@ -447,109 +534,6 @@
%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
%! assert (norm (x - X) / norm (X), 0, 1e-6);
%! assert (flag, 0);
-%! assert (iter, 1); # should converge in one iteration
+%! assert (iter, 1); # should converge in one iteration
%! assert (isnan (eigest), isnan ([NaN, NaN]));
-%!test
-%! ## pcg detect a non-Hermitian matrix, with a considerable imaginary part
-%! ## With this example, Matlab doesn't recognize the wrong type of matrix and
-%! ## makes iterations until it reaches the maxii
-%! N = 10;
-%! A = diag (1:N) + 1i*1e-04*rand (N);
-%! b = ones (N, 1);
-%! [x,flag] = pcg (A, b, []);
-%! assert (flag, 4)
-
-%!test
-%! ## The imaginary part is not influent (it is too small), so pcg doesn't stop
-%! N = 10;
-%! A = diag (1:N) + 1i*1e-10*rand (N);
-%! b = ones (N, 1);
-%! [x,flag] = pcg (A, b, [], N+1);
-%! assert (flag, 0)
-%! assert (x, A\b, -1e-6)
-
-%!test
-%! ## pcg solves linear system with A Hermitian positive definite
-%! N = 20;
-%! A = 2*rand (N)-1 + 1i*(2*rand (N)-1);
-%! A = A'*A;
-%! b = A * ones (N,1);
-%! Hermitian_A = ishermitian (A);
-%! [x,flag] = pcg (A, b, [], 2*N);
-%! assert (Hermitian_A, true)
-%! assert (flag, 0)
-%! assert (x, ones (N, 1), -1e-4)
-
-%!test
-%! ## pcg solves preconditioned linear system with A HPD
-%! N = 20;
-%! A = 2*rand (N)-1 + 1i*(2*rand (N)-1);
-%! A = A' * A;
-%! b = A * ones (N,1);
-%! M2 = chol (A + 0.1 * eye (N)); # factor of a perturbed matrix
-%! M = M2' * M2;
-%! Hermitian_A = ishermitian (A);
-%! Hermitian_M = ishermitian (M);
-%! [x,flag] = pcg (A, b, [], 2*N, M);
-%! assert (Hermitian_A, true)
-%! assert (Hermitian_M, true)
-%! assert (flag, 0)
-%! assert (x, ones (N, 1), -1e-4)
-
-%!test
-%! ## pcg recognizes that the preconditioner matrix is singular
-%! N = 3;
-%! A = rand(3);
-%! A = A*A';
-%! M = [1 0 0; 0 1 0; 0 0 0]; # the last rows is zero
-%! [x,flag] = pcg (A, ones(3,1), [], [], M);
-%! assert (flag, 2)
-
-%!test
-%! A = rand (4);
-%! A = A' * A;
-%! [x, flag] = pcg (A, zeros (4, 1), [], [], [], [], ones (4, 1));
-%! assert (x, zeros (4, 1))
-
-%!demo # simplest use
-%! n = 10;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
-%! b = A * ones (n, 1);
-%! M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
-%! M2 = M1';
-%! M = M1 * M2;
-%! x = pcg (A, b);
-%! Afun = @(x) A * x;
-%! x = pcg (Afun, b);
-%! x = pcg (A, b, 1e-6, 100, M);
-%! x = pcg (A, b, 1e-6, 100, M1, M2);
-%! Mfun = @(x) M \ x;
-%! x = pcg (Afun, b, 1e-6, 100, Mfun);
-%! M1fun = @(x) M1 \ x;
-%! M2fun = @(x) M2 \ x;
-%! x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun);
-%! function y = Ap (A, x, p) # compute A^p * x
-%! y = x;
-%! for i = 1:p
-%! y = A * y;
-%! endfor
-%! endfunction
-%! Afun = @(x, p) Ap (A, x, p);
-%! x = pcg (Afun, b, [], [], [], [], [], 2); # solution of A^2 * x = b
-
-%!demo
-%! n = 10;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
-%! b = A * ones (n, 1);
-%! M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
-%! M2 = M1';
-%! M = M1 * M2;
-%!
-%! ## reference solution computed by pcg after two iterations
-%! [x_ref, fl] = pcg (A, b, [], 2, M);
-%! x_ref
-%!
-%! ## split preconditioning
-%! [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2);
-%! x = M2 \ y # compare x and x_ref
diff --git a/scripts/sparse/private/__checkandstring__.m b/scripts/sparse/private/__checkandstring__.m
deleted file mode 100644
--- a/scripts/sparse/private/__checkandstring__.m
+++ /dev/null
@@ -1,105 +0,0 @@
-## Copyright (C) 2016 Cristiano Dorigo, Octave Arena
-##
-## 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} {[...] =} __checkandstring__ (@var{A}, @var{b}, @var{M1}, @var{M2})
-## @deftypefnx {Function File} {[@var{A}, @var{M1}, @var{M2}, @var{is_prec}, @var{type}] =} __checkandstring (...)
-##
-## Check if the parameters @var{A} (matrix of our linear system), @var{b}
-## (right hand side vector), @var{M1}, @var{M2} (preconditioner matrices) are
-## effectively matrices or functions handle, summarizing if they are void or not.
-##
-## The input parameters are:
-##
-## @itemize
-## @item @var{A} is the matrix of the linear system.
-##
-## @item @var{b} is the right hand side vector.
-##
-## @item @var{M1}, @var{M2} preconditioners. They can be [].
-##
-## @end itemize
-##
-## The output parameters are:
-## @itemize
-##
-## @item @var{A}, @var{M1}, @var{M2} can be strings in input, so they returned as the corresponding function of these strings
-##
-## @item @var{is_prec} boolean telling if a preconditioner has to be applied
-##
-## @item @var{is_matrix_A}, @var{is_matrix_M1}, and
-## @var{is_matrix_M2}, booleans telling if the corresponding input is
-## a square matrix.
-## @end itemize
-##
-## @end itemize
-## @end deftypefn
-
-function [A, M1, M2, is_prec, is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, M1, M2)
-
-## Check A and set its type
- if (isnumeric (A)) && (issquare (A))
- is_matrix_A = true;
- if (size (A, 2) != size (b, 1))
- error ("dimension of b is not consistent with A")
- endif
- elseif (isa (A, "function_handle"))
- is_matrix_A = false;
- elseif (ischar (A))
- A = str2func (A);
- is_matrix_A = false;
- else
- error("__checkandstring__: A must be a square matrix or a function handle")
- endif
-
-## Check if there exist preconditioners
-
- is_prec = true;
- is_matrix_M1 = is_matrix_M2 = is_matrix_A;
-
-## Check M1 and sets its type
- if (isempty (M1)) # M1 empty, then no preconditioners
- is_prec = false;
- else # M1 not empty, I check its type and consistency
- if (isnumeric (M1) && issquare (M1))
- is_matrix_M1 = true;
- elseif (isa (M1, "function_handle"))
- is_matrix_M1 = false;
- elseif (ischar (M1))
- M1 = str2func (M1);
- is_matrix_M1 = false;
- else
- error("__checkandstring__: M1 must be a square matrix or a function handle")
- endif
- if (isempty (M2)) # M2 empty, then I set its as is_matrix_A
- is_matrix_M2 = is_matrix_A;
- else # M2 not empty, then I check its type
- if (isnumeric (M2) && issquare (M2))
- is_matrix_M2 = true;
- elseif (isa (M2, "function_handle"))
- is_matrix_M2 = false;
- elseif (ischar (M2))
- M2 = str2func (M2);
- is_matrix_M2 = false;
- else
- error("__checkandstring__: M2 must be a square matrix or a function handle")
- endif
- endif
- endif
-endfunction
diff --git a/scripts/sparse/private/__default__input__.m b/scripts/sparse/private/__default__input__.m
deleted file mode 100644
--- a/scripts/sparse/private/__default__input__.m
+++ /dev/null
@@ -1,49 +0,0 @@
-## Copyright (C) 2016 Cristiano Dorigo, Octave Arena
-##
-## 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 {} {@var{[varargout]} =} __default__input__ (@var{def_val}, @var{varargin})
-## Check if the arguments in input of a function are empty or missing and in such cases sets up them in default values
-##
-## The input argoments are:
-## @itemize @minus
-## @item @var{def_val} is a cell array that contains the values to use as default
-## @item @var{varargin} are the input argument
-## @end itemize
-##
-## The output argoments:
-## @itemize @minus
-## @item @var{varargout} all the input argument with filled the empty or missing paramenters
-##
-## @end itemize
-##
-## @end deftypefn
-
-
-function [varargout] = __default__input__ (def_val, varargin)
-
- m = length (def_val);
- n = length (varargin);
-
- for i = 1:m
- if (n < i || isempty (varargin {i}))
- varargout {i} = def_val {i};
- else
- varargout {i} = varargin {i};
- endif
- endfor
diff --git a/scripts/sparse/private/bicg_function.m b/scripts/sparse/private/bicg_function.m
deleted file mode 100644
--- a/scripts/sparse/private/bicg_function.m
+++ /dev/null
@@ -1,105 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- bicg_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, ...
- M2, varargin)
- ## bicg: function handle case
- x = x_min = x_pr = x0;
- c = b;
- iter = iter_min = 0;
- flag = 1; # Default flag is "maximum number of iterations reached"
- resvec = zeros (maxit + 1, 1);
- r0 = b - feval (A, x, "notransp", varargin{:}); # Residual of the sytem
- s0 = c - feval (A, x, "transp", varargin{:}); # Residual of the "dual system"
- resvec (1) = norm (r0, 2);
-
- if (!is_prec) # No preconditioners
- p = r0; # Direction of the system
- q = s0; # Direction of the "dual system"
- while (iter < maxit) && (resvec (iter + 1) >= norm_b *tol)
- v = feval (A, p, "notransp", varargin{:});
- prod_qv = q' * v;
- if (prod_qv == 0)
- flag = 4;
- break
- endif
- alpha = (s0' * r0) / prod_qv;
- x += alpha * p;
- prod_rs = (s0' * r0); # Product between r0 and s0
- r0 -= alpha * v;
- s0 -= conj (alpha) * feval (A, q, "transp", varargin{:});
- iter += 1;
- resvec (iter + 1) = norm (r0);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm (x) * eps)
- flag = 3;
- break;
- endif
- if (prod_rs == 0)
- flag = 4;
- break;
- endif
- beta = (s0' * r0) / prod_rs;
- p = r0 + beta*p;
- q = s0 + conj (beta)*q;
- endwhile
- else # There are preconditioners
- exist_M2 = (!isempty (M2));
- try
- warning("error", "Octave:singular-matrix", "local")
- prec_r0 = feval (M1, r0, "notransp", varargin{:}); # r0 preconditioned
- prec_s0 = s0;
- if (exist_M2)
- prec_r0 = feval (M2, prec_r0, "notransp", varargin{:});
- prec_s0 = feval (M2, prec_s0, "transp", varargin{:});
- endif
- prec_s0 = feval (M1, prec_s0, "transp", varargin{:}); # s0 preconditioned
- p = prec_r0; # Direction of the system
- q = prec_s0; # Direction of the "dual system"
- catch
- lasterror
- flag = 2;
- end_try_catch
-
- while ((flag != 2) && (iter < maxit) && ...
- (resvec (iter + 1) >= norm_b * tol))
- v = feval (A, p, "notransp", varargin{:});
- prod_qv = q' * v;
- if (prod_qv == 0)
- flag = 4;
- break
- endif
- alpha = (s0' * prec_r0) / prod_qv;
- x += alpha * p;
- prod_rs = (s0' * prec_r0); # Product between r0 and s0
- r0 -= alpha * v;
- s0 -= conj (alpha) * feval (A, q, "transp", varargin{:});
- prec_r0 = feval (M1, r0, "notransp", varargin{:});
- prec_s0 = s0;
- if (exist_M2)
- prec_r0 = feval (M2, prec_r0, "notransp", varargin{:});
- prec_s0 = feval (M2, prec_s0, "transp", varargin{:});
- endif
- prec_s0 = feval (M1, prec_s0, "transp", varargin{:});
- iter += 1;
- resvec (iter + 1) = norm (r0);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm (x) * eps)
- flag = 3;
- break;
- endif
- if (prod_rs == 0)
- flag = 4;
- break;
- endif
- beta = (s0' * prec_r0) / prod_rs;
- p = prec_r0 + beta*p;
- q = prec_s0 + conj (beta) * q;
- endwhile
- endif
- resvec = resvec (1:iter+1,1);
-endfunction
diff --git a/scripts/sparse/private/bicg_matrix.m b/scripts/sparse/private/bicg_matrix.m
deleted file mode 100644
--- a/scripts/sparse/private/bicg_matrix.m
+++ /dev/null
@@ -1,110 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- bicg_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2)
- ## bicg: matrix case
-
- ## x approximation of the system
- ## y approximation of the "dual system"
- ## x_pr approximation at the previous iteration
- ## x_min iterated which has the minimum residual
- x = y = x_min = x0;
- c = b; ## rhs for the "dual system"
- x = x_min = x_pr = x0;
- c = b;
- iter = iter_min = 0;
- flag = 1; # Default flag is "maximum number of iterations reached"
- resvec = zeros (maxit + 1, 1);
-
- r0 = b - A * x; # Residual of the sytem
- s0 = c - A' * x; # Residual of the "dual system"
- resvec (1) = norm (r0, 2);
-
- if (!is_prec) # No preconditioners
- p = r0; # Direction of the system
- q = s0; # Direction of the "dual system"
- while (iter < maxit) && (resvec (iter + 1) >= norm_b * tol)
- v = A * p;
- prod_qv = q' * v;
- if (prod_qv == 0)
- flag = 4;
- break
- endif
- alpha = (s0' * r0) / prod_qv;
- x += alpha * p;
- prod_rs = (s0' * r0); # Product between r0 and s0
- r0 -= alpha * v;
- s0 -= conj (alpha) * A' * q;
- iter += 1;
- resvec (iter + 1) = norm (r0);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm (x) * eps)
- flag = 3;
- break;
- endif
- if (prod_rs == 0)
- flag = 4;
- break;
- endif
- beta = (s0' * r0) / prod_rs;
- p = r0 + beta * p;
- q = s0 + conj (beta) * q;
- endwhile
- else # There are preconditioners
- exist_M2 = (!isempty (M2));
- try
- warning("error", "Octave:singular-matrix", "local")
- prec_r0 = M1 \ r0; # r0 preconditioned
- prec_s0 = s0;
- if (exist_M2)
- prec_r0 = M2 \ prec_r0;
- prec_s0 = M2' \ prec_s0;
- endif
- prec_s0 = M1' \ prec_s0; # s0 preconditioned
- p = prec_r0; # Direction of the system
- q = prec_s0; # Direction of the "dual system"
- catch
- flag = 2;
- end_try_catch
- while ((flag != 2) && (iter < maxit) && ...
- (resvec (iter + 1) >= norm_b * tol))
- v = A * p;
- prod_qv = q' * v;
- if (prod_qv == 0)
- flag = 4;
- break
- endif
- alpha = (s0' * prec_r0) / prod_qv;
- x += alpha * p;
- prod_rs = (s0' * prec_r0); # Product between r0 and s0
- r0 -= alpha * v;
- s0 -= conj (alpha) * A' * q;
- prec_r0 = M1 \ r0;
- prec_s0 = s0;
- if (exist_M2)
- prec_r0 = M2 \ prec_r0;
- prec_s0 = M2' \ prec_s0;
- endif
- prec_s0 = M1' \ prec_s0;
- iter += 1;
- resvec (iter + 1) = norm (r0);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm (x) * eps)
- flag = 3;
- break;
- endif
- if (prod_rs == 0)
- flag = 4;
- break;
- endif
- beta = (s0' * prec_r0) / prod_rs;
- p = prec_r0 + beta * p;
- q = prec_s0 + conj (beta)*q;
- endwhile
- endif
- resvec = resvec (1:iter+1,1);
-endfunction
diff --git a/scripts/sparse/private/bicgstab_function.m b/scripts/sparse/private/bicgstab_function.m
deleted file mode 100644
--- a/scripts/sparse/private/bicgstab_function.m
+++ /dev/null
@@ -1,140 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- bicgstab_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, ...
- varargin);
-
-## bicgstab: Function handle case
-
- ## Double maxit to mind also the "half iterations"
- d_maxit = 2 * maxit;
- iter = iter_min = 0;
- resvec = zeros (d_maxit,1);
- x = x_min = x_pr = x0;
- iter = iter_min = 0;
- ## deafult setting of flag is 1 (i.e. max number of iterations reached)
- flag = 1;
-
- res = b - feval (A, x, varargin{:});
- rr = p = res; ## rr is r_star
- rho_1 = rr' * res;
- resvec (1) = norm (res,2);
- real_tol = norm_b * tol;
-
- if (!is_prec) # No preconditioners
- while (iter < d_maxit) && (resvec (iter + 1) >= real_tol)
- v = feval (A, p, varargin{:});
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break;
- endif
- alpha = rho_1 / (prod_tmp);
- x += alpha * p;
- s = res - alpha * v;
- iter += 1;
- resvec (iter + 1) = norm (s, 2);
- if (resvec (iter + 1) <= real_tol) # reached the tol
- x_min = x;
- iter_min = iter;
- break
- elseif (resvec (iter + 1) <= resvec (iter_min + 1)) # Found min residual
- x_min = x;
- iter_min = iter;
- endif
- t = feval (A, s, varargin{:});
- omega = (t' * s) / (t' * t);
- if (omega == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- x += omega * s;
- res = s - omega * t;
- iter += 1;
- resvec (iter + 1) = norm (res);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm(x) * eps)
- flag = 3;
- break
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * res;
- if (rho_1 == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- beta = (rho_1 / rho_2) * (alpha / omega);
- p = res + beta * (p - omega * v);
- endwhile
- else # there are preconditioners
- exist_M2 = (!isempty (M2)); # I set a flag for the existence of M2
- ## To check if the preconditioners are singular or tehey have some NaN
- warning("error", "Octave:singular-matrix", "local");
- try
- p_hat = feval (M1, p, varargin{:});
- if (exist_M2)
- p_hat = feval (M2, p_hat, varargin{:});
- endif
- catch
- flag = 2;
- end_try_catch
- while (flag !=2) && (iter < d_maxit) && (resvec (iter + 1) >= real_tol)
- v = feval (A, p_hat, varargin{:});
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break;
- endif
- alpha = rho_1 / (prod_tmp);
- x += alpha * p_hat;
- s = res - alpha * v;
- iter += 1;
- resvec (iter+1) = norm (s,2);
- if (resvec (iter + 1) <= real_tol) # reached the tol
- x_min = x;
- iter_min = iter;
- break
- elseif (resvec (iter + 1) <= resvec (iter_min + 1))# Found min residual
- x_min = x;
- iter_min = iter;
- endif
- s_hat = feval (M1, s, varargin{:});
- if (exist_M2)
- s_hat = feval (M2, s_hat, varargin{:});
- endif
- t = feval (A, s_hat, varargin{:});
- omega = (t' * s) / (t' * t);
- if (omega == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- x += omega * s_hat;
- res = s - omega * t;
- iter += 1;
- resvec (iter + 1) = norm (res);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm (x) * eps)
- flag = 3;
- break
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * res;
- if (rho_1 == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- beta = (rho_1 / rho_2) * (alpha / omega);
- p = res + beta * (p - omega* v);
- p_hat = feval (M1, p, varargin{:});
- if (exist_M2)
- p_hat = feval (M2, p_hat, varargin{:});
- endif
- endwhile
- endif
-endfunction
diff --git a/scripts/sparse/private/bicgstab_matrix.m b/scripts/sparse/private/bicgstab_matrix.m
deleted file mode 100644
--- a/scripts/sparse/private/bicgstab_matrix.m
+++ /dev/null
@@ -1,140 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- bicgstab_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
-
- ## bicgstab: Matrix case
-
- ## Double maxit to mind also the "half iterations"
- d_maxit = 2 * maxit;
- iter = iter_min = 0;
- resvec = zeros (d_maxit,1);
- x = x_min = x_pr = x0;
- iter = iter_min = 0;
- ## deafult setting of flag is 1 (i.e. max number of iterations reached)
- flag = 1;
-
- res = b - A * x;
- rr = p = res; # rr is r_star
- rho_1 = rr' * res;
- resvec (1) = norm (res,2);
- real_tol = norm_b * tol;
-
- if (!is_prec) # No preconditioners
- while (iter < d_maxit) && (resvec (iter + 1) >= real_tol)
- v = A * p;
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break;
- endif
- alpha = rho_1 / (prod_tmp);
- x += alpha * p;
- s = res - alpha * v;
- iter += 1;
- resvec (iter + 1) = norm (s, 2);
- if (resvec (iter + 1) <= real_tol) # reached the tol
- x_min = x;
- iter_min = iter;
- break
- elseif ((resvec (iter + 1) <= resvec (iter_min + 1))) # Found min residual
- x_min = x;
- iter_min = iter;
- endif
- t = A * s;
- omega = (t' * s) / (t' * t);
- if (omega == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- x += omega * s;
- res = s - omega * t;
- iter += 1;
- resvec (iter + 1) = norm (res);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm(x) * eps)
- flag = 3;
- break
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * res;
- if (rho_1 == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- beta = (rho_1 / rho_2) * (alpha / omega);
- p = res + beta * (p - omega * v);
- endwhile
- else # there are preconditioners
- exist_M2 = !isempty (M2); # I set a flag for the existence of M2
- # To check if the preconditioners are singular or tehey have some NaN
- warning("error", "Octave:singular-matrix", "local");
- try
- p_hat = M1 \ p;
- if (exist_M2)
- p_hat = M2 \ p_hat;
- endif
- catch
- flag = 2;
- end_try_catch
- while ((flag != 2) && (iter < d_maxit) && (resvec (iter + 1) >= real_tol))
- v = A * p_hat;
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break;
- endif
- alpha = rho_1 / (prod_tmp);
- x += alpha * p_hat;
- s = res - alpha * v;
- iter += 1;
- resvec (iter+1) = norm (s,2);
- if (resvec (iter + 1) <= real_tol) # reached the tol
- x_min = x;
- iter_min = iter;
- break
- elseif (resvec (iter + 1) <= resvec (iter_min + 1)) # Found min residual
- x_min = x;
- iter_min = iter;
- endif
- s_hat = M1 \ s;
- if (exist_M2)
- s_hat = M2 \ s_hat;
- endif
- t = A * s_hat;
- omega = (t' * s) / (t' * t);
- if (omega == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- x += omega * s_hat;
- res = s - omega * t;
- iter += 1;
- resvec (iter + 1) = norm (res);
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x - x_pr) <= norm (x) * eps)
- flag = 3;
- break
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * res;
- if (rho_1 == 0) # x and residual don't change and the next it will be NaN
- flag = 4;
- break;
- endif
- beta = (rho_1 / rho_2) * (alpha / omega);
- p = res + beta * (p - omega * v);
- p_hat = M1 \ p;
- if (exist_M2)
- p_hat = M2 \ p_hat;
- endif
- endwhile
- endif
- resvec = resvec (1:iter+1,1);
-endfunction
diff --git a/scripts/sparse/private/cgs_function.m b/scripts/sparse/private/cgs_function.m
deleted file mode 100644
--- a/scripts/sparse/private/cgs_function.m
+++ /dev/null
@@ -1,105 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- cgs_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, ...
- varargin);
- ## cgs: function handle case
-
- flag = 1; # Default flag is 1, i.e. maximum number of iterations reached
- iter = iter_min = 0;
- x = x_min = x_pr = x0;
- ## x approximation at the actual iteration
- ## x_min approximation with the minimum residual
- ## x_pr approximation at the previous iteration (to check stagnation)
-
- r0 = rr = u = p = b - feval (A, x, varargin{:});
- resvec (1) = norm (r0, 2);
- rho_1 = rr' * r0;
-
- if (!is_prec) # No preconditioners
- while ((iter < maxit) && (resvec (iter + 1) >= tol * norm_b))
- v = feval (A, p, varargin{:});
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break
- endif
- alpha = rho_1 / prod_tmp;
- q = u - alpha * v;
- x += alpha * (u + q);
- r0 -= alpha * feval (A, u + q, varargin {:});
- iter += 1;
- resvec (iter + 1) = norm (r0, 2);
- if (norm (x - x_pr, 2) <= norm(x, 2) * eps) # Stagnation
- flag = 3;
- break
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1)) # Check min residual
- x_min = x;
- iter_min = iter;
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * r0;
- if (rho_1 == 0)
- ## if rho_1 = 0, at the actual iteration x doesn't change, and in the
- ## next will be a division by zero
- flag = 4;
- break
- endif
- beta = rho_1 / rho_2;
- u = r0 + beta * q;
- p = u + beta * (q + beta * p);
- endwhile
- else # There are preconditioners
- exist_M2 = !isempty (M2);
- try
- warning ("error","Octave:singular-matrix","local")
- p_hat = feval (M1, p, varargin{:});
- if (exist_M2)
- p_hat = feval (M2, p_hat, varargin {:});
- endif
- catch
- flag = 2;
- end_try_catch
- while ((flag != 2) && (iter < maxit) && (resvec (iter + 1) >= tol * norm_b))
- v = feval (A, p_hat, varargin{:});
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break;
- endif
- alpha = rho_1 / prod_tmp;
- q = u - alpha * v;
- u_hat = feval(M1, u + q, varargin{:});
- if (exist_M2)
- u_hat = feval (M2, u_hat, varargin{:});
- endif
- x += alpha*u_hat;
- r0 -= alpha* feval (A, u_hat, varargin{:});
- iter += 1;
- resvec (iter + 1) = norm (r0, 2);
- if (norm (x - x_pr, 2) <= norm(x, 2) * eps) # Stagnation
- flag = 3;
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1)) # Check min residual
- x_min = x;
- iter_min = iter;
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * r0;
- if (rho_1 == 0)
- flag = 4;
- break;
- endif
- beta = rho_1 / rho_2;
- u = r0 + beta * q;
- p = u + beta * (q + beta * p);
- p_hat = feval (M1, p, varargin {:});
- if (exist_M2)
- p_hat = feval (M2, p_hat, varargin{:});
- endif
- endwhile
- endif
- resvec = resvec (1: (iter + 1));
-endfunction
diff --git a/scripts/sparse/private/cgs_matrix.m b/scripts/sparse/private/cgs_matrix.m
deleted file mode 100644
--- a/scripts/sparse/private/cgs_matrix.m
+++ /dev/null
@@ -1,104 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- cgs_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
- ## cgs: matrix case
-
- flag = 1; # Default flag is 1, i.e. maximum number of iterations reached
- iter = iter_min = 0;
- x = x_min = x_pr = x0;
- ## x approximation at the actual iteration
- ## x_min approximation with the minimum residual
- ## x_pr approximation at the previous iteration (to check stagnation)
-
- r0 = rr = u = p = b - A * x0;
- resvec (1) = norm (r0, 2);
- rho_1 = rr' * r0;
-
- if (!is_prec) # No preconditioners
- while ((iter < maxit) && (resvec (iter + 1) >= tol * resvec (1)))
- v = A * p;
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break
- endif
- alpha = rho_1 / prod_tmp;
- q = u - alpha * v;
- x += alpha * (u + q);
- r0 -= alpha * A * (u + q);
- iter += 1;
- resvec (iter + 1) = norm (r0, 2);
- if (norm (x - x_pr, 2) <= norm(x, 2) * eps) # Stagnation
- flag = 3;
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1)) # Check min residual
- x_min = x;
- iter_min = iter;
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * r0;
- if (rho_1 == 0)
- ## if rho_1 = 0, at the actual iteration x doesn't change,
- ## and in the next will be a division by zero
- flag = 4;
- break;
- endif
- beta = rho_1 / rho_2;
- u = r0 + beta * q;
- p = u + beta * (q + beta * p);
- endwhile
- else # There are preconditioners
- exist_M2 = !isempty (M2);
- try
- warning ("error","Octave:singular-matrix","local")
- p_hat = M1 \ p;
- if (exist_M2)
- p_hat = M2 \ p_hat;
- endif
- catch
- flag = 2;
- end_try_catch
- while ((flag != 2) && (iter < maxit) && (resvec (iter + 1) >= tol * norm_b))
- v = A*p_hat;
- prod_tmp = (rr' * v);
- if (prod_tmp == 0)
- flag = 4;
- break;
- endif
- alpha = rho_1 / prod_tmp;
- q = u - alpha * v;
- u_hat = M1 \ (u + q);
- if (exist_M2)
- u_hat = M2 \ u_hat;
- endif
- x += alpha * u_hat;
- r0 -= alpha *A * u_hat;
- iter += 1;
- resvec (iter + 1) = norm (r0);
- if (norm (x - x_pr) <= norm(x) * eps) # Stagnation
- flag = 3;
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1)) # Check min residual
- x_min = x;
- iter_min = iter;
- endif
- x_pr = x;
- rho_2 = rho_1;
- rho_1 = rr' * r0;
- if (rho_1 == 0)
- flag = 4;
- break
- endif
- beta = rho_1 / rho_2;
- u = r0 + beta * q;
- p = u + beta * (q + beta * p);
- p_hat = M1 \ p;
- if (exist_M2)
- p_hat = M2 \ p_hat;
- endif
- endwhile
- endif
- resvec = resvec (1: (iter + 1));
-endfunction
diff --git a/scripts/sparse/private/gmres_function.m b/scripts/sparse/private/gmres_function.m
deleted file mode 100644
--- a/scripts/sparse/private/gmres_function.m
+++ /dev/null
@@ -1,132 +0,0 @@
-function [x_min, flag, it, iter, iter_min, resvec] = ...
- gmres_function (A, b, prec_b_norm, x0, restart, rtol, maxit, ...
- max_iter_number, is_prec, M1, M2, varargin)
-
- ## gmres: function handle case
-
- x_old = x_pr = x_min = x = x0;
-
- B = zeros (restart + 1, 1);
- V = zeros (rows (x), restart);
- H = zeros (restart + 1, restart);
-
- iter = 1; # total number of iterations
- iter_min = 0; # iteration with minimum residual
- outer_it = 1; # number of outer iterations
- restart_it = 1; # number of inner iterations
- it = zeros(1, 2);
- resvec = zeros (max_iter_number + 1, 1);
- flag = 1; # Default flag is maximum # of iterations exceeded
-
- if (!is_prec) # there are no preconditioners
- ## begin loop
- prec_res = b - feval (A, x_old, varargin{:});
- presn = norm (prec_res, 2);
- resvec (1) = presn;
- B (1) = presn;
- V (:, 1) = prec_res / presn;
- while (iter <= max_iter_number && presn > rtol * prec_b_norm)
- ## restart
- if (restart_it > restart)
- restart_it = 1;
- outer_it += 1;
- x_old = x;
- prec_res = b - feval (A, x_old, varargin{:});
- presn = norm (prec_res, 2);
- B (1) = presn;
- H (:) = 0;
- V (:, 1) = prec_res / presn;
- endif
- ## basic iteration
- tmp = feval (A, V (:, restart_it), varargin{:});
- [V(:, restart_it + 1), H(1:restart_it + 1, restart_it)] = ...
- mgorth (tmp, V(:,1:restart_it));
- Y = (H(1:restart_it + 1, 1:restart_it) \ B(1:restart_it + 1));
- little_res = B(1:restart_it + 1) - ...
- H(1:restart_it + 1, 1:restart_it) * Y(1:restart_it);
- presn = norm (little_res, 2);
- x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
- resvec (iter + 1) = presn;
- if (norm (x - x_pr) <= eps * norm (x))
- flag = 3; # Stagnation: little change between iterations
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- it = [outer_it, restart_it];
- endif
- x_pr = x;
- restart_it += 1;
- iter += 1;
- endwhile
-
- else # there are preconditioners
- ## begin loop
- exist_M2 = (!isempty (M2));
- u = feval (A, x_old, varargin{:});
- try
- warning("error", "Octave:singular-matrix", "local")
- prec_res = feval (M1, b - u, varargin{:}); # M1*(b-u)
- if (exist_M2)
- prec_res = feval (M2, prec_res, varargin{:});
- endif
- presn = norm (prec_res, 2);
- resvec(1) = presn;
- z = feval (M1, b, varargin{:});
- if (exist_M2)
- z = feval (M2, z, varargin{:});
- endif
- prec_b_norm = norm (z, 2);
- B (1) = presn;
- V(:, 1) = prec_res / presn;
- catch
- flag = 2;
- end_try_catch
-
- while (flag != 2) && (iter <= max_iter_number) && ...
- (presn > rtol * prec_b_norm)
- ## restart
- if (restart_it > restart)
- restart_it = 1;
- outer_it += 1;
- x_old = x;
- u = feval (A, x_old, varargin{:});
- prec_res = feval (M1, b - u, varargin{:});
- if (exist_M2)
- prec_res = feval (M2, prec_res, varargin{:});
- endif
- presn = norm (prec_res, 2);
- B(1) = presn;
- H(:) = 0;
- V(:, 1) = prec_res / presn;
- endif
- ## basic iteration
- u = feval (A, V(:, restart_it), varargin{:});
- tmp = feval (M1, u, varargin{:});
- if (exist_M2)
- tmp = feval (M2, tmp, varargin{:});
- endif
- [V(:,restart_it + 1), H(1:restart_it + 1, restart_it)] = ...
- mgorth (tmp, V(:,1:restart_it));
- Y = (H(1:restart_it + 1, 1:restart_it) \ B(1:restart_it + 1));
- little_res = B(1:restart_it + 1) - ...
- H(1:restart_it + 1, 1:restart_it) * Y(1:restart_it);
- presn = norm (little_res, 2);
- x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
- resvec(iter + 1) = presn;
- if (norm (x - x_pr) <= eps*norm (x))
- flag = 3; # Stagnation: little change between iterations
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- it = [outer_it, restart_it];
- endif
- x_pr = x;
- restart_it += 1;
- iter += 1;
- endwhile
- endif
-endfunction
diff --git a/scripts/sparse/private/gmres_matrix.m b/scripts/sparse/private/gmres_matrix.m
deleted file mode 100644
--- a/scripts/sparse/private/gmres_matrix.m
+++ /dev/null
@@ -1,133 +0,0 @@
-function [x_min, flag, it, iter, iter_min, resvec] = ...
- gmres_matrix (A, b, prec_b_norm, x0, restart, rtol, maxit, max_iter_number, is_prec, M1, M2)
-
-## gmres: matrix case
-
- x_old = x_pr = x_min = x = x0;
-
- B = zeros (restart + 1, 1);
- V = zeros (rows (x), restart);
- H = zeros (restart + 1, restart);
-
- iter = 1; # total number of iterations
- iter_min = 0; # iteration with minimum residual
- outer_it = 1; # number of outer iterations
- restart_it = 1; # number of inner iterations
- it = zeros(1, 2);
- resvec = zeros (max_iter_number + 1, 1);
- flag = 1; # Default flag is maximum # of iterations exceeded
-
- if (!is_prec)
- ## begin loop
- prec_res = b - A * x_old;
- presn = norm (prec_res, 2);
- resvec(1) = presn;
- prec_b_norm = norm (b, 2);
- B (1) = presn;
- V (:, 1) = prec_res / presn;
-
- while (iter <= max_iter_number && presn > rtol * prec_b_norm)
- ## restart
- if (restart_it > restart)
- restart_it = 1;
- outer_it += 1;
- x_old = x;
- prec_res = b - A * x_old;
- presn = norm (prec_res, 2);
- B(1) = presn;
- H(:) = 0;
- V(:, 1) = prec_res / presn;
- endif
-
- ## basic iteration
- tmp = A* V(:, restart_it);
- [V(:,restart_it + 1), H(1:restart_it + 1, restart_it)] = ...
- mgorth (tmp, V(:,1:restart_it));
- Y = (H(1:restart_it + 1, 1:restart_it) \ B(1:restart_it + 1));
- little_res = B(1:restart_it + 1) - ...
- H(1:restart_it + 1, 1:restart_it) * Y(1:restart_it);
- presn = norm (little_res, 2);
- x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
- resvec(iter + 1) = presn;
-
- if (norm (x - x_pr) <= eps * norm (x))
- flag = 3; # Stagnation: little change between iterations
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- it = [outer_it, restart_it];
- endif
-
- x_pr = x;
- restart_it += 1;
- iter += 1;
- endwhile
-
- else # there are preconditioners
- ## begin loop
- exist_M2 = (!isempty (M2));
- try
- warning("error", "Octave:singular-matrix", "local");
- prec_res = M1 \ (b - A * x_old);
- if (exist_M2)
- prec_res = M2 \ prec_res;
- endif
- presn = norm (prec_res, 2);
- resvec (1) = presn;
- z = M1 \ b;
- if (exist_M2)
- z = M2 \ z;
- endif
- prec_b_norm = norm (z, 2);
- B(1) = presn;
- V(:, 1) = prec_res / presn;
- catch
- flag = 2;
- end_try_catch
- while (flag != 2) && (iter <= max_iter_number && ...
- presn > rtol * prec_b_norm)
- ## restart
- if (restart_it > restart)
- restart_it = 1;
- outer_it += 1;
- x_old = x;
- prec_res = M1 \ (b - A * x_old);
- if (exist_M2)
- prec_res = M2 \ prec_res;
- endif
- presn = norm (prec_res, 2);
- B (1) = presn;
- H (:) = 0;
- V (:, 1) = prec_res / presn;
- endif
- ## basic iteration
- tmp = M1 \ (A * V(:, restart_it));
- if (exist_M2)
- tmp = M2 \ tmp;
- endif
- [V(:,restart_it + 1), H(1:restart_it + 1, restart_it)] = ...
- mgorth (tmp, V(:,1:restart_it));
- Y = (H(1:restart_it + 1, 1:restart_it) \ B(1:restart_it + 1));
- little_res = B(1:restart_it + 1) - ...
- H(1:restart_it + 1, 1:restart_it) * Y(1:restart_it);
- presn = norm (little_res, 2);
- x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
- resvec(iter + 1) = presn;
- if (norm (x - x_pr) <= eps * norm (x))
- flag = 3; # Stagnation: little change between iterations
- break;
- endif
- if (resvec (iter + 1) <= resvec (iter_min + 1))
- x_min = x;
- iter_min = iter;
- it = [outer_it, restart_it];
- endif
- x_pr = x;
- restart_it += 1;
- iter += 1;
- endwhile
- endif
-
-endfunction
diff --git a/scripts/sparse/private/pcg_function.m b/scripts/sparse/private/pcg_function.m
deleted file mode 100644
--- a/scripts/sparse/private/pcg_function.m
+++ /dev/null
@@ -1,166 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec, eigest] = ...
- pcg_function (A, b, b_norm, x0, tol, maxit, is_prec, m1, m2,...
- n_arg_out, varargin)
-
-## pcg: function handle case
-
- x = x_pr = x_min = x0;
- ## x_pr (x previous) needs to check the stagnation
- ## x_min needs to save the iterated with minimum residual
-
- r = b - feval (A, x, varargin{:});
- iter = 2;
- iter_min = 0;
- flag = 1;
- resvec = zeros (maxit + 1, 2);
- resvec(1, 1) = norm (r);
- p = zeros (size (b));
- alpha = old_tau = 1;
-
- if (n_arg_out > 5)
- T = zeros (maxit, maxit);
- else
- T = [];
- endif
-
- if (!is_prec) # There are no preconditioners
- while (resvec(iter - 1, 1) > tol * b_norm && iter < maxit)
- tau = r' * r;
- resvec(iter - 1, 2) = sqrt (tau);
- beta = tau / old_tau;
- old_tau = tau;
- p = r + beta * p;
- w = feval (A, p, varargin{:});
- ## Needed only for eigest
- old_alpha = alpha;
- den = p' * w; # denominator of alpha
- alpha = tau / den;
- ## Check if alpha is negative and/or if it has a consistent
- ## imaginary part: if yes then A probably is not positive definite
- if ((tau <= 0) || (abs (imag (den)) >= eps (real (den)) * tol / eps) || ...
- (real (den) <= 0))
- flag = 4;
- break;
- endif
-
- x += alpha * p;
- r -= alpha * w;
- resvec(iter, 1) = norm (r);
- ## Chek if the iterated has minimum residual
- if (resvec (iter,1) <= resvec (iter_min + 1,1))
- x_min = x;
- iter_min = iter - 1;
- endif
- if (n_arg_out > 5 && iter > 2)
- T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
- [1, sqrt(beta); sqrt(beta), beta] ./ ...
- old_alpha;
- endif
- iter += 1;
- if (norm (x - x_pr) <= eps * norm (x)) # Check the stagnation
- flag = 3;
- break;
- endif
- x_pr = x;
- endwhile
-
- if (n_arg_out > 5)
- z = r; # Needs only for eig
- endif
-
- else # There are the preconditioners
- exist_m2 = (!isempty (m2));
- while (resvec(iter-1,1) > tol * b_norm && iter < maxit)
- if (iter == 2) # Check whether M1 or M2 are singular
- try
- warning("error","Octave:singular-matrix","local")
- z = feval (m1, r, varargin{:});
- if (exist_m2)
- z = feval (m2, z, varargin{:});
- endif
- catch
- flag = 2;
- break;
- end_try_catch
- else
- z = feval (m1, r, varargin{:});
- if (exist_m2)
- z = feval (m2, z, varargin{:});
- endif
- endif
-
- tau = z' * r;
- resvec(iter - 1, 2) = sqrt (tau);
- beta = tau / old_tau;
- old_tau = tau;
- p = z + beta * p;
- w = feval (A, p, varargin{:});
- ## Needed only for eigest.
- old_alpha = alpha;
- den = p' * w;
- alpha = tau / den;
-
- ## Check if alpha is negative and/or if it has a consistent
- ## imaginary part: if yes then A probably is not positive definite
- if ((abs (imag (tau)) >= eps (real (tau)) * tol / eps) || ...
- real (tau) <= 0 || ...
- (abs (imag (den)) >= eps (real (den)) * tol / eps) || ...
- (real (den) <= 0))
- flag = 4;
- break;
- endif
-
- x += alpha * p;
- r -= alpha * w;
- resvec(iter, 1) = norm (r);
- ## Chek if the iterated has minimum residual
- if (resvec (iter,1) <= resvec (iter_min + 1,1))
- x_min = x;
- iter_min = iter - 1;
- endif
- if (n_arg_out > 5 && iter > 2)
- T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
- [1, sqrt(beta); sqrt(beta), beta] ./ ...
- old_alpha;
- endif
- iter += 1;
- if (norm (x - x_pr) <= eps * norm (x)) # Check the stagnation
- flag = 3;
- break;
- endif
- x_pr = x;
- endwhile
-
- if (n_arg_out > 5)
- ## Apply the preconditioner once more and finish with the precond
- ## residual.
- z = feval (m1, r, varargin{:});
- if (exist_m2)
- z = feval (m2, z, varargin{:});
- endif
- endif
- endif
-
- ## (Eventually) computes the eigenvalue of inv(m2)*inv(m1)*A
- if (n_arg_out > 5)
- if (flag != 4)
- if (iter > 3)
- T = T(2:iter-2,2:iter-2);
- l = eig (T);
- eigest = [min(l), max(l)];
- else
- eigest = [NaN, NaN];
- warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
- endif
- else
- eigest = [NaN, NaN];
- warning ('pcg: eigenvalue estimate failed: matrix not positive definite?')
- endif
- resvec(iter - 1, 2) = sqrt (r' * z);
- resvec = resvec (1:(iter-1), :);
- else
- eigest = [NaN, NaN];
- resvec = resvec(1:(iter-1),1);
- endif
-
-endfunction
diff --git a/scripts/sparse/private/pcg_matrix.m b/scripts/sparse/private/pcg_matrix.m
deleted file mode 100644
--- a/scripts/sparse/private/pcg_matrix.m
+++ /dev/null
@@ -1,162 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec, eigest] = ...
- pcg_matrix (A, b, b_norm, x0, tol, maxit, is_prec, m1, m2, n_arg_out)
-
-## pcg: matrix case
-
- x = x_pr = x_min = x0;
- ## x_pr (x previous) needs to check the stagnation
- ## x_min needs to save the iterated with minimum residual
-
- r = b - A * x;
- iter = 2;
- iter_min = 0;
- flag = 1;
- resvec = zeros (maxit + 1, 2);
- resvec(1, 1) = norm (r);
- p = zeros (size (b));
- alpha = old_tau = 1;
-
- if (n_arg_out > 5)
- T = zeros (maxit, maxit);
- else
- T = [];
- endif
-
- if (!is_prec) # There are no preconditioners
- while (resvec(iter - 1, 1) > tol * b_norm && iter < maxit)
- tau = r' * r;
- resvec(iter - 1, 2) = sqrt (tau);
- beta = tau / old_tau;
- old_tau = tau;
- p = r + beta * p;
- w = A * p;
- ## Needed only for eigest.
- old_alpha = alpha;
- den = p' * w; # denominator of alpha
- alpha = tau / den;
- ## Check if alpha is negative and/or if it has a not small
- ## imaginary part: if yes then A probably is not positive definite
- if ((tau <= 0) || (abs (imag (den)) >= eps (real (den)) * tol / eps) ...
- || (real (den) <= 0))
- flag = 4;
- break;
- endif
-
- x += alpha * p;
- r -= alpha * w;
- resvec(iter, 1) = norm (r);
- ## Check if the iterated has minimum residual
- if (resvec (iter, 1) <= resvec (iter_min + 1, 1))
- x_min = x;
- iter_min = iter - 1;
- endif
- if (n_arg_out > 5 && iter > 2)
- T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
- [1, sqrt(beta); sqrt(beta), beta] ./ ...
- old_alpha;
- endif
- iter += 1;
- if (norm (x - x_pr) <= eps * norm (x)) # Check for stagnation
- flag = 3;
- break;
- endif
- x_pr = x;
- endwhile
-
- if (n_arg_out > 5)
- z = r; # Needs only for eig
- endif
-
- ## Else there are preconditioners
- else
- exist_m2 = (!isempty (m2));
- while (resvec(iter - 1, 1) > tol * b_norm && iter < maxit)
- if (iter == 2) # check whether M1 or M2 are singular
- try
- warning("error","Octave:singular-matrix","local");
- z = m1 \ r;
- if (exist_m2)
- z = m2 \ z;
- endif
- catch
- flag = 2;
- break;
- end_try_catch
- else
- z = m1 \ r;
- if (exist_m2)
- z = m2 \ z;
- endif
- endif
- tau = z' * r;
- resvec(iter - 1, 2) = sqrt (tau);
- beta = tau / old_tau;
- old_tau = tau;
- p = z + beta * p;
- w = A * p;
- ## Needed only for eigest.
- old_alpha = alpha;
- den = p'*w;
- alpha = tau / den;
- ## Check if alpha is negative and/or if it has a consistent
- ## imaginary part: if yes then A probably is not positive definite
- if ((abs (imag (tau)) >= eps (real (tau)) * tol / eps) || ...
- real (tau) <= 0 || ...
- (abs (imag (den)) >= eps (real (den)) * tol / eps) || ...
- (real (den) <= 0))
- flag = 4;
- break;
- endif
-
- x += alpha * p;
- r -= alpha * w;
- resvec(iter, 1) = norm (r);
- ## Check if the iterated has minimum residual
- if (resvec (iter, 1) <= resvec (iter_min + 1,1))
- x_min = x;
- iter_min = iter - 1;
- endif
- if (n_arg_out > 5 && iter > 2)
- T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
- [1, sqrt(beta); sqrt(beta), beta] ./ ...
- old_alpha;
- endif
- iter += 1;
- if (norm (x - x_pr) <= eps * norm (x)) # Check the stagnation
- flag = 3;
- break;
- endif
- x_pr = x;
- endwhile
-
- if (n_arg_out > 5)
- z = m1 \ r;
- if (exist_m2)
- z = m2 \ z;
- endif
- endif
- endif
-
- ## (Eventually) computes the eigenvalue of inv(m2)*inv(m1)*A
- if (n_arg_out > 5)
- if (flag != 4)
- if (iter > 3)
- T = T(2:iter-2,2:iter-2);
- l = eig (T);
- eigest = [min(l), max(l)];
- else
- eigest = [NaN, NaN];
- warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
- endif
- else
- eigest = [NaN, NaN];
- warning ('pcg: eigenvalue estimate failed: matrix not positive definite?')
- endif
- resvec(iter - 1, 2) = sqrt (r' * z);
- resvec = resvec (1:(iter-1), :);
- else
- eigest = [NaN, NaN];
- resvec = resvec(1:(iter-1),1);
- endif
-
-endfunction
diff --git a/scripts/sparse/private/tfqmr_function.m b/scripts/sparse/private/tfqmr_function.m
deleted file mode 100644
--- a/scripts/sparse/private/tfqmr_function.m
+++ /dev/null
@@ -1,142 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- tfqmr_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, ...
- varargin);
-## tfqmr: function handle case
-
- x = x_pr = x_min = x0;
- iter = iter_min = m = 0;
- resvec = zeros (maxit, 1);
- flag = 1;
-
- w = u = r = r_star = b - feval (A, x0, varargin{:});
- rho_1 = (r_star' * r);
- d = 0;
- tau = norm (r, 2);
- theta = eta = 0;
- resvec (1, 1) = norm (r, 2);
- it = 1;
-
- if (!is_prec) # No preconditioners
- v = feval (A, u, varargin{:});
- while ((iter < maxit) && (resvec (iter + 1, 1) >= resvec (1, 1)*tol))
- if (it > 0) # iter is even
- v_r = r_star' * v; # inner prod between r_star and v
- if (v_r == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break
- endif
- alpha = rho_1 / v_r;
- u_1 = u - alpha * v; # u at the after iteration
- endif
- w -= alpha * feval (A, u, varargin{:});
- d = u + ((theta * theta) / alpha) * eta * d;
- theta = norm (w, 2) / tau;
- c = 1 / sqrt (1 + theta * theta);
- tau *= theta * c;
- eta = (c * c) * alpha;
- x += eta * d;
- r -= eta * feval (A, d, varargin{:});
- if (it < 0) # iter is odd
- rho_2 = rho_1;
- rho_1 = (r_star' * w);
- if (rho_1 == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break;
- endif
- beta = rho_1 / rho_2;
- u_1 = w + beta * u; # u at the after iteration
- v = feval (A, u_1, varargin{:}) + ...
- beta * (feval (A, u, varargin{:}) + beta * v);
- endif
- u = u_1;
- iter += 1;
- resvec(iter + 1, 1) = norm (r, 2);
- if (resvec (iter + 1, 1) <= resvec (iter_min + 1, 1))
- ## iter with min residual
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x_pr - x) <= norm (x) * eps)
- flag = 3; # Stagnation
- break
- endif
- x_pr = x;
- it = -it;
- endwhile
- else # There are Preconditioners
- exist_M2 = !isempty (M2);
- try
- warning("error", "Octave:singular-matrix", "local");
- u_hat = feval (M1, u, varargin{:});
- if (exist_M2)
- u_hat = feval (M2, u_hat, varargin{:});
- endif
- v = feval (A, u_hat, varargin{:});
- catch
- flag = 2;
- end_try_catch
- while ((flag != 2) && (iter < maxit) && ...
- (resvec (iter + 1, 1) >= norm_b * tol))
- if (it > 0) # iter is even
- v_r = r_star' * v; # inner prod between r_star and v
- if (v_r == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break
- endif
- alpha = rho_1 / v_r;
- u_1 = u - alpha * v; # u at the after iteration
- endif
- u_hat = feval (M1, u, varargin{:});
- if (exist_M2)
- u_hat = feval (M2, u_hat, varargin{:});
- endif
- w -= alpha * feval (A, u_hat, varargin{:});
- d = u_hat + ((theta * theta) / alpha) * eta * d;
- theta = norm (w, 2) / tau;
- c = 1 / sqrt (1 + theta * theta);
- tau *= theta * c;
- eta = (c * c) * alpha;
- x += eta * d;
- r -= eta * feval (A, d, varargin{:});
- if (it < 0) # iter is odd
- rho_2 = rho_1;
- rho_1 = (r_star' * w);
- if (rho_1 == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break
- endif
- beta = rho_1 / rho_2;
- u_1 = w + beta * u; # u at the after iteration
- u1_hat = feval (M1, u_1, varargin{:});
- if (exist_M2)
- u1_hat = feval (M2, u1_hat, varargin{:});
- endif
- v = feval (A, u1_hat, varargin{:}) + ...
- beta * (feval (A, u_hat, varargin{:}) + beta * v);
- endif
- u = u_1;
- iter += 1;
- resvec (iter + 1, 1) = norm (r, 2);
- if (resvec (iter + 1, 1) <= resvec (iter_min + 1, 1))
- ## iter with min residual
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x_pr - x) <= norm (x) * eps)
- flag = 3; # Stagnation
- break
- endif
- x_pr = x;
- it = -it;
- endwhile
- endif
- resvec = resvec (1: (iter + 1));
-endfunction
diff --git a/scripts/sparse/private/tfqmr_matrix.m b/scripts/sparse/private/tfqmr_matrix.m
deleted file mode 100644
--- a/scripts/sparse/private/tfqmr_matrix.m
+++ /dev/null
@@ -1,139 +0,0 @@
-function [x_min, flag, iter, iter_min, resvec] = ...
- tfqmr_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
- ## tfqmr: matrix case
-
- x = x_pr = x_min = x0;
- iter = iter_min = m = 0;
- resvec = zeros (maxit, 1);
- flag = 1;
-
- w = u = r = r_star = b - A * x0;
- rho_1 = (r_star' * r);
- d = 0;
- tau = norm (r, 2);
- theta = eta = 0;
- resvec (1, 1) = norm (r, 2);
-
- it = 1; # it will detect even or odd iterations
- if (!is_prec) # No preconditioners
- v = A * u;
- while ((iter < maxit) && (resvec (iter + 1, 1) >= resvec (1, 1) * tol))
- if (it > 0) # iter is even
- v_r = r_star' * v; # inner prod between r_star and v
- if (v_r == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break
- endif
- alpha = rho_1 / v_r;
- u_1 = u - alpha * v; # u at the after iteration
- endif
- w -= alpha * A * u;
- d = u + ((theta * theta) / alpha) * eta * d;
- theta = norm (w, 2) / tau;
- c = 1 / sqrt (1 + theta * theta);
- tau *= theta * c;
- eta = (c * c) * alpha;
- x += eta * d;
- r -= eta * A * d;
- if (it < 0) # iter is odd
- rho_2 = rho_1;
- rho_1 = (r_star' * w);
- if (rho_1 == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break;
- endif
- beta = rho_1 / rho_2;
- u_1 = w + beta * u; # u at the after iteration
- v = A * u_1 + beta * (A * u + beta * v);
- endif
- u = u_1;
- iter += 1;
- resvec(iter + 1, 1) = norm (r, 2);
- if (resvec (iter + 1, 1) <= resvec (iter_min + 1, 1))
- ## iter with min residual
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x_pr - x) <= norm (x) * eps)
- flag = 3; # Stagnation
- break
- endif
- x_pr = x;
- it = -it;
- endwhile
- else # There are Preconditioners
- exist_M2 = !isempty (M2);
- try
- warning ("error", "Octave:singular-matrix", "local");
- u_hat = M1 \ u;
- if (exist_M2)
- u_hat = M2 \ u_hat;
- endif
- v = A * u_hat;
- catch
- flag = 2;
- end_try_catch
- while ((flag != 2) && (iter < maxit) && ...
- (resvec (iter + 1, 1) >= norm_b*tol))
- if (it > 0) # iter is even
- v_r = r_star' * v; # inner prod between r_star and v
- if (v_r == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break
- endif
- alpha = rho_1 / v_r;
- u_1 = u - alpha * v; # u at the after iteration
- endif
- u_hat = M1 \ u;
- if (exist_M2)
- u_hat = M2 \ u_hat;
- endif
- w -= alpha * A * u_hat;
- d = u_hat + ((theta * theta) / alpha) * eta * d;
- theta = norm (w, 2) / tau;
- c = 1 / sqrt (1 + theta * theta);
- tau *= theta * c;
- eta = (c * c) * alpha;
- x += eta * d;
- r -= eta * A * d;
- if (it < 0) # iter is odd
- rho_2 = rho_1;
- rho_1 = (r_star' * w);
- if (rho_1 == 0)
- ## Essentially the next iteration doesn't change x,
- ## and the iter after this will have a division by zero
- flag = 4;
- break;
- endif
- beta = rho_1 / rho_2;
- u_1 = w + beta * u; # u at the after iteration
- u1_hat = M1 \ u_1;
- if (exist_M2)
- u1_hat = M2 \ u1_hat;
- endif
- v = A * u1_hat + beta * (A * u_hat + beta * v);
- endif
- u = u_1;
- iter += 1;
- resvec(iter + 1, 1) = norm(r, 2);
- if (resvec (iter + 1, 1) <= resvec (iter_min + 1, 1))
- ## iter with min residual
- x_min = x;
- iter_min = iter;
- endif
- if (norm (x_pr - x) <= norm (x) * eps)
- flag = 3; # Stagnation
- break
- endif
- x_pr = x;
- it = -it;
- endwhile
- endif
- resvec = resvec (1: (iter + 1));
-endfunction
diff --git a/scripts/sparse/tfqmr.m b/scripts/sparse/tfqmr.m
deleted file mode 100644
--- a/scripts/sparse/tfqmr.m
+++ /dev/null
@@ -1,416 +0,0 @@
-## Copyright (C) 2016 Cristiano Dorigo, Octave Arena
-##
-## 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} {@var{x} =} tfqmr (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
-## @deftypefnx {Function File} {@var{x} =} tfqmr (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
-## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} tfqmr (@var{A}, @var{b}, @dots{})
-## Solve @code{A x = b} using the Transpose-Tree qmr method, based on the cgs
-##
-## The input parameters are:
-##
-## @itemize @minus
-##
-## @item @var{A} is the square matrix of the linear system. @var{A}
-## can be passed
-## as a matrix or as a function handle or inline function @code{f} such that
-## @code{f(x) = A*x}.
-##
-## @item @var{b} is the right hand side vector. It must be a column vector
-## with the same number of rows as @var{A}.
-##
-## @item @var{rtol} is the relative tolerance, if not given or set to [] the
-## default value 1e-6 is used.
-##
-## @item @var{maxit} the maximum number of outer iterations, if not given or
-## set to [] the default value @code{min (20, numel (b))} is used. To be
-## compatible, since the method as different behaviours in the iteration
-## number is odd or even, is considered as iteration in @code{tfqmr} the
-## entire odd-even cycle. That is, to make an entire iteration, the algorithm
-## performs two sub-iterations: the odd one and the even one.
-##
-## @item @var{M1}, @var{M2} are the preconditioners. The preconditioner
-## @var{P} is given as @code{P = M1 * M2}.
-## Both @var{M1} and @var{M2} can be passed as a matrix or as a function
-## handle or inline function @code{g} such that @code{g(x) = M1 \ x} or
-## @code {g(x) = M2 \ x}.
-## The technique used is the rigth-preconditioning, i.e. it is solved
-## @code{A*inv(P)*y = b} and then @code{x = inv(P)*y}, instead of
-## @code{A x = b}.
-##
-## @item @var{x0} the initial guess, if not given or set to [] the default
-## value @code{zeros (size (b))} is used.
-##
-## @end itemize
-##
-## The arguments which follow @var{x0} are treated as parameters, and passed in
-## a proper way to any of the functions (@var{A} or @var{P}) which are passed
-## to @code{tfqmr}.
-##
-## The output parameters are:
-##
-## @itemize @minus
-##
-## @item @var{x} is the approximation computed. If the method doesn't
-## converge then it is the iterated with the minimum residual.
-##
-## @item @var{flag} indicates the exit status:
-##
-## @itemize @minus
-## @item 0: iteration converged to the within the chosen tolerance
-##
-## @item 1: the maximum number of iterations was reached before convergence
-##
-## @item 2: the preconditioner matrix is singular
-##
-## @item 3: the algorithm reached stagnation
-##
-## @item 4: the algorithm can't continue due to a division by zero
-## @end itemize
-##
-## @item @var{relres} is the relative residual obtained as
-## @code{(@var{A}*@var{x}-@var{b}) / @code{norm(@var{b})}}.
-##
-## @item @var{iter} is the iteration which @var{x} is
-## computed.
-##
-## @item @var{resvec} is a vector containing the residual at each iteration
-## (including @code{norm(b - A x0)}).
-## Doing @code{length(@var{resvec}) - 1} is possible to see the
-## total number of iterations performed.
-##
-## @end itemize
-##
-## Let us consider a trivial problem with a tridiagonal matrix
-##
-## @example
-## @group
-## n = 20;
-## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-## toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-## sparse (1, 2, 1, 1, n) * n / 2);
-## b = A * ones (n, 1);
-## restart = 5;
-## [M1, M2] = ilu (A); # in this tridiagonal case it corresponds to chol (A)'
-## M = M1 * M2;
-## Afun = @@(x) A * x;
-## Mfun = @@(x) M \ x;
-## M1fun = @@(x) M1 \ x;
-## M2fun = @@(x) M2 \ x;
-## @end group
-## @end example
-##
-## @sc{Example 1:} simplest usage of @code{tfqmr}
-##
-## @example
-## x = tfqmr (A, b, [], n)
-## @end example
-##
-## @sc{Example 2:} @code{tfqmr} with a function which computes
-## @code{@var{A} * @var{x}}
-##
-## @example
-## x = tfqmr (Afun, b, [], n)
-## @end example
-##
-## @sc{Example 3:} @code{tfqmr} with a preconditioner matrix @var{M}
-##
-## @example
-## x = tfqmr (A, b, [], 1e-06, n, M)
-## @end example
-##
-## @sc{Example 4:} @code{tfqmr} with a function as preconditioner
-##
-## @example
-## x = tfqmr (Afun, b, 1e-6, n, Mfun)
-## @end example
-##
-## @sc{Example 5:} @code{tfqmr} with preconditioner matrices @var{M1}
-## and @var{M2}
-##
-## @example
-## x = tfqmr (A, b, [], 1e-6, n, M1, M2)
-## @end example
-##
-## @sc{Example 6:} @code{tfmqr} with functions as preconditioners
-##
-## @example
-## x = tfqmr (Afun, b, 1e-6, n, M1fun, M2fun)
-## @end example
-##
-## @sc{Example 7:} @code {tfqmr} with as input a function requiring an argument
-##
-## @example
-## @group
-## function y = Ap (A, x, z) # compute A^z * x
-## y = x;
-## for i = 1:z
-## y = A * y;
-## endfor
-## endfunction
-## Apfun = @(x, string, p) Ap (A, x, string, p);
-## x = tfqmr (Apfun, b, [], [], [], [], [], 2);
-## @end group
-## @end example
-##
-## @sc{Example 8:} explicit example to show that @code{tfqmr} uses a
-## right preconditioner
-##
-## @example
-## @group
-## [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
-## M = M1 * M2;
-##
-## ## reference solution computed by tfqmr after one iteration
-## [x_ref, fl] = tfqmr (A, b, [], 1, M)
-##
-## ## rigth preconditioning
-## [y, fl] = tfqmr (A / M, b, [], 1)
-## x = M \ y # compare x and x_ref
-##
-## @end group
-## @end example
-##
-## References:
-##
-## @enumerate
-##
-## @item @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear Systems},
-## Second edition, 2003, SIAM
-##
-## @end enumerate
-##
-## @seealso{bicg, bicgstab, cgs, gmres, pcg, qmr, pcr}
-##
-## @end deftypefn
-
-function [x_min, flag, relres, iter_min, resvec] = ...
- tfqmr (A, b, tol = [], maxit = [], M1 = [], M2 = [], ...
- x0 = [], varargin)
-
- [A, M1, M2, is_prec, is_matrix_A, is_matrix_M1, is_matrix_M2] = ...
- __checkandstring__ (A, b, M1, M2);
-
- [tol, maxit, x0] = __default__input__ ({1e-06, 2 * min(20, rows (b)), ...
- zeros(rows (b), 1)}, tol, ...
- maxit, x0);
-
- maxit = 2 * maxit; # To be compatible, since iteration = odd+even ones
-
- norm_b = norm (b, 2);
- if (norm_b == 0)
- if (nargout < 2)
- printf("The right hand side vector is all zero so tfqmr \n")
- printf ("returned an all zero solution without iterating.\n")
- endif
- x_min = zeros (numel (b), 1);
- iter_min = 0;
- flag = 0;
- resvec = 0;
- relres = 0;
- return
- endif
-
- if (is_matrix_A && is_matrix_M1 && is_matrix_M2)
- [x_min, flag, iter, iter_min, resvec] = ...
- tfqmr_matrix (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2);
- else
- if (is_matrix_A)
- A = @(x) A * x;
- endif
- if (is_prec)
- if (is_matrix_M1)
- M1 = @(x) M1 \ x;
- endif
- if (!isempty (M2) && is_matrix_M2)
- M2 = @(x) M2 \ x;
- endif
- endif
- [x_min, flag, iter, iter_min, resvec] = ...
- tfqmr_function (A, b, norm_b, x0, tol, maxit, is_prec, M1, M2, ...
- varargin{:});
- endif
- relres = resvec (iter_min + 1) / norm (b);
-
- iter_min = floor(iter_min / 2); # compatibility, since it
- # makes two times the effective iterations
-
- if (relres <= tol)
- flag = 0;
- endif
-
- if (nargout < 2) # Output strings
- switch (flag)
- case {0}
- printf ("tfqmr converged at iteration %i ", iter_min);
- printf ("to a solution with relative residual %e\n", relres);
- case {1}
- printf ("tfqmr stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the maximum number of iterations was reached.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {2}
- printf ("tfqmr stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the preconditioner matrix is singular.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {3}
- printf ("tfqmr stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method stagnated.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- case {4}
- printf ("tfqmr stopped at iteration %i ", iter);
- printf ("without converging to the desired tolerance %e\n", tol);
- printf ("because the method can't continue.\n");
- printf ("The iterate returned (number %i) ", iter_min);
- printf ("has relative residual %e\n", relres);
- endswitch
- endif
-endfunction
-%!test
-%! ## Check that all the subscripts works
-%! A = magic (5);
-%! b = rand (5, 1);
-%! M1 = diag (sqrt (diag(A)));
-%! M2 = M1;
-%! maxit = 10;
-%! Afun = @(z) A*z;
-%! M1_fun = @(z) M1 \ z;
-%! M2_fun = @(z) M2 \ z;
-%! [x, flag] = tfqmr (A,b);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (A, b, [], maxit, M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (A, b, [], maxit, M1_fun, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (A, b, [], maxit, M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (A, b, [], maxit, M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (Afun, b);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (Afun, b, [], maxit, M1, M2);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (Afun, b, [], maxit, M1_fun, M2);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (Afun, b, [], maxit, M1, M2_fun);
-%! assert(flag, 0);
-%! [x, flag] = tfqmr (Afun, b, [], maxit, M1_fun, M2_fun);
-%! assert(flag, 0);
-
-%!shared A, b, n, M1, M2
-%!
-%!test
-%! n = 100;
-%! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
-%! b = sum (A, 2);
-%! tol = 1e-8;
-%! maxit = 15;
-%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
-%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
-%! [x, flag, relres, iter, resvec] = tfqmr (A, b, tol, maxit, M1, M2);
-%! assert (x, ones (size (b)), 1e-7);
-%!
-%!test
-%!function y = afun (x, a)
-%! y = a * x;
-%!endfunction
-%!
-%! tol = 1e-8;
-%! maxit = 15;
-%!
-%! [x, flag, relres, iter, resvec] = tfqmr (@(x) afun (x, A), b,
-%! tol, maxit, M1, M2);
-%! assert (x, ones (size (b)), 1e-7);
-
-%!test
-%! n = 10;
-%! tol = 1e-8;
-%! a = (2 * sprand (n, n, .1) - 1) + 1i * (2 * sprand (n, n, .1) - 1);
-%! A = a + 2 * eye (n);
-%! b = sum (A, 2);
-%! [x, flag, relres, iter, resvec] = tfqmr (A, b, tol, [], diag (diag (A)));
-%! assert (x, ones (size (b)), 1e-7);
-
-%!test
-%! ## Solve complex linear system
-%! A = [1 + 1i, 1 + 1i; 2 - 1i, 2 + 1i];
-%! b = A * [1; 1];
-%! [x, flag, relres, iter, resvec] = tfqmr (A, b, [], 3);
-%! assert (x, [1; 1], 1e-6);
-
-%!test
-%! A = diag(1:50);
-%! A (1,50) = 10000;
-%! b = ones (50,1);
-%! [x, flag, relres, iter, resvec] = tfqmr (A, b, [], 100);
-%! assert (flag, 0)
-%! assert (x, A\b, 1e-05)
-%! ## Detects a singular preconditioner
-%! M = ones(50);
-%! M(1,1) = 0;
-%! [x, flag] = tfqmr (A, b, [], 100, M);
-%! assert(flag, 2)
-
-%!demo # simplest use
-%! n = 20;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.1 * eye (n));
-%! M = M1 * M2;
-%! x = tfqmr (A, b, [], n);
-%! Afun = @(x) A * x;
-%! x = tfqmr (Afun, b, [], n);
-%! x = tfqmr (A, b, 1e-6, n, M);
-%! x = tfqmr (A, b, 1e-6, n, M1, M2);
-%! Mfun = @(z) M \ z;
-%! x = tfqmr (Afun, b, 1e-6, n, Mfun);
-%! M1fun = @(z) M1 \ z;
-%! M2fun = @(z) M2 \ z;
-%! x = tfqmr (Afun, b, 1e-6, n, M1fun, M2fun);
-%! function y = Ap (A, x, z) # compute A^z * x or (A^z)' * x
-%! y = x;
-%! for i = 1:z
-%! y = A * y;
-%! endfor
-%! endfunction
-%! Afun = @(x, p) Ap (A, x, p);
-%! x = tfqmr (Afun, b, [], 2*n, [], [], [], 2); # solution of A^2 * x = b
-
-%!demo
-%! n = 10;
-%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
-%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
-%! sparse (1, 2, 1, 1, n) * n / 2);
-%! b = A * ones (n, 1);
-%! [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
-%! M = M1 * M2;
-%!
-%! ## reference solution computed by tfqmr after one iteration
-%! [x_ref, fl] = tfqmr (A, b, [], 1, M);
-%! x_ref
-%!
-%! ## right preconditioning
-%! [y, fl] = tfqmr (A / M, b, [], 1);
-%! x = M \ y # compare x and x_ref