Mercurial > hg > octave-lyh

comparison scripts/sparse/bicg.m @ 13073:51bc892d5cf8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move bicg.m to scripts/sparse for uniformity * bicg.m : Move to scripts/sparse where the other iterative solvers are.
author Carlo de Falco <kingcrimson@tiscali.it>
date Sat, 03 Sep 2011 21:51:26 +0200
parents scripts/linear-algebra/bicg.m@a0d854f079d2
children e5aaba072d2b
comparison
equal deleted inserted replaced
13072:39814ef3278b 13073:51bc892d5cf8
1 ## Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com>
2 ## Copyright (C) 2011 Carlo de Falco
3 ##
4 ## This program is free software; you can redistribute it and/or modify
5 ## it under the terms of the GNU General Public License as published by
6 ## the Free Software Foundation; either version 2 of the License, or
7 ## (at your option) any later version.
8 ##
9 ## This program is distributed in the hope that it will be useful,
10 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ## GNU General Public License for more details.
13 ##
14 ## You should have received a copy of the GNU General Public License
15 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
16
17 ## -*- texinfo -*-
18 ##
19 ## @deftypefn {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
20 ## @deftypefnx {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
21 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, ...)
22 ##
23 ## Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
24 ##
25 ## @itemize @minus
26 ## @item @var{rtol} is the relative tolerance, if not given or set to [] the default value 1e-6 is used.
27 ## @item @var{maxit} the maximum number of outer iterations, if not given or set to [] the default value @code{min (20, numel (b))} is used.
28 ## @item @var{x0} the initial guess, if not given or set to [] the default value @code{zeros (size (b))} is used.
29 ## @end itemize
30 ##
31 ## @var{A} can be passed as a matrix or as a function handle or
32 ## inline function @code{f} such that @code{f(x, "notransp") = A*x} and @code{f(x, "transp") = A'*x}.
33 ##
34 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
35 ## Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or
36 ## inline function @code{g} such that @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and
37 ## @code{g(x, 'transp') = M1' \ x} or @code{g(x, 'transp') = M2' \ x}.
38 ##
39 ## If colled with more than one output parameter
40 ##
41 ## @itemize @minus
42 ## @item @var{flag} indicates the exit status:
43 ## @itemize @minus
44 ## @item 0: iteration converged to the within the chosen tolerance
45 ## @item 1: the maximum number of iterations was reached before convergence
46 ## @item 3: the algorithm reached stagnation
47 ## @end itemize
48 ## (the value 2 is unused but skipped for compatibility).
49 ## @item @var{relres} is the final value of the relative residual.
50 ## @item @var{iter} is the number of iterations performed.
51 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
52 ## @end itemize
53 ##
54 ## @seealso{pcg,cgs,bicgstab,gmres}
55 ##
56 ## @end deftypefn
57
58
59 function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0)
60
61 if ((nargin >= 2) && isvector (full (b)))
62
63 if (ischar (A))
64 fun = str2func (A);
65 Ax = @(x) feval (fun, x, "notransp");
66 Atx = @(x) feval (fun, x, "transp");
67 elseif (ismatrix (A))
68 Ax = @(x) A * x;
69 Atx = @(x) A' * x;
70 elseif (isa (A, "function_handle"))
71 Ax = @(x) feval (A, x, "notransp");
72 Atx = @(x) feval (A, x, "transp");
73 else
74 error (["bicg: first argument is expected to " ...
75 "be a function or a square matrix"]);
76 endif
77
78 if ((nargin < 3) || (isempty (tol)))
79 tol = 1e-6;
80 endif
81
82 if ((nargin < 4) || (isempty (maxit)))
83 maxit = min (rows (b), 20);
84 endif
85
86 if ((nargin < 5) || isempty (M1))
87 M1m1x = @(x, ignore) x;
88 M1tm1x = M1m1x;
89 elseif (ischar (M1))
90 fun = str2func (M1);
91 M1m1x = @(x) feval (fun, x, "notransp");
92 M1tm1x = @(x) feval (fun, x, "transp");
93 elseif (ismatrix (M1))
94 M1m1x = @(x) M1 \ x;
95 M1tm1x = @(x) M1' \ x;
96 elseif (isa (M1, "function_handle"))
97 M1m1x = @(x) feval (M1, x, "notransp");
98 M1tm1x = @(x) feval (M1, x, "transp");
99 else
100 error (["bicg: preconditioner is expected to " ...
101 "be a function or matrix"]);
102 endif
103
104 if ((nargin < 6) || isempty (M2))
105 M2m1x = @(x, ignore) x;
106 M2tm1x = M2m1x;
107 elseif (ischar (M2))
108 fun = str2func (M2);
109 M2m1x = @(x) feval (fun, x, "notransp");
110 M2tm1x = @(x) feval (fun, x, "transp");
111 elseif (ismatrix (M2))
112 M2m1x = @(x) M2 \ x;
113 M2tm1x = @(x) M2' \ x;
114 elseif (isa (M2, "function_handle"))
115 M2m1x = @(x) feval (M2, x, "notransp");
116 M2tm1x = @(x) feval (M2, x, "transp");
117 else
118 error (["bicg: preconditioner is expected to " ...
119 "be a function or matrix"]);
120 endif
121
122 Pm1x = @(x) M2m1x (M1m1x (x));
123 Ptm1x = @(x) M1tm1x (M2tm1x (x));
124
125 if ((nargin < 7) || (isempty (x0)))
126 x0 = zeros (size (b));
127 endif
128
129 y = x = x0;
130 c = b;
131
132 r0 = b - Ax (x);
133 s0 = c - Atx (y);
134
135 d = Pm1x (r0);
136 f = Ptm1x (s0);
137
138 bnorm = norm (b);
139 res0 = Inf;
140
141 if (any (r0 != 0))
142
143 for k = 1:maxit
144
145 a = (s0' * Pm1x (r0)) ./ (f' * Ax (d));
146
147 x += a * d;
148 y += conj (a) * f;
149
150 r1 = r0 - a * Ax (d);
151 s1 = s0 - conj (a) * Atx (f);
152
153 beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0));
154
155 d = Pm1x (r1) + beta * d;
156 f = Ptm1x (s1) + conj (beta) * f;
157
158 r0 = r1;
159 s0 = s1;
160
161 res1 = norm (b - Ax (x)) / bnorm;
162 if (res1 < tol)
163 flag = 0;
164 if (nargout < 2)
165 printf ("bicg converged at iteration %i ", k);
166 printf ("to a solution with relative residual %e\n", res1);
167 endif
168 break;
169 endif
170
171 if (res0 <= res1)
172 flag = 3;
173 printf ("bicg stopped at iteration %i ", k);
174 printf ("without converging to the desired tolerance %e\n", tol);
175 printf ("because the method stagnated.\n");
176 printf ("The iterate returned (number %i) ", k-1);
177 printf ("has relative residual %e\n", res0);
178 break
179 endif
180 res0 = res1;
181 if (nargout > 4)
182 resvec(k) = res0;
183 endif
184 endfor
185
186 if (k == maxit)
187 flag = 1;
188 printf ("bicg stopped at iteration %i ", maxit);
189 printf ("without converging to the desired tolerance %e\n", tol);
190 printf ("because the maximum number of iterations was reached. ");
191 printf ("The iterate returned (number %i) has ", maxit);
192 printf ("relative residual %e\n", res1);
193 endif
194
195 else
196 flag = 0;
197 if (nargout < 2)
198 printf ("bicg converged after 0 interations\n");
199 endif
200 endif
201
202 else
203 print_usage ();
204 endif
205
206 endfunction;
207
208
209 %!test
210 %! n = 100;
211 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
212 %! b = sum (A, 2);
213 %! tol = 1e-8;
214 %! maxit = 15;
215 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
216 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
217 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2);
218 %! assert (x, ones (size (b)), 1e-7);
219 %!
220 %!test
221 %! n = 100;
222 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
223 %! b = sum (A, 2);
224 %! tol = 1e-8;
225 %! maxit = 15;
226 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
227 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
228 %!
229 %! function y = afun (x, t, a)
230 %! switch t
231 %! case "notransp"
232 %! y = a * x;
233 %! case "transp"
234 %! y = a' * x;
235 %! endswitch
236 %! endfunction
237 %!
238 %! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A),
239 %! b, tol, maxit, M1, M2);
240 %! assert (x, ones (size (b)), 1e-7);
241
242 %!test
243 %! n = 100;
244 %! tol = 1e-8;
245 %! a = sprand (n, n, .1);
246 %! A = a' * a + 100 * eye (n);
247 %! b = sum (A, 2);
248 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A)));
249 %! assert (x, ones (size (b)), 1e-7);