Mercurial > octave

annotate scripts/sparse/qmr.m @ 31205:b0e90ca8e679 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
quad2d: fix unintended complex conjugate return (bug #62972) quad2d: use .' instead of ' to avoid complex conjugate in q and qerr.
author Arun Giridhar <arungiridhar@gmail.com>
date Sun, 28 Aug 2022 12:21:17 -0400
parents 796f54d4ddbf
children e1788b1a315f
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ########################################################################
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 ##
30564
796f54d4ddbf update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 29359
diff changeset
3 ## Copyright (C) 2014-2022 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
4 ##
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
7 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
8 ## This file is part of Octave.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
9 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
11 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
13 ## (at your option) any later version.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
14 ##
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
15 ## Octave is distributed in the hope that it will be useful, but
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
18 ## GNU General Public License for more details.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
19 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
20 ## You should have received a copy of the GNU General Public License
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
21 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
22 ## <https://www.gnu.org/licenses/>.
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 ##
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ########################################################################
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
25
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
26 ## -*- texinfo -*-
20852
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20711
diff changeset
27 ## @deftypefn {} {@var{x} =} qmr (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20711
diff changeset
28 ## @deftypefnx {} {@var{x} =} qmr (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
516bb87ea72e 2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents: 20711
diff changeset
29 ## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} qmr (@var{A}, @var{b}, @dots{})
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
30 ## Solve @code{A x = b} using the Quasi-Minimal Residual iterative method
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
31 ## (without look-ahead).
19354
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 ## @itemize @minus
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
34 ## @item @var{rtol} is the relative tolerance, if not given or set to [] the
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
35 ## default value 1e-6 is used.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
36 ##
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
37 ## @item @var{maxit} the maximum number of outer iterations, if not given or
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
38 ## set to [] the default value @code{min (20, numel (b))} is used.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
39 ##
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
40 ## @item @var{x0} the initial guess, if not given or set to [] the default
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
41 ## value @code{zeros (size (b))} is used.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
42 ## @end itemize
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
43 ##
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
44 ## @var{A} can be passed as a matrix or as a function handle or inline
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
45 ## function @code{f} such that @code{f(x, "notransp") = A*x} and
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
46 ## @code{f(x, "transp") = A'*x}.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
47 ##
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
48 ## The preconditioner @var{P} is given as @code{P = M1 * M2}. Both @var{M1}
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
49 ## and @var{M2} can be passed as a matrix or as a function handle or inline
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
50 ## function @code{g} such that @code{g(x, "notransp") = M1 \ x} or
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
51 ## @code{g(x, "notransp") = M2 \ x} and @code{g(x, "transp") = M1' \ x} or
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 20101
diff changeset
52 ## @code{g(x, "transp") = M2' \ x}.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
53 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
54 ## If called with more than one output parameter
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
55 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
56 ## @itemize @minus
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
57 ## @item @var{flag} indicates the exit status:
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
58 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
59 ## @itemize @minus
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
60 ## @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
61 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
62 ## @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
63 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
64 ## @item 3: the algorithm reached stagnation
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
65 ## @end itemize
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
66 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
67 ## (the value 2 is unused but skipped for compatibility).
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
68 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
69 ## @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
70 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
71 ## @item @var{iter} is the number of iterations performed.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
72 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
73 ## @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
74 ## iteration.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
75 ## @end itemize
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
76 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
77 ## References:
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
78 ##
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
79 ## @enumerate
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
80 ## @item
20098
21dbaf135667 doc: Periodic spellcheck of documentation.
Rik <rik@octave.org>
parents: 19854
diff changeset
81 ## @nospell{R. Freund and N. Nachtigal}, @cite{QMR: a quasi-minimal residual
21dbaf135667 doc: Periodic spellcheck of documentation.
Rik <rik@octave.org>
parents: 19854
diff changeset
82 ## method for non-Hermitian linear systems}, @nospell{Numerische Mathematik},
27800
5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations.
Rik <rik@octave.org>
parents: 27585
diff changeset
83 ## 1991, 60, pp.@: 315--339.
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 ## @item
20101
e51473fdb622 doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents: 20098
diff changeset
86 ## @nospell{ R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra},
e51473fdb622 doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents: 20098
diff changeset
87 ## @nospell{ V. Eijkhour, R. Pozo, C. Romine, and H. van der Vorst},
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
88 ## @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
89 ## 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
90 ## @end enumerate
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
91 ##
19854
0a3ca546d7fc doc: Add qmr to documentation.
Rik <rik@octave.org>
parents: 19833
diff changeset
92 ## @seealso{bicg, bicgstab, cgs, gmres, pcg}
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
93 ## @end deftypefn
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
94
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
95 function [x, flag, relres, iter, resvec] = qmr (A, b, rtol, maxit, M1, M2, x0)
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
96
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
97 if (nargin >= 2 && isvector (full (b)))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
98
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
99 if (ischar (A))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
100 fun = str2func (A);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
101 Ax = @(x) feval (fun, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
102 Atx = @(x) feval (fun, x, "transp");
25803
23483673ba43 Use is_function_handle instead of isa (x, "function_handle").
Rik <rik@octave.org>
parents: 25054
diff changeset
103 elseif (is_function_handle (A))
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
104 Ax = @(x) feval (A, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
105 Atx = @(x) feval (A, x, "transp");
21171
2935d56203a4 Fix regressions caused by ismatrix definition change (partial fix bug #47036).
Rik <rik@octave.org>
parents: 20164
diff changeset
106 elseif (isnumeric (A) && issquare (A))
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
107 Ax = @(x) A * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
108 Atx = @(x) A' * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
109 else
20711
7b608fadc663 Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents: 20231
diff changeset
110 error ("qmr: A must be a square matrix or function");
19354
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
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
113 if (nargin < 3 || isempty (rtol))
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
114 rtol = 1e-6;
19354
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 < 4 || isempty (maxit))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
118 maxit = min (rows (b), 20);
19725
5f2c0ca0ef51 Ensure that numbers passed to integer *printf format codes are integers (bug #44245).
Rik <rik@octave.org>
parents: 19700
diff changeset
119 else
5f2c0ca0ef51 Ensure that numbers passed to integer *printf format codes are integers (bug #44245).
Rik <rik@octave.org>
parents: 19700
diff changeset
120 maxit = fix (maxit);
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
121 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
122
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
123 if (nargin < 5 || isempty (M1))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
124 M1m1x = @(x, ignore) x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
125 M1tm1x = M1m1x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
126 elseif (ischar (M1))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
127 fun = str2func (M1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
128 M1m1x = @(x) feval (fun, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
129 M1tm1x = @(x) feval (fun, x, "transp");
25803
23483673ba43 Use is_function_handle instead of isa (x, "function_handle").
Rik <rik@octave.org>
parents: 25054
diff changeset
130 elseif (is_function_handle (M1))
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
131 M1m1x = @(x) feval (M1, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
132 M1tm1x = @(x) feval (M1, x, "transp");
19700
00e31f316a3a Fix Matlab incompatibility of "ismatrix" (bug #42422).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents: 19697
diff changeset
133 elseif (isnumeric (M1) && ismatrix (M1))
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
134 M1m1x = @(x) M1 \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
135 M1tm1x = @(x) M1' \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
136 else
21171
2935d56203a4 Fix regressions caused by ismatrix definition change (partial fix bug #47036).
Rik <rik@octave.org>
parents: 20164
diff changeset
137 error ("qmr: preconditioner M1 must be a function or matrix");
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
138 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
139
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
140 if (nargin < 6 || isempty (M2))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
141 M2m1x = @(x, ignore) x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
142 M2tm1x = M2m1x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
143 elseif (ischar (M2))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
144 fun = str2func (M2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
145 M2m1x = @(x) feval (fun, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
146 M2tm1x = @(x) feval (fun, x, "transp");
25803
23483673ba43 Use is_function_handle instead of isa (x, "function_handle").
Rik <rik@octave.org>
parents: 25054
diff changeset
147 elseif (is_function_handle (M2))
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
148 M2m1x = @(x) feval (M2, x, "notransp");
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
149 M2tm1x = @(x) feval (M2, x, "transp");
19700
00e31f316a3a Fix Matlab incompatibility of "ismatrix" (bug #42422).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents: 19697
diff changeset
150 elseif (isnumeric (M2) && ismatrix (M2))
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
151 M2m1x = @(x) M2 \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
152 M2tm1x = @(x) M2' \ x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
153 else
21171
2935d56203a4 Fix regressions caused by ismatrix definition change (partial fix bug #47036).
Rik <rik@octave.org>
parents: 20164
diff changeset
154 error ("qmr: preconditioner M2 must be a function or matrix");
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
155 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
156
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
157 if (nargin < 7 || isempty (x0))
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
158 x = zeros (size (b));
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
159 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
160 x = x0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
161 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
162
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
163 r = b - Ax (x);
19597
db92e7e28e1f strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents: 19356
diff changeset
164
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
165 bnorm = norm (b);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
166 res0 = norm (r);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
167 if (nargout > 4)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
168 resvec(1) = res0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
169 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
170 vt = r;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
171
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
172 y = M1m1x (vt);
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 rho0 = norm (y);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
175 wt = r;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
176
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
177 z = M2tm1x (wt);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
178
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
179 xi1 = norm (z);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
180 gamma0 = 1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
181 eta0 = -1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
182 flag = 1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
183 for iter=1:1:maxit
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
184 ## If rho0 == 0 or xi1 == 0, method fails.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
185 v = vt / rho0;
20231
83792dd9bcc1 Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents: 20164
diff changeset
186 y /= rho0;
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
187 w = wt / xi1;
20231
83792dd9bcc1 Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents: 20164
diff changeset
188 z /= xi1;
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
189
19833
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19725
diff changeset
190 delta1 = z' * y; # If delta1 == 0, method fails.
19354
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 yt = M2m1x (y);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
193 zt = M1tm1x (z);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
194
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
195 if (iter == 1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
196 p = yt;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
197 q = zt;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
198 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
199 p = yt - (xi1*delta1/eps0) * p;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
200 q = zt - (rho0*delta1/eps0) * q;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
201 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
202 pt = Ax (p);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
203
19833
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19725
diff changeset
204 eps0 = q' * pt; # If eps0 == 0, method fails.
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19725
diff changeset
205 beta1 = eps0 / delta1; # If beta1 == 0, method fails.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
206 vt = pt - beta1 * v;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
207
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
208 y = M1m1x (vt);
28912
0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents: 27978
diff changeset
209 rho1 = norm (y);
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
210 wt = Atx (q) - beta1 * w;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
211 z = M2tm1x (wt);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
212
28912
0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents: 27978
diff changeset
213 xi1 = norm (z);
0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents: 27978
diff changeset
214 theta1 = rho1 / (gamma0 * abs (beta1));
0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents: 27978
diff changeset
215 gamma1 = 1 / sqrt (1 + theta1^2); # If gamma1 == 0, method fails.
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
216 eta1 = -eta0 * rho0 * gamma1^2 / (beta1 * gamma0^2);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
217
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
218 if (iter == 1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
219 d = eta1 * p;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
220 s = eta1 * pt;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
221 else
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
222 d = eta1 * p + (theta0*gamma1)^2 * d;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
223 s = eta1 * pt + (theta0 * gamma1)^2 * s;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
224 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
225 x += d;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
226 r -= s;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
227
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
228 res1 = norm (r) / bnorm;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
229 if (nargout > 4)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
230 resvec(iter + 1, 1) = norm (r);
21031
66a08c3cafe3 maint: Follow Octave coding conventions in m-files.
Rik <rik@octave.org>
parents: 20852
diff changeset
231 endif
19597
db92e7e28e1f strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents: 19356
diff changeset
232
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
233 if (res1 < rtol)
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
234 ## Convergence achieved.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
235 flag = 0;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
236 break;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
237 elseif (res0 <= res1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
238 ## Stagnation encountered.
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
239 flag = 3;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
240 break;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
241 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
242 theta0 = theta1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
243 eta0 = eta1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
244 gamma0 = gamma1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
245 rho0 = rho1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
246 endfor
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
247
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
248 relres = res1;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
249 if (flag == 1)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
250 if (nargout < 2)
21634
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
251 printf ("qmr stopped at iteration %i ", iter);
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
252 printf ("without converging to the desired tolerance %e\n", rtol);
21634
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
253 printf ("because the maximum number of iterations was reached. ");
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
254 printf ("The iterate returned (number %i) has ", maxit);
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
255 printf ("relative residual %e\n", res1);
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
256 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
257 elseif (flag == 3)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
258 if (nargout < 2)
21634
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
259 printf ("qmr stopped at iteration %i ", iter);
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
260 printf (" without converging to the desired tolerance %e\n", rtol);
21634
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
261 printf ("because the method stagnated.\n");
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
262 printf ("The iterate returned (number %i) ", iter);
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
263 printf ("has relative residual %e\n", res1);
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
264 endif
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
265 elseif (nargout < 2)
21634
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
266 printf ("qmr converged at iteration %i ", iter);
96518f623c91 Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents: 21633
diff changeset
267 printf ("to a solution with relative residual %e\n", res1);
19354
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 else
28912
0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir.
Rik <rik@octave.org>
parents: 27978
diff changeset
270 print_usage ();
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
271 endif
21758
ffad2baa90f7 maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents: 21634
diff changeset
272
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
273 endfunction
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
274
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
275
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
276 %!demo
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
277 %! % Solve system of A*x=b
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
278 %! 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
279 %! b = [7;-1;4];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
280 %! [x, flag, relres, iter, resvec] = qmr (A, b)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
281
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
282 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
283 %! n = 100;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
284 %! 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
285 %! b = sum (A, 2);
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
286 %! rtol = 1e-8;
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
287 %! maxit = 15;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
288 %! 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
289 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
290 %! [x, flag, relres, iter, resvec] = qmr (A, b, rtol, maxit, M1, M2);
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
291 %! assert (x, ones (size (b)), 1e-7);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
292
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
293 %!function y = afun (x, t, a)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
294 %! switch (t)
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
295 %! case "notransp"
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
296 %! y = a * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
297 %! case "transp"
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
298 %! y = a' * x;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
299 %! endswitch
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
300 %!endfunction
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
301 %!
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
302 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
303 %! n = 100;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
304 %! 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
305 %! b = sum (A, 2);
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
306 %! rtol = 1e-8;
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
307 %! maxit = 15;
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
308 %! 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
309 %! 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
310 %!
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
311 %! [x, flag, relres, iter, resvec] = qmr (@(x, t) afun (x, t, A),
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
312 %! b, rtol, maxit, M1, M2);
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
313 %! assert (x, ones (size (b)), 1e-7);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
314
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
315 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
316 %! n = 100;
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
317 %! rtol = 1e-8;
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
318 %! a = sprand (n, n, .1);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
319 %! A = a' * a + 100 * eye (n);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
320 %! b = sum (A, 2);
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
321 %! [x, flag, relres, iter, resvec] = qmr (A, b, rtol, [], diag (diag (A)));
19354
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
322 %! assert (x, ones (size (b)), 1e-7);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
323
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
324 %!test
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
325 %! A = [1 + 1i, 1 + 1i; 2 - 1i, 2 + 1i];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
326 %! b = A * [1; 1];
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
327 %! [x, flag, relres, iter, resvec] = qmr (A, b);
cc6000ccd627 New function qmr (patch #8569).
Nathan Podlich <nathan.podlich@gmail.com>
parents:
diff changeset
328 %! assert (x, [1; 1], 1e-6);