Mercurial > octave
view scripts/sparse/gmres.m @ 29358:0a5b15007766 stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 10 Feb 2021 09:52:15 -0500 |
parents | b09432b20a84 |
children | 7854d5752dd2 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2009-2021 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{x} =} gmres (@var{A}, @var{b}, @var{restart}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}, @dots{}) ## @deftypefnx {} {@var{x} =} gmres (@var{A}, @var{b}, @var{restart}, @var{tol}, @var{maxit}, @var{M}, [], @var{x0}, @dots{}) ## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@var{A}, @var{b}, @dots{}) ## Solve @code{A x = b} using the Preconditioned GMRES iterative method with ## restart, a.k.a. PGMRES(restart). ## ## 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{Afun} such that @code{Afun(x) = A * x}. Additional ## parameters to @code{Afun} are passed after @var{x0}. ## ## @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{tol} 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})) ## @leq{} @var{tol} * norm (inv (@var{M}) * @var{B})}. If @var{tol} is ## omitted or empty, then a tolerance of 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{M} is given as @code{M = 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{M}) * @var{A} * @var{x} = inv(@var{M}) * @var{b}} instead of ## @code{@var{A} * @var{x} = @var{b}}. ## ## @item @var{x0} is the initial guess, ## if not given or set to [], then the default value ## @code{zeros (size (@var{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{M} or ## @var{M1} or @var{M2}) which are passed to @code{gmres}. ## ## The 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 ## @item 0 : iteration converged to within the specified tolerance ## ## @item 1 : maximum number of iterations exceeded ## ## @item 2 : the preconditioner matrix is singular ## ## @item 3 : algorithm reached stagnation (the relative difference between two ## consecutive iterations is less than eps) ## @end table ## ## @item @var{relres} is the value of the relative preconditioned ## residual of the approximation @var{x}. ## ## @item @var{iter} is a vector containing the number of outer iterations and ## inner iterations performed to compute @var{x}. That is: ## ## @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 @leq{} @var{iter(1)} @leq{} @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 @leq{} @var{iter(2)} @leq{} @var{maxit}. ## @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 ## ## Reference: ## ## @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear ## Systems}, Second edition, 2003, SIAM ## ## @seealso{bicg, bicgstab, cgs, pcg, pcr, qmr, tfqmr} ## @end deftypefn function [x_min, flag, relres, it, resvec] = ... gmres (A, b, restart = [], tol = [], maxit = [], M1 = [], M2 = [], x0 = [], varargin) if (strcmp (class (A), "single") || strcmp (class (b), "single")) class_name = "single"; else class_name = "double"; endif [Afun, M1fun, M2fun] = __alltohandles__ (A, b, M1, M2, "gmres"); ## Check if the inputs are empty, and in case set them [tol, x0] = __default__input__ ({1e-06, zeros(size (b))}, tol, x0); empty_restart = isempty (restart); empty_maxit = isempty (maxit); size_b = rows (b); 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 gmres. \n Try to use a bigger tolerance."); 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; else max_iter_number = restart*maxit; 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 endif ## gmres: function handle case x_old = x_pr = x_min = x = x0; B = zeros (restart + 1, 1); V = zeros (rows (x), restart, class_name); 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 ## begin loop u = feval (Afun, x_old, varargin{:}); try warning("error", "Octave:singular-matrix", "local") prec_res = feval (M1fun, b - u, varargin{:}); # M1*(b-u) prec_res = feval (M2fun, prec_res, varargin{:}); presn = norm (prec_res, 2); resvec(1) = presn; z = feval (M1fun, b, varargin{:}); z = feval (M2fun, z, varargin{:}); 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 > tol * prec_b_norm) ## restart if (restart_it > restart) restart_it = 1; outer_it += 1; x_old = x; u = feval (Afun, x_old, varargin{:}); prec_res = feval (M1fun, b - u, varargin{:}); prec_res = feval (M2fun, prec_res, varargin{:}); presn = norm (prec_res, 2); B(1) = presn; H(:) = 0; V(:, 1) = prec_res / presn; endif ## basic iteration u = feval (Afun, V(:, restart_it), varargin{:}); tmp = feval (M1fun, u, varargin{:}); tmp = feval (M2fun, tmp, 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 if (flag == 2) resvec = norm (b); relres = 1; else resvec = resvec (1:iter); relres = resvec (iter) / prec_b_norm; endif if ((relres <= tol) && (flag == 1)) flag = 0; # Converged to solution within tolerance 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 ", tol); 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 ", tol); 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", tol); 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 ", tol); 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 ", tol); 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 ", tol); printf ("because it stagnates\n"); printf ("The iterated returned (number %d) ", it(2)); printf ("has relative residual %d \n", relres); endswitch 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) %!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 %!test %! ## Check that all type of inputs work %! A = toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])); %! b = sum (A, 2); %! 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); %! assert (x, A\b, 1e-9*norm (x, Inf)); %! [x, flag] = 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 = 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, [], []); %! assert (x, A\b, 1e-9*norm (x, Inf)); %! [x, flag] = 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 %! A = toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])) + ... %! 1i * toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])); %! b = sum (A, 2); %! [x, flag] = gmres(A, b, [], [], 5); %! assert (flag, 0); %! assert (x, ones (5, 1), -1e-06) %!test %! ## Maximum number of iteration reached %! A = hilb (100); %! b = sum (A, 2); %! [x, flag, relres, iter] = gmres (A, b, [], 1e-14); %! assert(flag, 1); %!test %! ## gmres recognizes that the preconditioner matrix is singular %! AA = 2 * eye (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) 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 = toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])); %! b = A * ones (5, 1); %! [x, flag, relres, iter] = gmres (A, b, [], [], [], [], [], ... %! ones (5, 1) + 1e-8); %! 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 = hilb (20); %! b = A * ones (20, 1); %! [x, flag, relres, iter, resvec] = gmres (A, b ,5, 1e-14); %! assert (iter, [4, 5]) %!test %! A = single (1); %! b = 1; %! [x, flag] = gmres (A, b); %! assert (class (x), "single") %!test %! A = 1; %! b = single (1); %! [x, flag] = gmres (A, b); %! assert (class (x), "single") %!test %! A = single (1); %! b = single (1); %! [x, flag] = gmres (A, b); %! assert (class (x), "single") %!test %!function y = Afun (x) %! A = toeplitz ([2, 1, 0, 0], [2, -1, 0, 0]); %! y = A * x; %!endfunction %! [x, flag] = gmres ("Afun", [1; 2; 2; 3]); %! assert (x, ones(4, 1), 1e-6) %!test # preconditioned residual %! A = toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])); %! b = sum (A, 2); %! M = magic (5); %! [x, flag, relres] = gmres (A, b, [], [], 2, M); %! assert (relres, norm (M \ (b - A * x)) / norm (M \ b), 8 * eps)