Mercurial > octave-nkf

changeset 12509:e742720c5e71

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add mgorth and gmres functions from Octave Forge
author Carlo de Falco <kingcrimson@tiscali.it>
date Wed, 16 Mar 2011 13:36:27 -0700
parents 919cadf334f8
children a1b2da4967ac
files doc/ChangeLog doc/interpreter/linalg.txi scripts/ChangeLog scripts/linear-algebra/gmres.m src/ChangeLog src/DLD-FUNCTIONS/mgorth.cc src/DLD-FUNCTIONS/module-files
diffstat 7 files changed, 388 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/doc/ChangeLog	Tue Mar 15 22:45:31 2011 -0700
+++ b/doc/ChangeLog	Wed Mar 16 13:36:27 2011 -0700
@@ -49,6 +49,11 @@
 	* interpreter/arith.txi, interpreter/io.txi, interpreter/oop.txi:
 	Remove functions which have no DOCSTRING entries.
 
+2011-02-10 Carlo de Falco  <kingcrimson@tiscali.it>
+
+	* interpreter/linalg.txi: Added gmres to the specialized solvers.
+
+
 2011-02-06  John W. Eaton  <jwe@octave.org>
 
 	* interpreter/contributors.in: Add Fotios Kasolis to the list.
--- a/doc/interpreter/linalg.txi	Tue Mar 15 22:45:31 2011 -0700
+++ b/doc/interpreter/linalg.txi	Wed Mar 16 13:36:27 2011 -0700
@@ -106,6 +106,8 @@
 
 @DOCSTRING(orth)
 
+@DOCSTRING(mgorth)
+
 @DOCSTRING(pinv)
 
 @DOCSTRING(rank)
@@ -188,3 +190,5 @@
 @DOCSTRING(bicgstab)
 
 @DOCSTRING(cgs)
+
+@DOCSTRING(gmres)
--- a/scripts/ChangeLog	Tue Mar 15 22:45:31 2011 -0700
+++ b/scripts/ChangeLog	Wed Mar 16 13:36:27 2011 -0700
@@ -192,6 +192,11 @@
 
 	* statistics/base/mean.m: Also accept logical values.
 
+2011-02-10 Carlo de Falco  <kingcrimson@tiscali.it>
+
+	* linear-algebra/gmres.m: New file implementing the GMRES
+	iterative method for solving linear systems.
+
 2011-02-08  Ben Abbott  <bpabbott@mac.com>
 
 	* plot/__go_draw_axes__.m: Properly set fontspec for legends.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/gmres.m	Wed Mar 16 13:36:27 2011 -0700
@@ -0,0 +1,214 @@
+## Copyright (C) 2009-2011 Carlo de Falco
+##  
+## 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} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## @deftypefnx {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{P})
+## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@dots{})
+## Solve @code{A x = b} using the Preconditioned GMRES iterative method
+## with restart, a.k.a. PGMRES(m).
+##
+## @itemize @minus
+## @item @var{rtol} is the relative tolerance, 
+## if not given or set to [] the default value 1e-6 is used.
+##
+## @item @var{maxit} is the maximum number of outer iterations,
+## if not given or set to [] the default value 
+## @code{min (10, numel (b) / restart)} is used.
+##
+## @item @var{x0} is the initial guess,
+## if not given or set to [] the default value @code{zeros(size (b))} is used. 
+##
+## @item @var{m} is the restart parameter,
+## if not given or set to [] the default value @code{numel (b)} is used.
+## @end itemize
+##
+## Argument @var{A} can be passed as a matrix, function handle, or 
+## inline function @code{f} such that @code{f(x) = A*x}.
+##
+## The preconditioner @var{P} is given as @code{P = 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}.
+##
+## Besides the vector @var{x}, additional outputs are:
+##
+## @itemize @minus
+## @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 : unused, but skipped for compatibility
+##
+##   @item 3 : algorithm reached stagnation
+## @end table
+##
+## @item @var{relres} is the final value of the relative residual.
+##
+## @item @var{iter} is a vector containing the number of outer iterations and
+## total iterations performed. 
+##
+## @item @var{resvec} is a vector containing the relative residual at each
+## iteration.
+## @end itemize
+##
+## @seealso{pcg, cgs, bigcstab}
+## @end deftypefn
+
+function [x, flag, prec_res_norm, itcnt] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
+
+  if (nargin < 2 || nargin > 8)
+    print_usage ();    
+  end
+      
+  if (ischar (A))
+    Ax = str2func (A);
+  elseif (ismatrix (A))
+    Ax = @(x) A*x;
+  elseif (isa (A, "function_handle"))
+    Ax = A;
+  else
+    error ("gmres: A must be a function or matrix");
+  endif
+
+  if (nargin < 3 || isempty (restart))
+    restart = rows (b);
+  endif
+
+  if (nargin < 4 || isempty (rtol))
+    rtol = 1e-6;
+  endif
+
+  if (nargin < 5 || isempty (maxit))
+    maxit = min (rows (b)/restart, 10);
+  endif
+
+  if (nargin < 6 || isempty (M1))
+    M1m1x = @(x) x;
+  elseif (ischar (M1))
+    M1m1x = str2func (M1);
+  elseif (ismatrix (M1))
+    M1m1x = @(x) M1 \ x;
+  elseif (isa (M1, "function_handle"))
+    M1m1x = M1;
+  else
+    error ("gmres: preconditioner M1 must be a function or matrix");
+  endif
+  
+  if (nargin < 7 || isempty (M2))
+    M2m1x = @(x) x;
+  elseif (ischar (M2))
+    M2m1x = str2func (M2);
+  elseif (ismatrix (M2))
+    M2m1x = @(x) M2 \ x;
+  elseif (isa (M2, "function_handle"))
+    M2m1x = M2;
+  else
+    error ("gmres: preconditioner M2 must be a function or matrix");
+  endif
+
+  Pm1x = @(x) M2m1x (M1m1x (x));
+
+  if (nargin < 8 || isempty (x0))
+    x0 = zeros (size (b));
+  endif
+
+  x_old = x0; 
+  x = x_old;
+  prec_res = Pm1x (b - Ax (x_old));
+  prec_res_norm = norm (prec_res, 2);
+  
+  B = zeros (restart + 1, 1);
+  V = zeros (rows (x), restart);
+  H = zeros (restart + 1, restart);
+
+  ## begin loop
+  iter = 1;
+  restart_it  = restart + 1; 
+  resids      = zeros (maxit, 1);
+  resids(1)   = prec_res_norm;
+  prec_b_norm = norm (Pm1x (b), 2);
+  flag        = 1;
+
+  while ((iter <= maxit * restart) && (prec_res_norm > rtol * prec_b_norm))
+  
+    ## restart
+    if (restart_it > restart)
+      restart_it = 1;
+      x_old = x;	      
+      prec_res = Pm1x (b - Ax (x_old));
+      prec_res_norm = norm (prec_res, 2);
+      B(1) = prec_res_norm;
+      H(:) = 0;
+      V(:, 1) = prec_res / prec_res_norm;
+    endif  
+    
+    ## basic iteration
+    tmp = Pm1x (Ax (V(:, restart_it)));
+    [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);
+    prec_res_norm = norm (little_res, 2);
+	      
+    x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
+    
+    resids(iter) = prec_res_norm ;
+    if (norm (x - x_old, inf) <= eps)
+      flag = 3;
+      break
+    endif
+
+    restart_it++ ; iter++;
+  endwhile
+
+  if (prec_res_norm > rtol * prec_b_norm)
+    flag = 0;
+  endif
+
+  resids = resids(1:iter-1);
+  itcnt = [floor(maxit/restart), rem(maxit, restart)];
+endfunction
+
+
+%!shared A, b, dim
+%! dim = 100;
+%!test
+%! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
+%! b = ones(dim, 1); 
+%! x = gmres (A, b, 10, 1e-10, dim, @(x) x./diag(A), [],  b);
+%! assert(x, A\b, 1e-9*norm(x,inf));
+%!
+%!test
+%! x = gmres (A, b, dim, 1e-10, 1e4, @(x) diag(diag(A))\x, [],  b);  
+%! assert(x, A\b, 1e-7*norm(x,inf));
+%!
+%!test
+%! 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, resids] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag(A), [],  []);
+%! assert(x, A\b, 1e-9*norm(x,inf))
+%! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag(diag(A))\x, [],  []);
+%! assert(x, A\b, 1e-9*norm(x,inf));
+%!test
+%! x =  gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [],  []);
+%! assert(x, A\b, 1e-7*norm(x,inf));
--- a/src/ChangeLog	Tue Mar 15 22:45:31 2011 -0700
+++ b/src/ChangeLog	Wed Mar 16 13:36:27 2011 -0700
@@ -185,6 +185,11 @@
 	* load-path.cc (strip_trailing_separators): Declare K as size_t
 	rather than octave_idx_type.
 
+2011-02-10 Carlo de Falco  <kingcrimson@tiscali.it>
+
+	* DLD-FUNCTIONS/mgorth.cc: New file implementing modified
+	Gram-Schmidt orthogonalization.
+
 2011-02-10  John W. Eaton  <jwe@octave.org>
 
 	* DLD-FUNCTIONS/regexp.cc (octregexp_list): Avoid deprecated
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/mgorth.cc	Wed Mar 16 13:36:27 2011 -0700
@@ -0,0 +1,154 @@
+/*
+
+Copyright (C) 2009-2011 Carlo de Falco
+Copyright (C) 2010 VZLU Prague
+
+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/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "oct-norm.h"
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+
+template <class ColumnVector, class Matrix, class RowVector>
+static void 
+do_mgorth (ColumnVector& x, const Matrix& V, RowVector& h) 
+{
+  octave_idx_type Vc = V.columns ();
+  h = RowVector (Vc + 1);
+  for (octave_idx_type j = 0; j < Vc; j++)
+    {
+      ColumnVector Vcj = V.column (j);
+      h(j) = RowVector (Vcj.hermitian ()) * x;
+      x -= h(j) * Vcj;
+    }
+
+  h(Vc) = xnorm (x);
+  if (real (h(Vc)) > 0)
+    x = x / h(Vc);
+}
+
+DEFUN_DLD (mgorth, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{y}, @var{h}] =} mgorth (@var{x}, @var{v})\n\
+Orthogonalize a given column vector @var{x} with respect to a given\n\
+orthonormal basis @var{v} using a modified Gram-Schmidt orthogonalization.  \n\
+On exit, @var{y} is a unit vector such that:\n\
+\n\
+@example\n\
+@group\n\
+  norm (@var{y}) = 1\n\
+  @var{v}' * @var{y} = 0\n\
+  @var{x} = @var{h}*[@var{v}, @var{y}]\n\
+@end group\n\
+@end example\n\
+\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length();
+
+  if (nargin != 2 || nargout > 2)
+  {
+    print_usage (); 
+    return retval;
+  }
+
+  octave_value arg_x = args(0);
+  octave_value arg_v = args(1);
+
+  if (arg_v.ndims () != 2 || arg_x.ndims () != 2 || arg_x.columns () != 1
+      || arg_v.rows () != arg_x.rows ())
+    {
+      error ("mgorth: V should me a matrix, and X a column vector with"
+             " the same number of rows as V.");
+      return retval;
+    }
+
+  if (! arg_x.is_numeric_type () && ! arg_v.is_numeric_type ())
+    {
+      error ("mgorth: X and V must be numeric");
+    }
+ 
+  bool iscomplex = (arg_x.is_complex_type () || arg_v.is_complex_type ());
+  if (arg_x.is_single_type () || arg_v.is_single_type ())
+    {
+      if (iscomplex)
+        {
+          FloatComplexColumnVector x = arg_x.float_complex_column_vector_value ();
+          FloatComplexMatrix V = arg_v.float_complex_matrix_value ();
+          FloatComplexRowVector h;
+          do_mgorth (x, V, h);
+          retval(1) = h;
+          retval(0) = x;
+        }
+      else
+        {
+          FloatColumnVector x = arg_x.float_column_vector_value ();
+          FloatMatrix V = arg_v.float_matrix_value ();
+          FloatRowVector h;
+          do_mgorth (x, V, h);
+          retval(1) = h;
+          retval(0) = x;
+        }
+    }
+  else
+    {
+      if (iscomplex)
+        {
+          ComplexColumnVector x = arg_x.complex_column_vector_value ();
+          ComplexMatrix V = arg_v.complex_matrix_value ();
+          ComplexRowVector h;
+          do_mgorth (x, V, h);
+          retval(1) = h;
+          retval(0) = x;
+        }
+      else
+        {
+          ColumnVector x = arg_x.column_vector_value ();
+          Matrix V = arg_v.matrix_value ();
+          RowVector h;
+          do_mgorth (x, V, h);
+          retval(1) = h;
+          retval(0) = x;
+        }
+    }
+
+  return retval;
+}
+
+/*
+  
+%!test
+%! for ii=1:100; assert (abs (mgorth (randn (5, 1), eye (5, 4))), [0 0 0 0 1]', eps); endfor
+
+%!test
+%! a = hilb (5);
+%! a(:, 1) /= norm (a(:, 1));
+%! for ii = 1:5
+%!   a(:, ii) = mgorth (a(:, ii), a(:, 1:ii-1));
+%! endfor
+%! assert (a' * a, eye (5), 1e10);
+
+*/
--- a/src/DLD-FUNCTIONS/module-files	Tue Mar 15 22:45:31 2011 -0700
+++ b/src/DLD-FUNCTIONS/module-files	Wed Mar 16 13:36:27 2011 -0700
@@ -54,6 +54,7 @@
 matrix_type.cc
 max.cc
 md5sum.cc
+mgorth.cc
 onCleanup.cc
 pinv.cc
 qr.cc