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