Mercurial > hg > octave-lyh

comparison scripts/sparse/gmres.m @ 14773:5c269b73f467

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gmres.m: Overhaul function. Fix bug #36568. * gmres.m: Return correct exit flag (bug #36568). Return relative residual, not norm, for Matlab compatibility. Use Octave coding conventions including same variable names in documentation as in function. Add %!demo and %!error tests.
author Rik <octave@nomad.inbox5.com>
date Fri, 15 Jun 2012 10:15:07 -0700
parents b76f0740940e
children 93e5ade3fda4
comparison
equal deleted inserted replaced
14772:39a2e91a246e 14773:5c269b73f467
30 ## @item @var{maxit} is the maximum number of outer iterations, 30 ## @item @var{maxit} is the maximum number of outer iterations,
31 ## if not given or set to [] the default value 31 ## if not given or set to [] the default value
32 ## @code{min (10, numel (b) / restart)} is used. 32 ## @code{min (10, numel (b) / restart)} is used.
33 ## 33 ##
34 ## @item @var{x0} is the initial guess, 34 ## @item @var{x0} is the initial guess,
35 ## if not given or set to [] the default value @code{zeros(size (b))} is used. 35 ## if not given or set to [] the default value @code{zeros (size (b))} is used.
36 ## 36 ##
37 ## @item @var{m} is the restart parameter, 37 ## @item @var{m} is the restart parameter,
38 ## if not given or set to [] the default value @code{numel (b)} is used. 38 ## if not given or set to [] the default value @code{numel (b)} is used.
39 ## @end itemize 39 ## @end itemize
40 ## 40 ##
55 ## 55 ##
56 ## @item 1 : maximum number of iterations exceeded 56 ## @item 1 : maximum number of iterations exceeded
57 ## 57 ##
58 ## @item 2 : unused, but skipped for compatibility 58 ## @item 2 : unused, but skipped for compatibility
59 ## 59 ##
60 ## @item 3 : algorithm reached stagnation 60 ## @item 3 : algorithm reached stagnation (no change between iterations)
61 ## @end table 61 ## @end table
62 ## 62 ##
63 ## @item @var{relres} is the final value of the relative residual. 63 ## @item @var{relres} is the final value of the relative residual.
64 ## 64 ##
65 ## @item @var{iter} is a vector containing the number of outer iterations and 65 ## @item @var{iter} is a vector containing the number of outer iterations and
70 ## @end itemize 70 ## @end itemize
71 ## 71 ##
72 ## @seealso{bicg, bicgstab, cgs, pcg} 72 ## @seealso{bicg, bicgstab, cgs, pcg}
73 ## @end deftypefn 73 ## @end deftypefn
74 74
75 function [x, flag, presn, it, resids] = gmres (A, b, restart, rtol, maxit, M1, M2, x0) 75 function [x, flag, relres, it, resvec] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)
76 76
77 if (nargin < 2 || nargin > 8) 77 if (nargin < 2 || nargin > 8)
78 print_usage (); 78 print_usage ();
79 endif 79 endif
80 80
140 H = zeros (restart + 1, restart); 140 H = zeros (restart + 1, restart);
141 141
142 ## begin loop 142 ## begin loop
143 iter = 1; 143 iter = 1;
144 restart_it = restart + 1; 144 restart_it = restart + 1;
145 resids = zeros (maxit, 1); 145 resvec = zeros (maxit, 1);
146 resids(1) = presn; 146 resvec(1) = presn;
147 prec_b_norm = norm (Pm1x (b), 2); 147 prec_b_norm = norm (Pm1x (b), 2);
148 flag = 1; 148 flag = 1; # Default flag is maximum # of iterations exceeded
149 149
150 while (iter <= maxit * restart && presn > rtol * prec_b_norm) 150 while (iter <= maxit * restart && presn > rtol * prec_b_norm)
151 151
152 ## restart 152 ## restart
153 if (restart_it > restart) 153 if (restart_it > restart)
172 172
173 presn = norm (little_res, 2); 173 presn = norm (little_res, 2);
174 174
175 x = x_old + V(:, 1:restart_it) * Y(1:restart_it); 175 x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
176 176
177 resids(iter) = presn; 177 resvec(iter) = presn;
178 if (norm (x - x_old, inf) <= eps) 178 if (norm (x - x_old, inf) <= eps)
179 flag = 3; 179 flag = 3; # Stagnation: no change between iterations
180 break 180 break;
181 endif 181 endif
182 182
183 restart_it++ ; 183 restart_it++ ;
184 iter++; 184 iter++;
185 endwhile 185 endwhile
186 186
187 if (presn > rtol * prec_b_norm) 187 if (nargout > 1)
188 flag = 0; 188 ## Calculate extra outputs as requested
189 endif 189 relres = presn / prec_b_norm;
190 190 if (relres <= rtol)
191 resids = resids(1:iter-1); 191 flag = 0; # Converged to solution within tolerance
192 it = [ceil(iter / restart), rem(iter, restart)]; 192 endif
193
194 resvec = resvec(1:iter-1);
195 it = [ceil(iter / restart), rem(iter, restart)];
196 endif
193 197
194 endfunction 198 endfunction
195 199
200
201 %!demo
202 %! dim = 20;
203 %! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
204 %! b = ones (dim, 1);
205 %! [x, flag, relres, iter, resvec] = gmres (A, b, 10, 1e-10, dim, @(x) x ./ diag (A), [], b)
196 206
197 %!shared A, b, dim 207 %!shared A, b, dim
198 %! dim = 100; 208 %! dim = 100;
199 %!test 209 %!test
200 %! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim); 210 %! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
208 %! 218 %!
209 %!test 219 %!test
210 %! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim); 220 %! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim);
211 %! A = A'*A; 221 %! A = A'*A;
212 %! b = rand (dim, 1); 222 %! b = rand (dim, 1);
213 %! [x, resids] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag (A), [], []); 223 %! [x, resvec] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag (A), [], []);
214 %! assert (x, A\b, 1e-9*norm (x, Inf)); 224 %! assert (x, A\b, 1e-9*norm (x, Inf));
215 %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag (diag (A)) \ x, [], []); 225 %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag (diag (A)) \ x, [], []);
216 %! assert (x, A\b, 1e-9*norm (x, Inf)); 226 %! assert (x, A\b, 1e-9*norm (x, Inf));
217 %!test 227 %!test
218 %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [], []); 228 %! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [], []);
219 %! assert (x, A\b, 1e-7*norm (x, Inf)); 229 %! assert (x, A\b, 1e-7*norm (x, Inf));
220 230
231
232 %!error gmres (1)
233 %!error gmres (1,2,3,4,5,6,7,8,9)
234 %!error <A must be> gmres ({1},2)
235 %!error <A must be a function or matrix> gmres ({1},2)
236 %!error <M1 must be a function or matrix> gmres (1,2,3,4,5,{6})
237 %!error <M2 must be a function or matrix> gmres (1,2,3,4,5,6,{7})