Mercurial > hg > octave-lyh

comparison scripts/sparse/bicgstab.m @ 13091:e5aaba072d2b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
13090:7f127e079a7c 13091:e5aaba072d2b
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## 21 ##
22 ## @deftypefn {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) 22 ## @deftypefn {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
23 ## @deftypefnx {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P}) 23 ## @deftypefnx {Function File} {@var{x} =} bicgstab (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
24 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, ...) 24 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicgstab (@var{A}, @var{b}, ...)
25 ## 25 ## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative method.
26 ## Solve @code{A x = b} using the stabilizied Bi-conjugate gradient iterative method. 26 ##
27 ## 27 ## @itemize @minus
28 ## @itemize @minus 28 ## @item @var{rtol} is the relative tolerance, if not given or set to
29 ## 29 ## [] the default value 1e-6 is used.
30 ## @item @var{rtol} is the relative tolerance, if not given or set to 30 ## @item @var{maxit} the maximum number of outer iterations, if not
31 ## [] the default value 1e-6 is used. 31 ## given or set to [] the default value @code{min (20, numel (b))} is
32 ## 32 ## used.
33 ## @item @var{maxit} the maximum number of outer iterations, if not 33 ## @item @var{x0} the initial guess, if not given or set to [] the
34 ## given or set to [] the default value @code{min (20, numel (b))} is 34 ## default value @code{zeros (size (b))} is used.
35 ## used. 35 ## @end itemize
36 ## 36 ##
37 ## @item @var{x0} the initial guess, if not given or set to [] the 37 ## @var{A} can be passed as a matrix or as a function handle or
38 ## default value @code{zeros (size (b))} is used. 38 ## inline function @code{f} such that @code{f(x) = A*x}.
39 ## 39 ##
40 ## @end itemize 40 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
41 ## 41 ## Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or
42 ## @var{A} can be passed as a matrix or as a function handle or 42 ## inline function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}.
43 ## inline function @code{f} such that @code{f(x) = A*x}. 43 ##
44 ## 44 ## If called with more than one output parameter
45 ## The preconditioner @var{P} is given as @code{P = M1 * M2}. 45 ##
46 ## Both @var{M1} and @var{M2} can be passed as a matrix or as a function handle or 46 ## @itemize @minus
47 ## inline function @code{g} such that @code{g(x) = M1 \ x} or @code{g(x) = M2 \ x}. 47 ## @item @var{flag} indicates the exit status:
48 ## 48 ## @itemize @minus
49 ## If called with more than one output parameter 49 ## @item 0: iteration converged to the within the chosen tolerance
50 ## 50 ## @item 1: the maximum number of iterations was reached before convergence
51 ## @itemize @minus 51 ## @item 3: the algorithm reached stagnation
52 ## @item @var{flag} indicates the exit status: 52 ## @end itemize
53 ## @itemize @minus 53 ## (the value 2 is unused but skipped for compatibility).
54 ## @item 0: iteration converged to the within the chosen tolerance 54 ## @item @var{relres} is the final value of the relative residual.
55 ## @item 1: the maximum number of iterations was reached before convergence 55 ## @item @var{iter} is the number of iterations performed.
56 ## @item 3: the algorithm reached stagnation 56 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
57 ## @end itemize 57 ## @end itemize
58 ## (the value 2 is unused but skipped for compatibility). 58 ##
59 ## @item @var{relres} is the final value of the relative residual. 59 ## @seealso{bicg,cgs,gmres,pcg}
60 ## @item @var{iter} is the number of iterations performed.
61 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
62 ## @end itemize
63 ##
64 ## @seealso{pcg,cgs,bicg,gmres}
65 ## 60 ##
66 ## @end deftypefn 61 ## @end deftypefn
67 62
68 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, 63 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit,
69 M1, M2, x0) 64 M1, M2, x0)
70 65
71 if ((nargin >= 2) && (nargin <= 7) && isvector (full (b))) 66 if (nargin >= 2 && nargin <= 7 && isvector (full (b)))
72 67
73 if (ischar (A)) 68 if (ischar (A))
74 A = str2func (A); 69 A = str2func (A);
75 elseif (ismatrix (A)) 70 elseif (ismatrix (A))
76 Ax = @(x) A * x; 71 Ax = @(x) A * x;
79 else 74 else
80 error (["bicgstab: first argument is expected " ... 75 error (["bicgstab: first argument is expected " ...
81 "to be a function or a square matrix"]); 76 "to be a function or a square matrix"]);
82 endif 77 endif
83 78
84 if ((nargin < 3) || (isempty (tol))) 79 if (nargin < 3 || isempty (tol))
85 tol = 1e-6; 80 tol = 1e-6;
86 endif 81 endif
87 82
88 if ((nargin < 4) || (isempty (maxit))) 83 if (nargin < 4 || isempty (maxit))
89 maxit = min (rows (b), 20); 84 maxit = min (rows (b), 20);
90 endif 85 endif
91 86
92 if ((nargin < 5) || isempty (M1)) 87 if (nargin < 5 || isempty (M1))
93 M1m1x = @(x) x; 88 M1m1x = @(x) x;
94 elseif (ischar (M1)) 89 elseif (ischar (M1))
95 M1m1x = str2func (M1); 90 M1m1x = str2func (M1);
96 elseif (ismatrix (M1)) 91 elseif (ismatrix (M1))
97 M1m1x = @(x) M1 \ x; 92 M1m1x = @(x) M1 \ x;
98 elseif (isa (M1, "function_handle")) 93 elseif (isa (M1, "function_handle"))
99 M1m1x = @(x) feval (M1, x); 94 M1m1x = @(x) feval (M1, x);
100 else 95 else
101 error (["bicgstab: preconditioner is " ... 96 error (["bicgstab: preconditioner is " ...
102 "expected to be a function or matrix"]); 97 "expected to be a function or matrix"]);
103 endif 98 endif
104 99
105 if ((nargin < 6) || isempty (M2)) 100 if (nargin < 6 || isempty (M2))
106 M2m1x = @(x) x; 101 M2m1x = @(x) x;
107 elseif (ischar (M2)) 102 elseif (ischar (M2))
108 M2m1x = str2func (M2); 103 M2m1x = str2func (M2);
109 elseif (ismatrix (M2)) 104 elseif (ismatrix (M2))
110 M2m1x = @(x) M2 \ x; 105 M2m1x = @(x) M2 \ x;
111 elseif (isa (M2, "function_handle")) 106 elseif (isa (M2, "function_handle"))
112 M2m1x = @(x) feval (M2, x); 107 M2m1x = @(x) feval (M2, x);
113 else 108 else
114 error (["bicgstab: preconditioner is "... 109 error (["bicgstab: preconditioner is "...
115 "expected to be a function or matrix"]); 110 "expected to be a function or matrix"]);
116 endif 111 endif
117 112
118 precon = @(x) M2m1x (M1m1x (x)); 113 precon = @(x) M2m1x (M1m1x (x));
119 114
120 if ((nargin < 7) || (isempty (x0))) 115 if (nargin < 7 || isempty (x0))
121 x0 = zeros (size (b)); 116 x0 = zeros (size (b));
122 endif 117 endif
123 118
124
125 ## specifies initial estimate x0 119 ## specifies initial estimate x0
126 if (nargin < 7) 120 if (nargin < 7)
127 x = zeros (rows (b), 1); 121 x = zeros (rows (b), 1);
128 else 122 else
129 x = x0; 123 x = x0;