1######################################################################## 2## 3## Copyright (C) 2004-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{x} =} pcr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m}, @var{x0}, @dots{}) 28## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} pcr (@dots{}) 29## 30## Solve the linear system of equations @code{@var{A} * @var{x} = @var{b}} by 31## means of the Preconditioned Conjugate Residuals iterative method. 32## 33## The input arguments are 34## 35## @itemize 36## @item 37## @var{A} can be either a square (preferably sparse) matrix or a function 38## handle, inline function or string containing the name of a function which 39## computes @code{@var{A} * @var{x}}. In principle @var{A} should be 40## symmetric and non-singular; if @code{pcr} finds @var{A} to be numerically 41## singular, you will get a warning message and the @var{flag} output 42## parameter will be set. 43## 44## @item 45## @var{b} is the right hand side vector. 46## 47## @item 48## @var{tol} is the required relative tolerance for the residual error, 49## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if 50## @code{norm (@var{b} - @var{A} * @var{x}) <= 51## @var{tol} * norm (@var{b} - @var{A} * @var{x0})}. 52## If @var{tol} is empty or is omitted, the function sets 53## @code{@var{tol} = 1e-6} by default. 54## 55## @item 56## @var{maxit} is the maximum allowable number of iterations; if @code{[]} is 57## supplied for @var{maxit}, or @code{pcr} has less arguments, a default 58## value equal to 20 is used. 59## 60## @item 61## @var{m} is the (left) preconditioning matrix, so that the iteration is 62## (theoretically) equivalent to solving by 63## @code{pcr} @code{@var{P} * @var{x} = @var{m} \ @var{b}}, with 64## @code{@var{P} = @var{m} \ @var{A}}. Note that a proper choice of the 65## preconditioner may dramatically improve the overall performance of the 66## method. Instead of matrix @var{m}, the user may pass a function which 67## returns the results of applying the inverse of @var{m} to a vector 68## (usually this is the preferred way of using the preconditioner). If 69## @code{[]} is supplied for @var{m}, or @var{m} is omitted, no 70## preconditioning is applied. 71## 72## @item 73## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the 74## function sets @var{x0} to a zero vector by default. 75## @end itemize 76## 77## The arguments which follow @var{x0} are treated as parameters, and passed 78## in a proper way to any of the functions (@var{A} or @var{m}) which are 79## passed to @code{pcr}. See the examples below for further details. 80## 81## The output arguments are 82## 83## @itemize 84## @item 85## @var{x} is the computed approximation to the solution of 86## @code{@var{A} * @var{x} = @var{b}}. 87## 88## @item 89## @var{flag} reports on the convergence. @code{@var{flag} = 0} means the 90## solution converged and the tolerance criterion given by @var{tol} is 91## satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit for the 92## iteration count was reached. @code{@var{flag} = 3} reports a @code{pcr} 93## breakdown, see [1] for details. 94## 95## @item 96## @var{relres} is the ratio of the final residual to its initial value, 97## measured in the Euclidean norm. 98## 99## @item 100## @var{iter} is the actual number of iterations performed. 101## 102## @item 103## @var{resvec} describes the convergence history of the method, so that 104## @code{@var{resvec} (i)} contains the Euclidean norms of the residual after 105## the (@var{i}-1)-th iteration, @code{@var{i} = 1,2, @dots{}, @var{iter}+1}. 106## @end itemize 107## 108## Let us consider a trivial problem with a diagonal matrix (we exploit the 109## sparsity of A) 110## 111## @example 112## @group 113## n = 10; 114## A = sparse (diag (1:n)); 115## b = rand (N, 1); 116## @end group 117## @end example 118## 119## @sc{Example 1:} Simplest use of @code{pcr} 120## 121## @example 122## x = pcr (A, b) 123## @end example 124## 125## @sc{Example 2:} @code{pcr} with a function which computes 126## @code{@var{A} * @var{x}}. 127## 128## @example 129## @group 130## function y = apply_a (x) 131## y = [1:10]' .* x; 132## endfunction 133## 134## x = pcr ("apply_a", b) 135## @end group 136## @end example 137## 138## @sc{Example 3:} Preconditioned iteration, with full diagnostics. The 139## preconditioner (quite strange, because even the original matrix 140## @var{A} is trivial) is defined as a function 141## 142## @example 143## @group 144## function y = apply_m (x) 145## k = floor (length (x) - 2); 146## y = x; 147## y(1:k) = x(1:k) ./ [1:k]'; 148## endfunction 149## 150## [x, flag, relres, iter, resvec] = ... 151## pcr (A, b, [], [], "apply_m") 152## semilogy ([1:iter+1], resvec); 153## @end group 154## @end example 155## 156## @sc{Example 4:} Finally, a preconditioner which depends on a 157## parameter @var{k}. 158## 159## @example 160## @group 161## function y = apply_m (x, varargin) 162## k = varargin@{1@}; 163## y = x; 164## y(1:k) = x(1:k) ./ [1:k]'; 165## endfunction 166## 167## [x, flag, relres, iter, resvec] = ... 168## pcr (A, b, [], [], "apply_m"', [], 3) 169## @end group 170## @end example 171## 172## Reference: 173## 174## @nospell{W. Hackbusch}, @cite{Iterative Solution of Large Sparse 175## Systems of Equations}, section 9.5.4; @nospell{Springer}, 1994 176## 177## @seealso{sparse, pcg} 178## @end deftypefn 179 180function [x, flag, relres, iter, resvec] = pcr (A, b, tol, maxit, m, x0, varargin) 181 182 breakdown = false; 183 184 if (nargin < 6 || isempty (x0)) 185 x = zeros (size (b)); 186 else 187 x = x0; 188 endif 189 190 if (nargin < 5) 191 m = []; 192 endif 193 194 if (nargin < 4 || isempty (maxit)) 195 maxit = 20; 196 endif 197 198 maxit += 2; 199 200 if (nargin < 3 || isempty (tol)) 201 tol = 1e-6; 202 endif 203 204 if (nargin < 2) 205 print_usage (); 206 endif 207 208 ## init 209 if (isnumeric (A)) # is A a matrix? 210 r = b - A*x; 211 else # then A should be a function! 212 r = b - feval (A, x, varargin{:}); 213 endif 214 215 if (isnumeric (m)) # is M a matrix? 216 if (isempty (m)) # if M is empty, use no precond 217 p = r; 218 else # otherwise, apply the precond 219 p = m \ r; 220 endif 221 else # then M should be a function! 222 p = feval (m, r, varargin{:}); 223 endif 224 225 iter = 2; 226 227 b_bot_old = 1; 228 q_old = p_old = s_old = zeros (size (x)); 229 230 if (isnumeric (A)) # is A a matrix? 231 q = A * p; 232 else # then A should be a function! 233 q = feval (A, p, varargin{:}); 234 endif 235 236 resvec(1) = abs (norm (r)); 237 238 ## iteration 239 while (resvec(iter-1) > tol*resvec(1) && iter < maxit) 240 241 if (isnumeric (m)) # is M a matrix? 242 if (isempty (m)) # if M is empty, use no precond 243 s = q; 244 else # otherwise, apply the precond 245 s = m \ q; 246 endif 247 else # then M should be a function! 248 s = feval (m, q, varargin{:}); 249 endif 250 b_top = r' * s; 251 b_bot = q' * s; 252 253 if (b_bot == 0.0) 254 breakdown = true; 255 break; 256 endif 257 lambda = b_top / b_bot; 258 259 x += lambda*p; 260 r -= lambda*q; 261 262 if (isnumeric (A)) # is A a matrix? 263 t = A*s; 264 else # then A should be a function! 265 t = feval (A, s, varargin{:}); 266 endif 267 268 alpha0 = (t'*s) / b_bot; 269 alpha1 = (t'*s_old) / b_bot_old; 270 271 p_temp = p; 272 q_temp = q; 273 274 p = s - alpha0*p - alpha1*p_old; 275 q = t - alpha0*q - alpha1*q_old; 276 277 s_old = s; 278 p_old = p_temp; 279 q_old = q_temp; 280 b_bot_old = b_bot; 281 282 resvec(iter) = abs (norm (r)); 283 iter += 1; 284 endwhile 285 286 flag = 0; 287 relres = resvec(iter-1) ./ resvec(1); 288 iter -= 2; 289 if (iter >= maxit-2) 290 flag = 1; 291 if (nargout < 2) 292 warning ("pcr: maximum number of iterations (%d) reached\n", iter); 293 warning ("pcr: the initial residual norm was reduced %g times\n", 294 1.0/relres); 295 endif 296 elseif (nargout < 2 && ! breakdown) 297 fprintf (stderr, "pcr: converged in %d iterations. \n", iter); 298 fprintf (stderr, "pcr: the initial residual norm was reduced %g times\n", 299 1.0 / relres); 300 endif 301 302 if (breakdown) 303 flag = 3; 304 if (nargout < 2) 305 warning ("pcr: breakdown occurred:\n"); 306 warning ("system matrix singular or preconditioner indefinite?\n"); 307 endif 308 endif 309 310endfunction 311 312 313%!demo 314%! ## Simplest usage of PCR (see also 'help pcr') 315%! 316%! N = 20; 317%! A = diag (linspace (-3.1,3,N)); b = rand (N,1); 318%! y = A \ b; # y is the true solution 319%! x = pcr (A,b); 320%! printf ("The solution relative error is %g\n", norm (x-y) / norm (y)); 321%! 322%! ## You shouldn't be afraid if PCR issues some warning messages in this 323%! ## example: watch out in the second example, why it takes N iterations 324%! ## of PCR to converge to (a very accurate, by the way) solution. 325 326%!demo 327%! ## Full output from PCR 328%! ## We use this output to plot the convergence history 329%! 330%! N = 20; 331%! A = diag (linspace (-3.1,30,N)); b = rand (N,1); 332%! X = A \ b; # X is the true solution 333%! [x, flag, relres, iter, resvec] = pcr (A,b); 334%! printf ("The solution relative error is %g\n", norm (x-X) / norm (X)); 335%! clf; 336%! title ("Convergence history"); 337%! xlabel ("Iteration"); ylabel ("log(||b-Ax||/||b||)"); 338%! semilogy ([0:iter], resvec/resvec(1), "o-g;relative residual;"); 339 340%!demo 341%! ## Full output from PCR 342%! ## We use indefinite matrix based on the Hilbert matrix, with one 343%! ## strongly negative eigenvalue 344%! ## Hilbert matrix is extremely ill conditioned, so is ours, 345%! ## and that's why PCR WILL have problems 346%! 347%! N = 10; 348%! A = hilb (N); A(1,1) = -A(1,1); b = rand (N,1); 349%! X = A \ b; # X is the true solution 350%! printf ("Condition number of A is %g\n", cond (A)); 351%! [x, flag, relres, iter, resvec] = pcr (A,b,[],200); 352%! if (flag == 3) 353%! printf ("PCR breakdown. System matrix is [close to] singular\n"); 354%! endif 355%! clf; 356%! title ("Convergence history"); 357%! xlabel ("Iteration"); ylabel ("log(||b-Ax||)"); 358%! semilogy ([0:iter], resvec, "o-g;absolute residual;"); 359 360%!demo 361%! ## Full output from PCR 362%! ## We use an indefinite matrix based on the 1-D Laplacian matrix for A, 363%! ## and here we have cond(A) = O(N^2) 364%! ## That's the reason we need some preconditioner; here we take 365%! ## a very simple and not powerful Jacobi preconditioner, 366%! ## which is the diagonal of A. 367%! 368%! ## Note that we use here indefinite preconditioners! 369%! 370%! N = 100; 371%! ## Form 1-D Laplacian matrix 372%! A = 2 * eye (N,N); 373%! A(2:(N+1):end) = -1; 374%! A((N+1):(N+1):end) = -1; 375%! 376%! A = [A, zeros(size(A)); zeros(size(A)), -A]; 377%! b = rand (2*N,1); 378%! X = A \ b; # X is the true solution 379%! maxit = 80; 380%! printf ("System condition number is %g\n", cond (A)); 381%! ## No preconditioner: the convergence is very slow! 382%! 383%! [x, flag, relres, iter, resvec] = pcr (A,b,[],maxit); 384%! clf; 385%! title ("Convergence history"); 386%! xlabel ("Iteration"); ylabel ("log(||b-Ax||)"); 387%! semilogy ([0:iter], resvec, "o-g;NO preconditioning: absolute residual;"); 388%! 389%! pause (1); 390%! ## Test Jacobi preconditioner: it will not help much!!! 391%! 392%! M = diag (diag (A)); # Jacobi preconditioner 393%! [x, flag, relres, iter, resvec] = pcr (A,b,[],maxit,M); 394%! hold on; 395%! semilogy ([0:iter],resvec,"o-r;JACOBI preconditioner: absolute residual;"); 396%! 397%! pause (1); 398%! ## Test nonoverlapping block Jacobi preconditioner: this one should give 399%! ## some convergence speedup! 400%! 401%! M = zeros (N,N); k = 4; 402%! for i=1:k:N # get k x k diagonal blocks of A 403%! M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1); 404%! endfor 405%! M = [M, zeros(size (M)); zeros(size(M)), -M]; 406%! [x, flag, relres, iter, resvec] = pcr (A,b,[],maxit,M); 407%! semilogy ([0:iter], resvec, "o-b;BLOCK JACOBI preconditioner: absolute residual;"); 408%! hold off; 409 410%!test 411%! ## solve small indefinite diagonal system 412%! 413%! N = 10; 414%! A = diag (linspace (-10.1,10,N)); b = ones (N,1); 415%! X = A \ b; # X is the true solution 416%! [x, flag] = pcr (A,b,[],N+1); 417%! assert (norm (x-X) / norm (X) < 1e-10); 418%! assert (flag, 0); 419 420%!test 421%! ## solve tridiagonal system, do not converge in default 20 iterations 422%! ## should perform max allowable default number of iterations 423%! 424%! N = 100; 425%! ## Form 1-D Laplacian matrix 426%! A = 2 * eye (N,N); 427%! A(2:(N+1):end) = -1; 428%! A((N+1):(N+1):end) = -1; 429%! b = ones (N,1); 430%! X = A \ b; # X is the true solution 431%! [x, flag, relres, iter, resvec] = pcr (A,b,1e-12); 432%! assert (flag, 1); 433%! assert (relres > 0.6); 434%! assert (iter, 20); 435 436%!test 437%! ## solve tridiagonal system with "perfect" preconditioner 438%! ## converges in one iteration 439%! 440%! N = 100; 441%! ## Form 1-D Laplacian matrix 442%! A = 2 * eye (N,N); 443%! A(2:(N+1):end) = -1; 444%! A((N+1):(N+1):end) = -1; 445%! b = ones (N,1); 446%! X = A \ b; # X is the true solution 447%! [x, flag, relres, iter] = pcr (A,b,[],[],A,b); 448%! assert (norm (x-X) / norm(X) < 1e-6); 449%! assert (relres < 1e-6); 450%! assert (flag, 0); 451%! assert (iter, 1); # should converge in one iteration 452