X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fm%2Fsparse%2Fpcr.m;fp=octave_packages%2Fm%2Fsparse%2Fpcr.m;h=1edb73ff7dbc180a50771d2e6a6b96d64553a28c;hb=1c0469ada9531828709108a4882a751d2816994a;hp=0000000000000000000000000000000000000000;hpb=63de9f36673d49121015e3695f2c336ea92bc278;p=CreaPhase.git

diff --git a/octave_packages/m/sparse/pcr.m b/octave_packages/m/sparse/pcr.m
new file mode 100644
index 0000000..1edb73f
--- /dev/null
+++ b/octave_packages/m/sparse/pcr.m
@@ -0,0 +1,432 @@
+## Copyright (C) 2004-2012 Piotr Krzyzanowski
+##
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{x} =} pcr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m}, @var{x0}, @dots{})
+## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} pcr (@dots{})
+##
+## Solve the linear system of equations @code{@var{A} * @var{x} = @var{b}}
+## by means of the Preconditioned Conjugate Residuals 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 non-singular; if @code{pcr}
+## finds @var{A} to be numerically singular, you will get a warning
+## message and the @var{flag} output parameter will be set.
+##
+## @item
+## @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.
+##
+## @item
+## @var{maxit} is the maximum allowable number of iterations; if
+## @code{[]} is supplied for @code{maxit}, or @code{pcr} has less
+## arguments, a default value equal to 20 is used.
+##
+## @item
+## @var{m} is the (left) preconditioning matrix, so that the iteration is
+## (theoretically) equivalent to solving by @code{pcr} @code{@var{P} *
+## @var{x} = @var{m} \ @var{b}}, with @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 matrix
+## @var{m}, the user may pass a function which returns the results of
+## applying the inverse of @var{m} to a vector (usually this is the
+## preferred way of using the preconditioner).  If @code{[]} is supplied
+## for @var{m}, or @var{m} is omitted, no preconditioning is applied.
+##
+## @item
+## @var{x0} is the initial guess.  If @var{x0} is empty or omitted, 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{pcr}.  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}}.
+##
+## @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 t
+## @code{pcr} breakdown, see [1] for details.
+##
+## @item
+## @var{relres} is the ratio of the final residual to its initial value,
+## measured in the Euclidean norm.
+##
+## @item
+## @var{iter} is the actual number of iterations performed.
+##
+## @item
+## @var{resvec} describes the convergence history of the method,
+## so that @code{@var{resvec} (i)} contains the Euclidean norms of the
+## residual after the (@var{i}-1)-th iteration, @code{@var{i} =
+## 1,2, @dots{}, @var{iter}+1}.
+## @end itemize
+##
+## Let us consider a trivial problem with a diagonal matrix (we exploit the
+## sparsity of A)
+##
+## @example
+## @group
+## n = 10;
+## A = sparse (diag (1:n));
+## b = rand (N, 1);
+## @end group
+## @end example
+##
+## @sc{Example 1:} Simplest use of @code{pcr}
+##
+## @example
+## x = pcr (A, b)
+## @end example
+##
+## @sc{Example 2:} @code{pcr} with a function which computes
+## @code{@var{A} * @var{x}}.
+##
+## @example
+## @group
+## function y = apply_a (x)
+##   y = [1:10]' .* x;
+## endfunction
+##
+## x = pcr ("apply_a", b)
+## @end group
+## @end example
+##
+## @sc{Example 3:}  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 = apply_m (x)
+##   k = floor (length (x) - 2);
+##   y = x;
+##   y(1:k) = x(1:k) ./ [1:k]';
+## endfunction
+##
+## [x, flag, relres, iter, resvec] = ...
+##                    pcr (A, b, [], [], "apply_m")
+## semilogy ([1:iter+1], resvec);
+## @end group
+## @end example
+##
+## @sc{Example 4:} Finally, a preconditioner which depends on a
+## parameter @var{k}.
+##
+## @example
+## @group
+## function y = apply_m (x, varargin)
+##   k = varargin@{1@};
+##   y = x;
+##   y(1:k) = x(1:k) ./ [1:k]';
+## endfunction
+##
+## [x, flag, relres, iter, resvec] = ...
+##                    pcr (A, b, [], [], "apply_m"', [], 3)
+## @end group
+## @end example
+##
+## References:
+##
+##      [1] W. Hackbusch, @cite{Iterative Solution of Large Sparse Systems of
+##      Equations}, section 9.5.4; Springer, 1994
+##
+## @seealso{sparse, pcg}
+## @end deftypefn
+
+## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
+
+function [x, flag, relres, iter, resvec] = pcr (A, b, tol, maxit, m, x0, varargin)
+
+  breakdown = false;
+
+  if (nargin < 6 || isempty (x0))
+    x = zeros (size (b));
+  else
+    x = x0;
+  endif
+
+  if (nargin < 5)
+    m = [];
+  endif
+
+  if (nargin < 4 || isempty (maxit))
+    maxit = 20;
+  endif
+
+  maxit += 2;
+
+  if (nargin < 3 || isempty (tol))
+    tol = 1e-6;
+  endif
+
+  if (nargin < 2)
+    print_usage ();
+  endif
+
+  ##  init
+  if (isnumeric (A))            # is A a matrix?
+    r = b - A*x;
+  else                          # then A should be a function!
+    r = b - feval (A, x, varargin{:});
+  endif
+
+  if (isnumeric (m))            # is M a matrix?
+    if (isempty (m))            # if M is empty, use no precond
+      p = r;
+    else                        # otherwise, apply the precond
+      p = m \ r;
+    endif
+  else                          # then M should be a function!
+    p = feval (m, r, varargin{:});
+  endif
+
+  iter = 2;
+
+  b_bot_old = 1;
+  q_old = p_old = s_old = zeros (size (x));
+
+  if (isnumeric (A))            # is A a matrix?
+    q = A * p;
+  else                          # then A should be a function!
+    q = feval (A, p, varargin{:});
+  endif
+
+  resvec(1) = abs (norm (r));
+
+  ## iteration
+  while (resvec(iter-1) > tol*resvec(1) && iter < maxit)
+
+    if (isnumeric (m))          # is M a matrix?
+      if (isempty (m))          # if M is empty, use no precond
+        s = q;
+      else                      # otherwise, apply the precond
+        s = m \ q;
+      endif
+    else                        # then M should be a function!
+      s = feval (m, q, varargin{:});
+    endif
+    b_top = r' * s;
+    b_bot = q' * s;
+
+    if (b_bot == 0.0)
+      breakdown = true;
+      break;
+    endif
+    lambda = b_top / b_bot;
+
+    x += lambda*p;
+    r -= lambda*q;
+
+    if (isnumeric(A))           # is A a matrix?
+      t = A*s;
+    else                        # then A should be a function!
+      t = feval (A, s, varargin{:});
+    endif
+
+    alpha0 = (t'*s) / b_bot;
+    alpha1 = (t'*s_old) / b_bot_old;
+
+    p_temp = p;
+    q_temp = q;
+
+    p = s - alpha0*p - alpha1*p_old;
+    q = t - alpha0*q - alpha1*q_old;
+
+    s_old = s;
+    p_old = p_temp;
+    q_old = q_temp;
+    b_bot_old = b_bot;
+
+    resvec(iter) = abs (norm (r));
+    iter++;
+  endwhile
+
+  flag = 0;
+  relres = resvec(iter-1) ./ resvec(1);
+  iter -= 2;
+  if (iter >= maxit-2)
+    flag = 1;
+    if (nargout < 2)
+      warning ("pcr: maximum number of iterations (%d) reached\n", iter);
+      warning ("the initial residual norm was reduced %g times.\n", 1.0/relres);
+    endif
+  elseif (nargout < 2 && ! breakdown)
+    fprintf (stderr, "pcr: converged in %d iterations. \n", iter);
+    fprintf (stderr, "the initial residual norm was reduced %g times.\n",
+             1.0 / relres);
+  endif
+
+  if (breakdown)
+    flag = 3;
+    if (nargout < 2)
+      warning ("pcr: breakdown occurred:\n");
+      warning ("system matrix singular or preconditioner indefinite?\n");
+    endif
+  endif
+
+endfunction
+
+%!demo
+%!
+%!      # Simplest usage of PCR (see also 'help pcr')
+%!
+%!      N = 20;
+%!      A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
+%!      x = pcr(A,b);
+%!      printf('The solution relative error is %g\n', norm(x-y)/norm(y));
+%!
+%!      # You shouldn't be afraid if PCR issues some warning messages in this
+%!      # example: watch out in the second example, why it takes N iterations
+%!      # of PCR to converge to (a very accurate, by the way) solution
+%!demo
+%!
+%!      # Full output from PCR
+%!      # We use this output to plot the convergence history
+%!
+%!      N = 20;
+%!      A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
+%!      [x, flag, relres, iter, resvec] = pcr(A,b);
+%!      printf('The solution relative error is %g\n', norm(x-X)/norm(X));
+%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
+%!      semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');
+%!demo
+%!
+%!      # Full output from PCR
+%!      # We use indefinite matrix based on the Hilbert matrix, with one
+%!      # strongly negative eigenvalue
+%!      # Hilbert matrix is extremely ill conditioned, so is ours,
+%!      # and that's why PCR WILL have problems
+%!
+%!      N = 10;
+%!      A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
+%!      printf('Condition number of A is   %g\n', cond(A));
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],200);
+%!      if (flag == 3)
+%!        printf('PCR breakdown. System matrix is [close to] singular\n');
+%!      end
+%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!      semilogy([0:iter],resvec,'o-g;absolute residual;');
+%!demo
+%!
+%!      # Full output from PCR
+%!      # We use an indefinite matrix based on the 1-D Laplacian matrix for A,
+%!      # and here we have cond(A) = O(N^2)
+%!      # 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
+%!
+%!      # Note that we use here indefinite preconditioners!
+%!
+%!      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
+%!      A = [A, zeros(size(A)); zeros(size(A)), -A];
+%!      b = rand(2*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] = pcr(A,b,[],maxit);
+%!      title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
+%!      semilogy([0:iter],resvec,'o-g;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] = pcr(A,b,[],maxit,M);
+%!      hold on;
+%!      semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');
+%!
+%!      pause(1);
+%!      # Test nonoverlapping block Jacobi preconditioner: this one should give
+%!      # some convergence speedup!
+%!
+%!      M = zeros(N,N);k=4;
+%!      for i=1:k:N # get k x k diagonal blocks of A
+%!              M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
+%!      endfor
+%!      M = [M, zeros(size(M)); zeros(size(M)), -M];
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
+%!      semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
+%!      hold off;
+%!test
+%!
+%!      #solve small indefinite diagonal system
+%!
+%!      N = 10;
+%!      A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution
+%!      [x, flag] = pcr(A,b,[],N+1);
+%!      assert(norm(x-X)/norm(X)<1e-10);
+%!      assert(flag,0);
+%!
+%!test
+%!
+%!      #solve tridiagonal system, do not converge in default 20 iterations
+%!      #should perform max allowable default number of iterations
+%!
+%!      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 = ones(N,1); X = A\b; #X is the true solution
+%!      [x, flag, relres, iter, resvec] = pcr(A,b,1e-12);
+%!      assert(flag,1);
+%!      assert(relres>0.6);
+%!      assert(iter,20);
+%!
+%!test
+%!
+%!      #solve tridiagonal system with 'prefect' preconditioner
+%!      #converges in one iteration
+%!
+%!      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 = ones(N,1); X = A\b; #X is the true solution
+%!      [x, flag, relres, iter] = pcr(A,b,[],[],A,b);
+%!      assert(norm(x-X)/norm(X)<1e-6);
+%!      assert(relres<1e-6);
+%!      assert(flag,0);
+%!      assert(iter,1); #should converge in one iteration
+%!