Mercurial > hg > octave-lojdl
changeset 12509:e742720c5e71
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 +++ b/doc/ChangeLog @@ -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 +++ b/doc/interpreter/linalg.txi @@ -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 +++ b/scripts/ChangeLog @@ -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.
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/gmres.m @@ -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 +++ b/src/ChangeLog @@ -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
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/mgorth.cc @@ -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); + +*/