# HG changeset patch # User Rik # Date 1382051105 25200 # Node ID f6fded8395133cc87b01bbe4d616d294bd6a66cf # Parent 4264c78951ec6beb990447da801480350f725940 pcg.m: Change success criterion to match reference text (bug #40057) * scripts/sparse/pcg.m: Change success criteria to be "norm (residual) <= tol * norm (b)" in accordance with reference text and Matlab. Redo docstring. Use Octave coding conventions regarding parentheses. diff -r 4264c78951ec -r f6fded839513 scripts/sparse/pcg.m --- a/scripts/sparse/pcg.m Tue Oct 08 13:38:37 2013 +0200 +++ b/scripts/sparse/pcg.m Thu Oct 17 16:05:05 2013 -0700 @@ -20,70 +20,68 @@ ## @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{}) ## -## Solve the linear system of equations @code{@var{A} * @var{x} = @var{b}} -## by means of the Preconditioned Conjugate Gradient iterative -## method. The input arguments are +## 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 ## ## @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 @code{@var{A} * @var{x}}. In principle -## @var{A} should be symmetric and positive definite; if @code{pcg} -## finds @var{A} to not be positive definite, you will get a warning -## message and the @var{flag} output parameter will be set. +## @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}}}. 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. +## @var{b} is the right-hand side vector. ## ## @item ## @var{tol} is the required relative tolerance for the residual error, -## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if -## @code{norm (@var{b} - @var{A} * @var{x}) <= -## @var{tol} * norm (@var{b} - @var{A} * @var{x0})}. -## If @var{tol} is empty or is omitted, the function sets -## @code{@var{tol} = 1e-6} by default. +## @w{@code{@var{b} - @var{A} * @var{x}}}. The iteration stops if +## @w{@code{norm (@var{b} - @var{A} * @var{x})} @leq{} +## @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 allowable number of iterations; if -## @code{[]} is supplied for @code{maxit}, or @code{pcg} has less -## arguments, a default value equal to 20 is used. +## @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} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that ## the iteration is (theoretically) equivalent to solving by @code{pcg} -## @code{@var{P} * -## @var{x} = @var{m} \ @var{b}}, with @code{@var{P} = @var{m} \ @var{A}}. +## @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 (usually this is the preferred way of using the preconditioner). -## If @code{[]} is supplied for @var{m1}, or @var{m1} is omitted, no -## preconditioning is applied. If @var{m2} is omitted, @var{m} = @var{m1} -## will be used as 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 empty or omitted, the +## @var{x0} is the initial guess. If @var{x0} is omitted or empty then the ## function sets @var{x0} to a zero vector by default. ## @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}) -## which are passed to @code{pcg}. See the examples below for further -## details. The output arguments are +## 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}) 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 -## @code{@var{A} * @var{x} = @var{b}}. +## @w{@code{@var{A} * @var{x} = @var{b}}}. ## ## @item -## @var{flag} reports on the convergence. @code{@var{flag} = 0} means -## the solution converged and the tolerance criterion given by @var{tol} -## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit -## for the iteration count was reached. @code{@var{flag} = 3} reports that -## the (preconditioned) matrix was found not positive definite. +## @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, @@ -94,29 +92,28 @@ ## ## @item ## @var{resvec} describes the convergence history of the method. -## @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{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{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)} and largest @code{@var{eigest} (2)} eigenvalues of the -## preconditioned matrix @code{@var{P} = @var{m} \ @var{A}}. In -## particular, if no preconditioning is used, the estimates for the -## extreme eigenvalues of @var{A} are returned. @code{@var{eigest} (1)} -## is an overestimate and @code{@var{eigest} (2)} is an underestimate, -## so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound -## for @code{cond (@var{P}, 2)}, which nevertheless in the limit should +## @var{eigest} returns the estimate for the smallest @code{@var{eigest}(1)} +## and largest @code{@var{eigest}(2)} eigenvalues of the preconditioned matrix +## @w{@code{@var{P} = @var{m} \ @var{A}}}. In particular, if no +## preconditioning is used, the estimates for the extreme eigenvalues of +## @var{A} are returned. @code{@var{eigest}(1)} is an overestimate and +## @code{@var{eigest}(2)} is an underestimate, so that +## @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. +## 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 diagonal matrix (we exploit the @@ -134,7 +131,7 @@ ## @sc{Example 1:} Simplest use of @code{pcg} ## ## @example -## x = pcg (A,b) +## x = pcg (A, b) ## @end example ## ## @sc{Example 2:} @code{pcg} with a function which computes @@ -270,11 +267,12 @@ r = b - feval (A, x, varargin{:}); endif + b_norm = norm (b); resvec(1,1) = norm (r); alpha = 1; iter = 2; - while (resvec (iter-1,1) > tol * resvec (1,1) && iter < maxit) + while (resvec(iter-1,1) > tol * b_norm && iter < maxit) if (exist_m1) if (isnumeric (m1)) y = m1 \ r; @@ -294,7 +292,7 @@ z = y; endif tau = z' * r; - resvec (iter-1,2) = sqrt (tau); + resvec(iter-1,2) = sqrt (tau); beta = tau / oldtau; oldtau = tau; p = z + beta * p; @@ -320,7 +318,7 @@ ## 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); + resvec(iter,1) = norm (r); iter++; endwhile @@ -360,13 +358,13 @@ z = y; endif - resvec (iter-1,2) = sqrt (r' * z); + resvec(iter-1,2) = sqrt (r' * z); else resvec = resvec(:,1); endif flag = 0; - relres = resvec (iter-1,1) ./ resvec(1,1); + relres = resvec(iter-1,1) ./ resvec(1,1); iter -= 2; if (iter >= maxit - 2) flag = 1;