octave_packages/m/sparse/pcr.m

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