Mercurial > octave
annotate scripts/sparse/bicgstab.m @ 23083:e9a0469dedd9 stable
maint: strip extra trailing newlines from files.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 20 Jan 2017 12:19:08 -0500 |
parents | bac0d6f07a3e |
children | ef4d915df748 3ac9f9ecfae5 |
rev | line source |
---|---|
22323
bac0d6f07a3e
maint: Update copyright notices for 2016.
John W. Eaton <jwe@octave.org>
parents:
21758
diff
changeset
|
1 ## Copyright (C) 2008-2016 Radek Salac |
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
13929
diff
changeset
|
2 ## Copyright (C) 2012 Carlo de Falco |
8416 | 3 ## |
4 ## This file is part of Octave. | |
5 ## | |
6 ## Octave is free software; you can redistribute it and/or modify it | |
7 ## under the terms of the GNU General Public License as published by | |
8 ## the Free Software Foundation; either version 3 of the License, or (at | |
9 ## your option) any later version. | |
10 ## | |
11 ## Octave is distributed in the hope that it will be useful, but | |
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 ## General Public License for more details. | |
15 ## | |
16 ## You should have received a copy of the GNU General Public License | |
17 ## along with Octave; see the file COPYING. If not, see | |
18 ## <http://www.gnu.org/licenses/>. | |
19 | |
20 ## -*- texinfo -*- | |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20711
diff
changeset
|
21 ## @deftypefn {} {@var{x} =} bicgstab (@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
|
22 ## @deftypefnx {} {@var{x} =} bicgstab (@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
|
23 ## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, @dots{}) |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
24 ## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative |
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
25 ## method. |
13070
a0d854f079d2
codesprint: Fix building of docs for new bicg functions
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
13061
diff
changeset
|
26 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
27 ## @itemize @minus |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19854
diff
changeset
|
28 ## @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:
19854
diff
changeset
|
29 ## default value 1e-6 is used. |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
30 ## |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19854
diff
changeset
|
31 ## @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:
19854
diff
changeset
|
32 ## set to [] the default value @code{min (20, numel (b))} is used. |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
33 ## |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19854
diff
changeset
|
34 ## @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:
19854
diff
changeset
|
35 ## value @code{zeros (size (b))} is used. |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
36 ## @end itemize |
13070
a0d854f079d2
codesprint: Fix building of docs for new bicg functions
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
13061
diff
changeset
|
37 ## |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19854
diff
changeset
|
38 ## @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:
19854
diff
changeset
|
39 ## function @code{f} such that @code{f(x) = A*x}. |
8416 | 40 ## |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19854
diff
changeset
|
41 ## 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:
19854
diff
changeset
|
42 ## 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:
19854
diff
changeset
|
43 ## function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}. |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
44 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
45 ## If called with more than one output parameter |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
46 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
47 ## @itemize @minus |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
48 ## @item @var{flag} indicates the exit status: |
14853
72b8b39e12be
doc: Periodic grammarcheck of documentation.
Rik <octave@nomad.inbox5.com>
parents:
14609
diff
changeset
|
49 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
50 ## @itemize @minus |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
51 ## @item 0: iteration converged to the within the chosen tolerance |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
52 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
53 ## @item 1: the maximum number of iterations was reached before convergence |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
54 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
55 ## @item 3: the algorithm reached stagnation |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
56 ## @end itemize |
16816
12005245b645
doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents:
14868
diff
changeset
|
57 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
58 ## (the value 2 is unused but skipped for compatibility). |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
59 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
60 ## @item @var{relres} is the final value of the relative residual. |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
61 ## |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
62 ## @item @var{iter} is the number of iterations performed. |
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13837
diff
changeset
|
63 ## |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19700
diff
changeset
|
64 ## @item @var{resvec} is a vector containing the relative residual at each |
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19700
diff
changeset
|
65 ## iteration. |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
66 ## @end itemize |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
67 ## |
19854 | 68 ## @seealso{bicg, cgs, gmres, pcg, qmr} |
8416 | 69 ## |
70 ## @end deftypefn | |
71 | |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
72 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
73 M1, M2, x0) |
8416 | 74 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
75 if (nargin < 2 || nargin > 7 || ! isvector (full (b))) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
76 print_usage (); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
77 endif |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
78 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
79 if (ischar (A)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
80 A = str2func (A); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
81 elseif (isnumeric(A) && issquare (A)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
82 Ax = @(x) A * x; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
83 elseif (isa (A, "function_handle")) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
84 Ax = @(x) feval (A, x); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
85 else |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
86 error ("bicgstab: A must be a square matrix or function"); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
87 endif |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
88 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
89 if (nargin < 3 || isempty (tol)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
90 tol = 1e-6; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
91 endif |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
92 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
93 if (nargin < 4 || isempty (maxit)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
94 maxit = min (rows (b), 20); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
95 endif |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
96 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
97 if (nargin < 5 || isempty (M1)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
98 M1m1x = @(x) x; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
99 elseif (ischar (M1)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
100 M1m1x = str2func (M1); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
101 elseif (isnumeric(M1) && ismatrix (M1)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
102 M1m1x = @(x) M1 \ x; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
103 elseif (isa (M1, "function_handle")) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
104 M1m1x = @(x) feval (M1, x); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
105 else |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
106 error ("bicgstab: preconditioner M1 must be a function or matrix"); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
107 endif |
8416 | 108 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
109 if (nargin < 6 || isempty (M2)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
110 M2m1x = @(x) x; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
111 elseif (ischar (M2)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
112 M2m1x = str2func (M2); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
113 elseif (isnumeric(M2) && ismatrix (M2)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
114 M2m1x = @(x) M2 \ x; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
115 elseif (isa (M2, "function_handle")) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
116 M2m1x = @(x) feval (M2, x); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
117 else |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
118 error ("bicgstab: preconditioner M2 must be a function or matrix"); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
119 endif |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
120 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
121 precon = @(x) M2m1x (M1m1x (x)); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
122 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
123 if (nargin < 7 || isempty (x0)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
124 x0 = zeros (size (b)); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
125 endif |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
126 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
127 ## specifies initial estimate x0 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
128 if (nargin < 7) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
129 x = zeros (rows (b), 1); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
130 else |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
131 x = x0; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
132 endif |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
133 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
134 norm_b = norm (b); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
135 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
136 res = b - Ax (x); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
137 rr = res; |
8416 | 138 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
139 ## Vector of the residual norms for each iteration. |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
140 resvec = norm (res) / norm_b; |
8416 | 141 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
142 ## Default behavior we don't reach tolerance tol within maxit iterations. |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
143 flag = 1; |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
144 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
145 for iter = 1:maxit |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
146 rho_1 = rr' * res; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
147 |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
148 if (iter == 1) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
149 p = res; |
8416 | 150 else |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
151 beta = (rho_1 / rho_2) * (alpha / omega); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
152 p = res + beta * (p - omega * v); |
8416 | 153 endif |
154 | |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
155 phat = precon (p); |
8416 | 156 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
157 v = Ax (phat); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
158 alpha = rho_1 / (rr' * v); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
159 s = res - alpha * v; |
8416 | 160 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
161 shat = precon (s); |
8416 | 162 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
163 t = Ax (shat); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
164 omega = (s' * t) / (t' * t); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
165 x += alpha * phat + omega * shat; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
166 res = s - omega * t; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
167 rho_2 = rho_1; |
8416 | 168 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
169 relres = norm (res) / norm_b; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
170 resvec = [resvec; relres]; |
8416 | 171 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
172 if (relres <= tol) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
173 ## We reach tolerance tol within maxit iterations. |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
174 flag = 0; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
175 break; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
176 elseif (resvec(end) == resvec(end - 1)) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
177 ## The method stagnates. |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
178 flag = 3; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
179 break; |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
180 endif |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
181 endfor |
8416 | 182 |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
183 if (nargout < 2) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
184 if (flag == 0) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
185 printf ("bicgstab converged at iteration %i ", iter); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
186 printf ("to a solution with relative residual %e\n", relres); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
187 elseif (flag == 3) |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
188 printf ("bicgstab stopped at iteration %i ", iter); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
189 printf ("without converging to the desired tolerance %e\n", tol); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
190 printf ("because the method stagnated.\n"); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
191 printf ("The iterate returned (number %i) ", iter); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
192 printf ("has relative residual %e\n", relres); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
193 else |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
194 printf ("bicgstab stopped at iteration %i ", iter); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
195 printf ("without converging to the desired toleranc %e\n", tol); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
196 printf ("because the maximum number of iterations was reached.\n"); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
197 printf ("The iterate returned (number %i) ", iter); |
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
21634
diff
changeset
|
198 printf ("has relative residual %e\n", relres); |
9278 | 199 endif |
200 endif | |
201 | |
8416 | 202 endfunction |
203 | |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
204 |
8416 | 205 %!demo |
206 %! % Solve system of A*x=b | |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
207 %! A = [5 -1 3;-1 2 -2;3 -2 3]; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
208 %! b = [7;-1;4]; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
209 %! [x, flag, relres, iter, resvec] = bicgstab (A, b) |
8416 | 210 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
211 %!shared A, b, n, M1, M2 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
212 %! |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
213 %!test |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
214 %! n = 100; |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
215 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
216 %! b = sum (A, 2); |
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
217 %! tol = 1e-8; |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
218 %! maxit = 15; |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
219 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); |
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
220 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
221 %! [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, M1, M2); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
222 %! assert (x, ones (size (b)), 1e-7); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
223 %! |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
224 %!test |
13837
2c80bbd87f5d
don't define functions in test and demo blocks
John W. Eaton <jwe@octave.org>
parents:
13141
diff
changeset
|
225 %!function y = afun (x, a) |
2c80bbd87f5d
don't define functions in test and demo blocks
John W. Eaton <jwe@octave.org>
parents:
13141
diff
changeset
|
226 %! y = a * x; |
2c80bbd87f5d
don't define functions in test and demo blocks
John W. Eaton <jwe@octave.org>
parents:
13141
diff
changeset
|
227 %!endfunction |
2c80bbd87f5d
don't define functions in test and demo blocks
John W. Eaton <jwe@octave.org>
parents:
13141
diff
changeset
|
228 %! |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
229 %! tol = 1e-8; |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
230 %! maxit = 15; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
231 %! |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
232 %! [x, flag, relres, iter, resvec] = bicgstab (@(x) afun (x, A), b, |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
233 %! tol, maxit, M1, M2); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
234 %! assert (x, ones (size (b)), 1e-7); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
235 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
236 %!test |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
237 %! n = 100; |
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
238 %! tol = 1e-8; |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
239 %! a = sprand (n, n, .1); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
240 %! A = a'*a + 100 * eye (n); |
13141
e81ddf9cacd5
maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
13091
diff
changeset
|
241 %! b = sum (A, 2); |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
242 %! [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, [], diag (diag (A))); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
243 %! assert (x, ones (size (b)), 1e-7); |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14237
diff
changeset
|
244 |
17422
7ad3b9ca66f5
fix behaviour of bicgstab for complex matrices (bug #40017)
Marco Caliari <marco.caliari@univr.it>
parents:
16816
diff
changeset
|
245 %!test |
7ad3b9ca66f5
fix behaviour of bicgstab for complex matrices (bug #40017)
Marco Caliari <marco.caliari@univr.it>
parents:
16816
diff
changeset
|
246 %! A = [1 + 1i, 1 + 1i; 2 - 1i, 2 + 1i]; |
7ad3b9ca66f5
fix behaviour of bicgstab for complex matrices (bug #40017)
Marco Caliari <marco.caliari@univr.it>
parents:
16816
diff
changeset
|
247 %! b = A * [1; 1]; |
7ad3b9ca66f5
fix behaviour of bicgstab for complex matrices (bug #40017)
Marco Caliari <marco.caliari@univr.it>
parents:
16816
diff
changeset
|
248 %! [x, flag, relres, iter, resvec] = bicgstab (A, b); |
7ad3b9ca66f5
fix behaviour of bicgstab for complex matrices (bug #40017)
Marco Caliari <marco.caliari@univr.it>
parents:
16816
diff
changeset
|
249 %! assert (x, [1; 1], 1e-6); |