Mercurial > octave

annotate scripts/sparse/qmr.m @ 19597:db92e7e28e1f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from most source files * NEWS, doc/interpreter/contributors.in, doc/interpreter/func.txi, doc/interpreter/genpropdoc.m, doc/interpreter/octave_logo.eps, doc/interpreter/plot.txi, doc/interpreter/stmt.txi, examples/data/Makefile.am, libinterp/corefcn/data.cc, libinterp/corefcn/debug.cc, libinterp/corefcn/error.cc, libinterp/corefcn/file-io.cc, libinterp/corefcn/gl-render.cc, libinterp/corefcn/graphics.cc, libinterp/corefcn/graphics.in.h, libinterp/corefcn/load-path.cc, libinterp/corefcn/pr-output.cc, libinterp/corefcn/pt-jit.cc, libinterp/corefcn/strfind.cc, libinterp/corefcn/toplev.cc, libinterp/corefcn/toplev.h, libinterp/corefcn/urlwrite.cc, libinterp/corefcn/variables.cc, libinterp/octave-value/ov-classdef.cc, libinterp/octave-value/ov-classdef.h, libinterp/octave.cc, libinterp/parse-tree/lex.h, libinterp/parse-tree/oct-parse.in.yy, libinterp/parse-tree/pt-classdef.h, liboctave/system/file-ops.cc, liboctave/system/oct-env.cc, m4/acinclude.m4, scripts/deprecated/finite.m, scripts/deprecated/fmod.m, scripts/deprecated/fnmatch.m, scripts/deprecated/luinc.m, scripts/deprecated/octave_tmp_file_name.m, scripts/deprecated/syl.m, scripts/deprecated/usage.m, scripts/general/inputParser.m, scripts/general/interp1.m, scripts/general/interp2.m, scripts/general/interp3.m, scripts/general/isequal.m, scripts/general/private/__isequal__.m, scripts/geometry/voronoi.m, scripts/image/image.m, scripts/image/imshow.m, scripts/image/ind2rgb.m, scripts/linear-algebra/bandwidth.m, scripts/linear-algebra/isbanded.m, scripts/miscellaneous/bzip2.m, scripts/miscellaneous/cast.m, scripts/miscellaneous/copyfile.m, scripts/miscellaneous/delete.m, scripts/miscellaneous/fullfile.m, scripts/miscellaneous/getappdata.m, scripts/miscellaneous/gunzip.m, scripts/miscellaneous/isappdata.m, scripts/miscellaneous/ls.m, scripts/miscellaneous/mex.m, scripts/miscellaneous/movefile.m, scripts/miscellaneous/orderfields.m, scripts/miscellaneous/recycle.m, scripts/miscellaneous/rmappdata.m, scripts/miscellaneous/setfield.m, scripts/miscellaneous/symvar.m, scripts/miscellaneous/tar.m, scripts/miscellaneous/tmpnam.m, scripts/miscellaneous/unpack.m, scripts/miscellaneous/ver.m, scripts/miscellaneous/what.m, scripts/miscellaneous/xor.m, scripts/miscellaneous/zip.m, scripts/optimization/fminbnd.m, scripts/optimization/sqp.m, scripts/path/private/getsavepath.m, scripts/path/savepath.m, scripts/pkg/pkg.m, scripts/pkg/private/installed_packages.m, scripts/plot/draw/plotyy.m, scripts/plot/draw/polar.m, scripts/plot/draw/private/__quiver__.m, scripts/plot/draw/private/__scatter__.m, scripts/plot/draw/private/__stem__.m, scripts/plot/draw/surface.m, scripts/plot/draw/surfnorm.m, scripts/plot/util/copyobj.m, scripts/plot/util/hgload.m, scripts/plot/util/hgsave.m, scripts/plot/util/isprop.m, scripts/plot/util/linkprop.m, scripts/plot/util/private/__go_draw_axes__.m, scripts/set/setdiff.m, scripts/set/union.m, scripts/signal/periodogram.m, scripts/sparse/eigs.m, scripts/sparse/ilu.m, scripts/sparse/qmr.m, scripts/sparse/sprand.m, scripts/sparse/sprandn.m, scripts/specfun/beta.m, scripts/specfun/ellipke.m, scripts/specfun/isprime.m, scripts/statistics/base/lscov.m, scripts/testfun/__run_test_suite__.m, scripts/testfun/test.m: Strip trailing whitespace.
author John W. Eaton <jwe@octave.org>
date Tue, 20 Jan 2015 10:29:54 -0500
parents 3ea57d22d9eb
children 4197fc428c7d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
1 ## Copyright (C) 2014 Nathan Podlich
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
2 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
3 ## This file is part of Octave.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
4 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
6 ## it under the terms of the GNU General Public License as published by
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
8 ## (at your option) any later version.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
9 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
10 ## Octave is distributed in the hope that it will be useful,
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
13 ## GNU General Public License for more details.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
14 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
16 ## along with Octave; If not, see <http://www.gnu.org/licenses/>.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
17
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
18 ## -*- texinfo -*-
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
19 ## @deftypefn {Function File} {@var{x} =} qmr (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
20 ## @deftypefnx {Function File} {@var{x} =} qmr (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
21 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} qmr (@var{A}, @var{b}, @dots{})
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
22 ## Solve @code{A x = b} using the Quasi-Minimal Residual iterative
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
23 ## method (without look-ahead).
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
24 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
25 ## @itemize @minus
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
26 ## @item @var{rtol} is the relative tolerance, if not given
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
27 ## or set to [] the default value 1e-6 is used.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
28 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
29 ## @item @var{maxit} the maximum number of outer iterations,
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
30 ## if not given or set to [] the default value
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
31 ## @code{min (20, numel (b))} is used.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
32 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
33 ## @item @var{x0} the initial guess, if not given or set to []
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
34 ## the default value @code{zeros (size (b))} is used.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
35 ## @end itemize
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
36 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
37 ## @var{A} can be passed as a matrix or as a function handle or
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
38 ## inline function @code{f} such that @code{f(x, "notransp") = A*x}
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
39 ## and @code{f(x, "transp") = A'*x}.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
40 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
41 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
42 ## Both @var{M1} and @var{M2} can be passed as a matrix or as
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
43 ## a function handle or inline function @code{g} such that
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
44 ## @code{g(x, "notransp") = M1 \ x} or @code{g(x, "notransp") = M2 \ x} and
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
45 ## @code{g(x, "transp") = M1' \ x} or @code{g(x, "transp") = M2' \ x}.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
46 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
47 ## If called with more than one output parameter
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
48 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
49 ## @itemize @minus
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
50 ## @item @var{flag} indicates the exit status:
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
51 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
52 ## @itemize @minus
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
53 ## @item 0: iteration converged to the within the chosen tolerance
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
54 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
55 ## @item 1: the maximum number of iterations was reached before convergence
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
56 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
57 ## @item 3: the algorithm reached stagnation
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
58 ## @end itemize
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
59 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
60 ## (the value 2 is unused but skipped for compatibility).
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
61 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
62 ## @item @var{relres} is the final value of the relative residual.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
63 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
64 ## @item @var{iter} is the number of iterations performed.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
65 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
66 ## @item @var{resvec} is a vector containing the residual norms at each
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
67 ## iteration.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
68 ## @end itemize
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
69 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
70 ## References:
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
71 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
72 ## @enumerate
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
73 ## @item
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
74 ## R. Freund and N. Nachtigal, @cite{QMR: a quasi-minimal residual
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
75 ## method for non-Hermitian linear systems}, Numerische Mathematik,
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
76 ## 1991, 60, pp. 315-339.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
77 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
78 ## @item
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
79 ## R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra,
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
80 ## V. Eijkhour, R. Pozo, C. Romine, and H. van der Vorst,
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
81 ## @cite{Templates for the solution of linear systems: Building blocks
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
82 ## for iterative methods}, SIAM, 2nd ed., 1994.
19356
3ea57d22d9eb qmr.m: Add missing "@end enumerate" in docstring.
Rik <rik@octave.org>
parents: 19354
diff changeset
83 ## @end enumerate
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
84 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
85 ## @seealso{bicg, bicgstab, cgs, gmres}
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
86 ## @end deftypefn
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
87
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
88 ## Author: Nathan Podlich <nathan.podlich@gmail.com>
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
89
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
90 function [x, flag, relres, iter, resvec] = qmr (A, b, tol, maxit, M1, M2, x0)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
91
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
92 if (nargin >= 2 && isvector (full (b)))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
93
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
94 if (ischar (A))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
95 fun = str2func (A);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
96 Ax = @(x) feval (fun, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
97 Atx = @(x) feval (fun, x, "transp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
98 elseif (isa (A, "function_handle"))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
99 Ax = @(x) feval (A, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
100 Atx = @(x) feval (A, x, "transp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
101 elseif (ismatrix (A))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
102 Ax = @(x) A * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
103 Atx = @(x) A' * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
104 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
105 error (["qmr: first argument is expected to " ...
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
106 "be a function or a square matrix"]);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
107 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
108
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
109 if (nargin < 3 || isempty (tol))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
110 tol = 1e-6;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
111 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
112
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
113 if (nargin < 4 || isempty (maxit))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
114 maxit = min (rows (b), 20);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
115 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
116
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
117 if (nargin < 5 || isempty (M1))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
118 M1m1x = @(x, ignore) x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
119 M1tm1x = M1m1x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
120 elseif (ischar (M1))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
121 fun = str2func (M1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
122 M1m1x = @(x) feval (fun, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
123 M1tm1x = @(x) feval (fun, x, "transp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
124 elseif (isa (M1, "function_handle"))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
125 M1m1x = @(x) feval (M1, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
126 M1tm1x = @(x) feval (M1, x, "transp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
127 elseif (ismatrix (M1))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
128 M1m1x = @(x) M1 \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
129 M1tm1x = @(x) M1' \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
130 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
131 error (["qmr: preconditioner is expected to " ...
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
132 "be a function or matrix"]);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
133 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
134
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
135 if (nargin < 6 || isempty (M2))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
136 M2m1x = @(x, ignore) x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
137 M2tm1x = M2m1x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
138 elseif (ischar (M2))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
139 fun = str2func (M2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
140 M2m1x = @(x) feval (fun, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
141 M2tm1x = @(x) feval (fun, x, "transp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
142 elseif (isa (M2, "function_handle"))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
143 M2m1x = @(x) feval (M2, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
144 M2tm1x = @(x) feval (M2, x, "transp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
145 elseif (ismatrix (M2))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
146 M2m1x = @(x) M2 \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
147 M2tm1x = @(x) M2' \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
148 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
149 error (["qmr: preconditioner is expected to " ...
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
150 "be a function or matrix"]);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
151 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
152
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
153
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
154 if (nargin < 7 || isempty (x0))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
155 x = zeros (size (b));
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
156 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
157 x = x0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
158 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
159
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
160 r = b - Ax (x);
19597
db92e7e28e1f strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents: 19356
diff changeset
161
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
162 bnorm = norm (b);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
163 res0 = norm (r);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
164 if (nargout > 4)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
165 resvec(1) = res0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
166 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
167 vt = r;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
168
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
169 y = M1m1x (vt);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
170
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
171 rho0 = norm (y);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
172 wt = r;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
173
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
174 z = M2tm1x (wt);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
175
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
176 xi1 = norm (z);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
177 gamma0 = 1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
178 eta0 = -1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
179 flag = 1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
180 for iter=1:1:maxit
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
181 ## If rho0 == 0 or xi1 == 0, method fails.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
182 v = vt / rho0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
183 y = y / rho0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
184 w = wt / xi1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
185 z = z / xi1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
186
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
187 delta1 = z' * y; ## If delta1 == 0, method fails.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
188
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
189 yt = M2m1x (y);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
190 zt = M1tm1x (z);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
191
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
192 if (iter == 1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
193 p = yt;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
194 q = zt;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
195 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
196 p = yt - (xi1*delta1/eps0) * p;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
197 q = zt - (rho0*delta1/eps0) * q;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
198 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
199 pt = Ax (p);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
200
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
201 eps0 = q' * pt; ## If eps0 == 0, method fails.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
202 beta1 = eps0 / delta1; ## If beta1 == 0, method fails.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
203 vt = pt - beta1 * v;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
204
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
205 y = M1m1x (vt);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
206 rho1 = norm(y);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
207 wt = Atx (q) - beta1 * w;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
208 z = M2tm1x (wt);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
209
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
210 xi1 = norm(z);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
211 theta1 = rho1 / (gamma0 * abs(beta1));
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
212 gamma1 = 1 / sqrt(1 + theta1^2); ## If gamma1 == 0, method fails.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
213 eta1 = -eta0 * rho0 * gamma1^2 / (beta1 * gamma0^2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
214
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
215 if (iter == 1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
216 d = eta1 * p;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
217 s = eta1 * pt;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
218 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
219 d = eta1 * p + (theta0*gamma1)^2 * d;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
220 s = eta1 * pt + (theta0 * gamma1)^2 * s;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
221 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
222 x += d;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
223 r -= s;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
224
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
225 res1 = norm (r) / bnorm;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
226 if (nargout > 4)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
227 resvec(iter + 1, 1) = norm (r);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
228 end
19597
db92e7e28e1f strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents: 19356
diff changeset
229
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
230 if (res1 < tol)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
231 ## Convergence achieved.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
232 flag = 0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
233 break;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
234 elseif (res0 <= res1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
235 ## Stagnation encountered.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
236 flag = 3;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
237 break;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
238 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
239 theta0 = theta1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
240 eta0 = eta1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
241 gamma0 = gamma1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
242 rho0 = rho1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
243 endfor
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
244
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
245 relres = res1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
246 if (flag == 1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
247 if (nargout < 2)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
248 printf ("qmr stopped at iteration %i ", iter);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
249 printf ("without converging to the desired tolerance %e\n", tol);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
250 printf ("because the maximum number of iterations was reached. ");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
251 printf ("The iterate returned (number %i) has ", maxit);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
252 printf ("relative residual %e\n", res1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
253 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
254 elseif (flag == 3)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
255 if (nargout < 2)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
256 printf ("qmr stopped at iteration %i ", iter);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
257 printf (" without converging to the desired tolerance %e\n", tol);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
258 printf ("because the method stagnated.\n");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
259 printf ("The iterate returned (number %i) ", iter);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
260 printf ("has relative residual %e\n", res1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
261 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
262 elseif (nargout < 2)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
263 printf ("qmr converged at iteration %i ", iter);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
264 printf ("to a solution with relative residual %e\n", res1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
265 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
266 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
267 print usage();
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
268 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
269 endfunction
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
270
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
271
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
272 %!demo
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
273 %! % Solve system of A*x=b
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
274 %! A = [5 -1 3;-1 2 -2;3 -2 3];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
275 %! b = [7;-1;4];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
276 %! [x, flag, relres, iter, resvec] = qmr (A, b)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
277
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
278 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
279 %! n = 100;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
280 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
281 %! b = sum (A, 2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
282 %! tol = 1e-8;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
283 %! maxit = 15;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
284 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
285 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
286 %! [x, flag, relres, iter, resvec] = qmr (A, b, tol, maxit, M1, M2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
287 %! assert (x, ones (size (b)), 1e-7);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
288
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
289 %!function y = afun (x, t, a)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
290 %! switch (t)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
291 %! case "notransp"
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
292 %! y = a * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
293 %! case "transp"
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
294 %! y = a' * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
295 %! endswitch
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
296 %!endfunction
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
297 %!
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
298 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
299 %! n = 100;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
300 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
301 %! b = sum (A, 2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
302 %! tol = 1e-8;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
303 %! maxit = 15;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
304 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
305 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
306 %!
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
307 %! [x, flag, relres, iter, resvec] = qmr (@(x, t) afun (x, t, A),
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
308 %! b, tol, maxit, M1, M2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
309 %! assert (x, ones (size (b)), 1e-7);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
310
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
311 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
312 %! n = 100;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
313 %! tol = 1e-8;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
314 %! a = sprand (n, n, .1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
315 %! A = a' * a + 100 * eye (n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
316 %! b = sum (A, 2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
317 %! [x, flag, relres, iter, resvec] = qmr (A, b, tol, [], diag (diag (A)));
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
318 %! assert (x, ones (size (b)), 1e-7);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
319
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
320 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
321 %! A = [1 + 1i, 1 + 1i; 2 - 1i, 2 + 1i];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
322 %! b = A * [1; 1];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
323 %! [x, flag, relres, iter, resvec] = qmr (A, b);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
324 %! assert (x, [1; 1], 1e-6);