Mercurial > octave-antonio
annotate scripts/sparse/bicgstab.m @ 13091:e5aaba072d2b
maint: style fixes in iterative linear solvers
* bicg.m, bicgstab.m, cgs.m, gmres.m: Style fixes.
author | Carlo de Falco <kingcrimson@tiscali.it> |
---|---|
date | Sun, 04 Sep 2011 01:21:43 +0200 |
parents | a0d854f079d2 |
children | e81ddf9cacd5 |
rev | line source |
---|---|
11523 | 1 ## Copyright (C) 2008-2011 Radek Salac |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
2 ## Copyright (C) 2011 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 -*- | |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
21 ## |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
22 ## @deftypefn {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
23 ## @deftypefnx {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P}) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
24 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, ...) |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
25 ## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative 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 |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
28 ## @item @var{rtol} is the relative tolerance, if not given or set to |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
29 ## [] the default value 1e-6 is used. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
30 ## @item @var{maxit} the maximum number of outer iterations, if not |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
31 ## given or set to [] the default value @code{min (20, numel (b))} is |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
32 ## used. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
33 ## @item @var{x0} the initial guess, if not given or set to [] the |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
34 ## default value @code{zeros (size (b))} is used. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
35 ## @end itemize |
13070
a0d854f079d2
codesprint: Fix building of docs for new bicg functions
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
13061
diff
changeset
|
36 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
37 ## @var{A} can be passed as a matrix or as a function handle or |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
38 ## inline function @code{f} such that @code{f(x) = A*x}. |
8416 | 39 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
40 ## The preconditioner @var{P} is given as @code{P = M1 * M2}. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
41 ## Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
42 ## inline 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
|
43 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
44 ## 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
|
45 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
46 ## @itemize @minus |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
47 ## @item @var{flag} indicates the exit status: |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
48 ## @itemize @minus |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
49 ## @item 0: iteration converged to the within the chosen tolerance |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
50 ## @item 1: the maximum number of iterations was reached before convergence |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
51 ## @item 3: the algorithm reached stagnation |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
52 ## @end itemize |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
53 ## (the value 2 is unused but skipped for compatibility). |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
54 ## @item @var{relres} is the final value of the relative residual. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
55 ## @item @var{iter} is the number of iterations performed. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
56 ## @item @var{resvec} is a vector containing the relative residual at each iteration. |
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
57 ## @end itemize |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
58 ## |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
59 ## @seealso{bicg,cgs,gmres,pcg} |
8416 | 60 ## |
61 ## @end deftypefn | |
62 | |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
63 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
64 M1, M2, x0) |
8416 | 65 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
66 if (nargin >= 2 && nargin <= 7 && isvector (full (b))) |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
67 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
68 if (ischar (A)) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
69 A = str2func (A); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
70 elseif (ismatrix (A)) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
71 Ax = @(x) A * x; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
72 elseif (isa (A, "function_handle")) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
73 Ax = @(x) feval (A, x); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
74 else |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
75 error (["bicgstab: first argument is expected " ... |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
76 "to be a function or a square matrix"]); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
77 endif |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
78 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
79 if (nargin < 3 || isempty (tol)) |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
80 tol = 1e-6; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
81 endif |
8416 | 82 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
83 if (nargin < 4 || isempty (maxit)) |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
84 maxit = min (rows (b), 20); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
85 endif |
8416 | 86 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
87 if (nargin < 5 || isempty (M1)) |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
88 M1m1x = @(x) x; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
89 elseif (ischar (M1)) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
90 M1m1x = str2func (M1); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
91 elseif (ismatrix (M1)) |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
92 M1m1x = @(x) M1 \ x; |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
93 elseif (isa (M1, "function_handle")) |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
94 M1m1x = @(x) feval (M1, x); |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
95 else |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
96 error (["bicgstab: preconditioner is " ... |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
97 "expected to be a function or matrix"]); |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
98 endif |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
99 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
100 if (nargin < 6 || isempty (M2)) |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
101 M2m1x = @(x) x; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
102 elseif (ischar (M2)) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
103 M2m1x = str2func (M2); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
104 elseif (ismatrix (M2)) |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
105 M2m1x = @(x) M2 \ x; |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
106 elseif (isa (M2, "function_handle")) |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
107 M2m1x = @(x) feval (M2, x); |
8416 | 108 else |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
109 error (["bicgstab: preconditioner is "... |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
110 "expected to be a function or matrix"]); |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
111 endif |
8416 | 112 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
113 precon = @(x) M2m1x (M1m1x (x)); |
8416 | 114 |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
115 if (nargin < 7 || isempty (x0)) |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
116 x0 = zeros (size (b)); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
117 endif |
13091
e5aaba072d2b
maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents:
13070
diff
changeset
|
118 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
119 ## specifies initial estimate x0 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
120 if (nargin < 7) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
121 x = zeros (rows (b), 1); |
8416 | 122 else |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
123 x = x0; |
8416 | 124 endif |
125 | |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
126 norm_b = norm (b); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
127 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
128 res = b - Ax (x); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
129 rr = res; |
8416 | 130 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
131 ## Vector of the residual norms for each iteration. |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
132 resvec = norm(res) / norm_b; |
8416 | 133 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
134 ## Default behaviour we don't reach tolerance tol within maxit iterations. |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
135 flag = 1; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
136 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
137 for iter = 1:maxit |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
138 rho_1 = res' * rr; |
8416 | 139 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
140 if (iter == 1) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
141 p = res; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
142 else |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
143 beta = (rho_1 / rho_2) * (alpha / omega); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
144 p = res + beta * (p - omega * v); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
145 endif |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
146 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
147 phat = precon (p); |
8416 | 148 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
149 v = Ax (phat); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
150 alpha = rho_1 / (rr' * v); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
151 s = res - alpha * v; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
152 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
153 shat = precon (s); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
154 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
155 t = Ax (shat); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
156 omega = (t' * s) / (t' * t); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
157 x = x + alpha * phat + omega * shat; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
158 res = s - omega * t; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
159 rho_2 = rho_1; |
8416 | 160 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
161 relres = norm (res) / norm_b; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
162 resvec = [resvec; relres]; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
163 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
164 if (relres <= tol) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
165 ## We reach tolerance tol within maxit iterations. |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
166 flag = 0; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
167 break; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
168 elseif (resvec(end) == resvec(end - 1)) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
169 ## The method stagnates. |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
170 flag = 3; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
171 break; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
172 endif |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
173 endfor |
8416 | 174 |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
175 if (nargout < 2) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
176 if (flag == 0) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
177 printf ("bicgstab converged at iteration %i ", iter); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
178 printf ("to a solution with relative residual %e\n", relres); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
179 elseif (flag == 3) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
180 printf ("bicgstab stopped at iteration %i ", iter); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
181 printf ("without converging to the desired tolerance %e\n", tol); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
182 printf ("because the method stagnated.\n"); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
183 printf ("The iterate returned (number %i) ", iter); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
184 printf ("has relative residual %e\n", relres); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
185 else |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
186 printf ("bicgstab stopped at iteration %i ", iter); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
187 printf ("without converging to the desired toleranc %e\n", tol); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
188 printf ("because the maximum number of iterations was reached.\n"); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
189 printf ("The iterate returned (number %i) ", iter); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
190 printf ("has relative residual %e\n", relres); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
191 endif |
9278 | 192 endif |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
193 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
194 else |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
195 print_usage (); |
9278 | 196 endif |
197 | |
8416 | 198 endfunction |
199 | |
200 %!demo | |
201 %! % Solve system of A*x=b | |
202 %! A = [5 -1 3;-1 2 -2;3 -2 3] | |
203 %! b = [7;-1;4] | |
204 %! [x, flag, relres, iter, resvec] = bicgstab(A, b) | |
205 | |
13061
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
206 %!shared A, b, n, M1, M2 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
207 %! |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
208 %!test |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
209 %! n = 100; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
210 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
211 %! b = sum (A, 2); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
212 %! tol = 1e-8; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
213 %! maxit = 15; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
214 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
215 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
216 %! [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
|
217 %! assert (x, ones (size (b)), 1e-7); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
218 %! |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
219 %!test |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
220 %! tol = 1e-8; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
221 %! maxit = 15; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
222 %! |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
223 %! function y = afun (x, a) |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
224 %! y = a * x; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
225 %! endfunction |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
226 %! |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
227 %! [x, flag, relres, iter, resvec] = bicgstab (@(x) afun (x, A), b, |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
228 %! tol, maxit, M1, M2); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
229 %! assert (x, ones (size (b)), 1e-7); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
230 |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
231 %!test |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
232 %! n = 100; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
233 %! tol = 1e-8; |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
234 %! a = sprand (n, n, .1); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
235 %! A = a'*a + 100 * eye (n); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
236 %! b = sum (A, 2); |
addfc0ae69c0
Make bicgstab interface more compatible
Carlo de Falco <kingcrimson@tiscali.it>
parents:
11587
diff
changeset
|
237 %! [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
|
238 %! assert (x, ones (size (b)), 1e-7); |