Mercurial > octave-antonio

comparison scripts/sparse/bicg.m @ 13141:e81ddf9cacd5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: untabify and remove trailing whitespace from source files * bicg.m, gmres.m, pkg.m: Untabify and remove trailing whitespace. * libcruft/Makefile.am, libcruft/blas-xtra/cdotc3.f, libcruft/blas-xtra/cmatm3.f, libcruft/blas-xtra/ddot3.f, libcruft/blas-xtra/dmatm3.f, libcruft/blas-xtra/sdot3.f, libcruft/blas-xtra/smatm3.f, libcruft/blas-xtra/zdotc3.f, libcruft/blas-xtra/zmatm3.f, libcruft/lapack-xtra/crsf2csf.f, libcruft/lapack-xtra/zrsf2csf.f, liboctave/Array.cc, liboctave/DASPK-opts.in, liboctave/DASRT-opts.in, liboctave/DASSL-opts.in, liboctave/LSODE-opts.in, liboctave/Makefile.a,mliboctave/Quad-opts.in, liboctave/Sparse-perm-op-defs.h, scripts/Makefile.a,mscripts/deprecated/glpkmex.m, scripts/general/blkdiag.m, scripts/general/interp1.m, scripts/general/profshow.m, scripts/general/quadl.m, scripts/general/triplequad.m, scripts/help/__makeinfo__.m, scripts/io/strread.m, scripts/io/textread.m, scripts/io/textscan.m, scripts/linear-algebra/rank.m, scripts/miscellaneous/gzip.m, scripts/miscellaneous/private/__xzip__.m, scripts/miscellaneous/tempdir.m, scripts/miscellaneous/unpack.m, scripts/pkg/pkg.m, scripts/plot/allchild.m, scripts/plot/ancestor.m, scripts/plot/cla.m, scripts/plot/clf.m, scripts/plot/findall.m, scripts/plot/findobj.m, scripts/plot/gca.m, scripts/plot/gcf.m, scripts/plot/hggroup.m, scripts/plot/isfigure.m, scripts/plot/ishghandle.m, scripts/plot/legend.m, scripts/plot/line.m, scripts/plot/loglog.m, scripts/plot/patch.m, scripts/plot/print.m, scripts/plot/private/__quiver__.m, scripts/plot/private/__scatter__.m, scripts/plot/rectangle.m, scripts/plot/semilogx.m, scripts/plot/semilogy.m, scripts/plot/surface.m, scripts/plot/text.m, scripts/plot/title.m, scripts/plot/trisurf.m, scripts/plot/view.m, scripts/plot/whitebg.m, scripts/plot/xlabel.m, scripts/plot/xlim.m, scripts/plot/ylabel.m, scripts/plot/ylim.m, scripts/plot/zlabel.m, scripts/plot/zlim.m, scripts/polynomial/mkpp.m, scripts/polynomial/polygcd.m, scripts/polynomial/ppint.m, scripts/polynomial/ppjumps.m, scripts/polynomial/ppval.m, scripts/set/setxor.m, scripts/sparse/bicgstab.m, scripts/sparse/cgs.m, scripts/sparse/spconvert.m, scripts/specfun/nthroot.m, scripts/strings/strmatch.m, scripts/strings/untabify.m, scripts/testfun/demo.m, scripts/testfun/example.m, src/DLD-FUNCTIONS/filter.cc, src/DLD-FUNCTIONS/mgorth.cc, src/DLD-FUNCTIONS/quadcc.cc, src/DLD-FUNCTIONS/str2double.cc, src/Makefile.a,msrc/gl-render.cc, src/gl2ps-renderer.cc, src/graphics.cc, src/octave-config.cc.in, src/octave-config.in, src/ov-class.h, src/ov-fcn.h, src/profiler.cc, src/profiler.h, src/pt-binop.cc, src/pt-unop.cc, src/symtab.cc, src/txt-eng-ft.cc: Remove trailing whitespace.
author John W. Eaton <jwe@octave.org>
date Thu, 15 Sep 2011 12:51:10 -0400
parents e5aaba072d2b
children 53674ceb9133
comparison
equal deleted inserted replaced
13140:98d23b0f16e1 13141:e81ddf9cacd5
13 ## 13 ##
14 ## You should have received a copy of the GNU General Public License 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/>. 15 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
16 16
17 ## -*- texinfo -*- 17 ## -*- texinfo -*-
18 ## 18 ##
19 ## @deftypefn {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) 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}) 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}, ...) 21 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, ...)
22 ## Solve @code{A x = b} using the Bi-conjugate gradient iterative method. 22 ## Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
23 ## 23 ##
24 ## @itemize @minus 24 ## @itemize @minus
25 ## @item @var{rtol} is the relative tolerance, if not given 25 ## @item @var{rtol} is the relative tolerance, if not given
26 ## or set to [] the default value 1e-6 is used. 26 ## or set to [] the default value 1e-6 is used.
27 ## @item @var{maxit} the maximum number of outer iterations, 27 ## @item @var{maxit} the maximum number of outer iterations,
28 ## if not given or set to [] the default value 28 ## if not given or set to [] the default value
29 ## @code{min (20, numel (b))} is used. 29 ## @code{min (20, numel (b))} is used.
30 ## @item @var{x0} the initial guess, if not given or set to [] 30 ## @item @var{x0} the initial guess, if not given or set to []
31 ## the default value @code{zeros (size (b))} is used. 31 ## the default value @code{zeros (size (b))} is used.
32 ## @end itemize 32 ## @end itemize
33 ## 33 ##
34 ## @var{A} can be passed as a matrix or as a function handle or 34 ## @var{A} can be passed as a matrix or as a function handle or
35 ## inline function @code{f} such that @code{f(x, "notransp") = A*x} 35 ## inline function @code{f} such that @code{f(x, "notransp") = A*x}
36 ## and @code{f(x, "transp") = A'*x}. 36 ## and @code{f(x, "transp") = A'*x}.
37 ## 37 ##
38 ## The preconditioner @var{P} is given as @code{P = M1 * M2}. 38 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
39 ## Both @var{M1} and @var{M2} can be passed as a matrix or as 39 ## Both @var{M1} and @var{M2} can be passed as a matrix or as
40 ## a function handle or inline function @code{g} such that 40 ## a function handle or inline function @code{g} such that
41 ## @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and 41 ## @code{g(x, 'notransp') = M1 \ x} or @code{g(x, 'notransp') = M2 \ x} and
42 ## @code{g(x, 'transp') = M1' \ x} or @code{g(x, 'transp') = M2' \ x}. 42 ## @code{g(x, 'transp') = M1' \ x} or @code{g(x, 'transp') = M2' \ x}.
43 ## 43 ##
44 ## If colled with more than one output parameter 44 ## If colled with more than one output parameter
45 ## 45 ##
46 ## @itemize @minus 46 ## @itemize @minus
50 ## @item 1: the maximum number of iterations was reached before convergence 50 ## @item 1: the maximum number of iterations was reached before convergence
51 ## @item 3: the algorithm reached stagnation 51 ## @item 3: the algorithm reached stagnation
52 ## @end itemize 52 ## @end itemize
53 ## (the value 2 is unused but skipped for compatibility). 53 ## (the value 2 is unused but skipped for compatibility).
54 ## @item @var{relres} is the final value of the relative residual. 54 ## @item @var{relres} is the final value of the relative residual.
55 ## @item @var{iter} is the number of iterations performed. 55 ## @item @var{iter} is the number of iterations performed.
56 ## @item @var{resvec} is a vector containing the relative residual at each iteration. 56 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
57 ## @end itemize 57 ## @end itemize
58 ## 58 ##
59 ## @seealso{bicgstab,cgs,gmres,pcg} 59 ## @seealso{bicgstab,cgs,gmres,pcg}
60 ## 60 ##
62 62
63 63
64 function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0) 64 function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0)
65 65
66 if (nargin >= 2 && isvector (full (b))) 66 if (nargin >= 2 && isvector (full (b)))
67 67
68 if (ischar (A)) 68 if (ischar (A))
69 fun = str2func (A); 69 fun = str2func (A);
70 Ax = @(x) feval (fun, x, "notransp"); 70 Ax = @(x) feval (fun, x, "notransp");
71 Atx = @(x) feval (fun, x, "transp"); 71 Atx = @(x) feval (fun, x, "transp");
72 elseif (ismatrix (A)) 72 elseif (ismatrix (A))
103 M1tm1x = @(x) feval (M1, x, "transp"); 103 M1tm1x = @(x) feval (M1, x, "transp");
104 else 104 else
105 error (["bicg: preconditioner is expected to " ... 105 error (["bicg: preconditioner is expected to " ...
106 "be a function or matrix"]); 106 "be a function or matrix"]);
107 endif 107 endif
108 108
109 if (nargin < 6 || isempty (M2)) 109 if (nargin < 6 || isempty (M2))
110 M2m1x = @(x, ignore) x; 110 M2m1x = @(x, ignore) x;
111 M2tm1x = M2m1x; 111 M2tm1x = M2m1x;
112 elseif (ischar (M2)) 112 elseif (ischar (M2))
113 fun = str2func (M2); 113 fun = str2func (M2);
142 142
143 bnorm = norm (b); 143 bnorm = norm (b);
144 res0 = Inf; 144 res0 = Inf;
145 145
146 if (any (r0 != 0)) 146 if (any (r0 != 0))
147 147
148 for k = 1:maxit 148 for k = 1:maxit
149 149
150 a = (s0' * Pm1x (r0)) ./ (f' * Ax (d)); 150 a = (s0' * Pm1x (r0)) ./ (f' * Ax (d));
151 151
152 x += a * d; 152 x += a * d;
153 y += conj (a) * f; 153 y += conj (a) * f;
154 154
155 r1 = r0 - a * Ax (d); 155 r1 = r0 - a * Ax (d);
156 s1 = s0 - conj (a) * Atx (f); 156 s1 = s0 - conj (a) * Atx (f);
157 157
158 beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0)); 158 beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0));
159 159
160 d = Pm1x (r1) + beta * d; 160 d = Pm1x (r1) + beta * d;
161 f = Ptm1x (s1) + conj (beta) * f; 161 f = Ptm1x (s1) + conj (beta) * f;
162 162
163 r0 = r1; 163 r0 = r1;
164 s0 = s1; 164 s0 = s1;
165 165
166 res1 = norm (b - Ax (x)) / bnorm; 166 res1 = norm (b - Ax (x)) / bnorm;
167 if (res1 < tol) 167 if (res1 < tol)
168 flag = 0; 168 flag = 0;
169 if (nargout < 2) 169 if (nargout < 2)
170 printf ("bicg converged at iteration %i ", k); 170 printf ("bicg converged at iteration %i ", k);
171 printf ("to a solution with relative residual %e\n", res1); 171 printf ("to a solution with relative residual %e\n", res1);
172 endif 172 endif
173 break; 173 break;
174 endif 174 endif
175 175
176 if (res0 <= res1) 176 if (res0 <= res1)
177 flag = 3; 177 flag = 3;
178 printf ("bicg stopped at iteration %i ", k); 178 printf ("bicg stopped at iteration %i ", k);
179 printf ("without converging to the desired tolerance %e\n", tol); 179 printf ("without converging to the desired tolerance %e\n", tol);
180 printf ("because the method stagnated.\n"); 180 printf ("because the method stagnated.\n");
181 printf ("The iterate returned (number %i) ", k-1); 181 printf ("The iterate returned (number %i) ", k-1);
182 printf ("has relative residual %e\n", res0); 182 printf ("has relative residual %e\n", res0);
183 break 183 break
184 endif 184 endif
185 res0 = res1; 185 res0 = res1;
186 if (nargout > 4) 186 if (nargout > 4)
187 resvec(k) = res0; 187 resvec(k) = res0;
188 endif 188 endif
189 endfor 189 endfor
190 190
191 if (k == maxit) 191 if (k == maxit)
192 flag = 1; 192 flag = 1;
193 printf ("bicg stopped at iteration %i ", maxit); 193 printf ("bicg stopped at iteration %i ", maxit);
194 printf ("without converging to the desired tolerance %e\n", tol); 194 printf ("without converging to the desired tolerance %e\n", tol);
195 printf ("because the maximum number of iterations was reached. "); 195 printf ("because the maximum number of iterations was reached. ");
196 printf ("The iterate returned (number %i) has ", maxit); 196 printf ("The iterate returned (number %i) has ", maxit);
197 printf ("relative residual %e\n", res1); 197 printf ("relative residual %e\n", res1);
198 endif 198 endif
201 flag = 0; 201 flag = 0;
202 if (nargout < 2) 202 if (nargout < 2)
203 printf ("bicg converged after 0 interations\n"); 203 printf ("bicg converged after 0 interations\n");
204 endif 204 endif
205 endif 205 endif
206 206
207 else 207 else
208 print_usage (); 208 print_usage ();
209 endif 209 endif
210 210
211 endfunction; 211 endfunction;
212 212
213 213
214 %!test 214 %!test
215 %! n = 100; 215 %! n = 100;
216 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); 216 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
217 %! b = sum (A, 2); 217 %! b = sum (A, 2);
218 %! tol = 1e-8; 218 %! tol = 1e-8;
219 %! maxit = 15; 219 %! maxit = 15;
220 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); 220 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
221 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); 221 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
222 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2); 222 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2);
223 %! assert (x, ones (size (b)), 1e-7); 223 %! assert (x, ones (size (b)), 1e-7);
224 %! 224 %!
225 %!test 225 %!test
226 %! n = 100; 226 %! n = 100;
227 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); 227 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
228 %! b = sum (A, 2); 228 %! b = sum (A, 2);
229 %! tol = 1e-8; 229 %! tol = 1e-8;
230 %! maxit = 15; 230 %! maxit = 15;
231 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); 231 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
232 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); 232 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
233 %! 233 %!
234 %! function y = afun (x, t, a) 234 %! function y = afun (x, t, a)
235 %! switch t 235 %! switch t
236 %! case "notransp" 236 %! case "notransp"
237 %! y = a * x; 237 %! y = a * x;
238 %! case "transp" 238 %! case "transp"
239 %! y = a' * x; 239 %! y = a' * x;
240 %! endswitch 240 %! endswitch
241 %! endfunction 241 %! endfunction
242 %! 242 %!
243 %! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A), 243 %! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A),
244 %! b, tol, maxit, M1, M2); 244 %! b, tol, maxit, M1, M2);
245 %! assert (x, ones (size (b)), 1e-7); 245 %! assert (x, ones (size (b)), 1e-7);
246 246
247 %!test 247 %!test
248 %! n = 100; 248 %! n = 100;
249 %! tol = 1e-8; 249 %! tol = 1e-8;
250 %! a = sprand (n, n, .1); 250 %! a = sprand (n, n, .1);
251 %! A = a' * a + 100 * eye (n); 251 %! A = a' * a + 100 * eye (n);
252 %! b = sum (A, 2); 252 %! b = sum (A, 2);
253 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A))); 253 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A)));
254 %! assert (x, ones (size (b)), 1e-7); 254 %! assert (x, ones (size (b)), 1e-7);