|
5837
|
1 ## Copyright (C) 2004 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
7016
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
5837
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
5837
|
18 |
|
|
19 ## -*- texinfo -*- |
|
5838
|
20 ## @deftypefn {Function File} {@var{x} =} pcr (@var{a}, @var{b}, @var{tol}, @var{maxit}, @var{m}, @var{x0}, @dots{}) |
|
5837
|
21 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} pcr (@dots{}) |
|
|
22 ## |
|
5838
|
23 ## Solves the linear system of equations @code{@var{a} * @var{x} = |
|
5837
|
24 ## @var{b}} by means of the Preconditioned Conjugate Residuals iterative |
|
|
25 ## method. The input arguments are |
|
|
26 ## |
|
|
27 ## @itemize |
|
|
28 ## @item |
|
5838
|
29 ## @var{a} can be either a square (preferably sparse) matrix or a |
|
5837
|
30 ## function handle, inline function or string containing the name |
|
5838
|
31 ## of a function which computes @code{@var{a} * @var{x}}. In principle |
|
|
32 ## @var{a} should be symmetric and non-singular; if @code{pcr} |
|
|
33 ## finds @var{a} to be numerically singular, you will get a warning |
|
5837
|
34 ## message and the @var{flag} output parameter will be set. |
|
|
35 ## |
|
|
36 ## @item |
|
|
37 ## @var{b} is the right hand side vector. |
|
|
38 ## |
|
|
39 ## @item |
|
|
40 ## @var{tol} is the required relative tolerance for the residual error, |
|
5838
|
41 ## @code{@var{b} - @var{a} * @var{x}}. The iteration stops if @code{norm |
|
|
42 ## (@var{b} - @var{a} * @var{x}) <= @var{tol} * norm (@var{b} - @var{a} * |
|
5837
|
43 ## @var{x0})}. If @var{tol} is empty or is omitted, the function sets |
|
|
44 ## @code{@var{tol} = 1e-6} by default. |
|
|
45 ## |
|
|
46 ## @item |
|
|
47 ## @var{maxit} is the maximum allowable number of iterations; if |
|
|
48 ## @code{[]} is supplied for @code{maxit}, or @code{pcr} has less |
|
|
49 ## arguments, a default value equal to 20 is used. |
|
|
50 ## |
|
|
51 ## @item |
|
5838
|
52 ## @var{m} is the (left) preconditioning matrix, so that the iteration is |
|
5837
|
53 ## (theoretically) equivalent to solving by @code{pcr} @code{@var{P} * |
|
5838
|
54 ## @var{x} = @var{m} \ @var{b}}, with @code{@var{P} = @var{m} \ @var{a}}. |
|
5837
|
55 ## Note that a proper choice of the preconditioner may dramatically |
|
|
56 ## improve the overall performance of the method. Instead of matrix |
|
5838
|
57 ## @var{m}, the user may pass a function which returns the results of |
|
|
58 ## applying the inverse of @var{m} to a vector (usually this is the |
|
5837
|
59 ## preferred way of using the preconditioner). If @code{[]} is supplied |
|
5838
|
60 ## for @var{m}, or @var{m} is omitted, no preconditioning is applied. |
|
5837
|
61 ## |
|
|
62 ## @item |
|
|
63 ## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the |
|
|
64 ## function sets @var{x0} to a zero vector by default. |
|
|
65 ## @end itemize |
|
|
66 ## |
|
|
67 ## The arguments which follow @var{x0} are treated as parameters, and |
|
5838
|
68 ## passed in a proper way to any of the functions (@var{a} or @var{m}) |
|
5837
|
69 ## which are passed to @code{pcr}. See the examples below for further |
|
|
70 ## details. The output arguments are |
|
|
71 ## |
|
|
72 ## @itemize |
|
|
73 ## @item |
|
|
74 ## @var{x} is the computed approximation to the solution of |
|
5838
|
75 ## @code{@var{a} * @var{x} = @var{b}}. |
|
5837
|
76 ## |
|
|
77 ## @item |
|
|
78 ## @var{flag} reports on the convergence. @code{@var{flag} = 0} means |
|
|
79 ## the solution converged and the tolerance criterion given by @var{tol} |
|
|
80 ## is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit |
|
|
81 ## for the iteration count was reached. @code{@var{flag} = 3} reports t |
|
|
82 ## @code{pcr} breakdown, see [1] for details. |
|
|
83 ## |
|
|
84 ## @item |
|
|
85 ## @var{relres} is the ratio of the final residual to its initial value, |
|
|
86 ## measured in the Euclidean norm. |
|
|
87 ## |
|
|
88 ## @item |
|
|
89 ## @var{iter} is the actual number of iterations performed. |
|
|
90 ## |
|
|
91 ## @item |
|
|
92 ## @var{resvec} describes the convergence history of the method, |
|
|
93 ## so that @code{@var{resvec} (i)} contains the Euclidean norms of the |
|
7001
|
94 ## residual after the (@var{i}-1)-th iteration, @code{@var{i} = |
|
5838
|
95 ## 1,2, @dots{}, @var{iter}+1}. |
|
5837
|
96 ## @end itemize |
|
|
97 ## |
|
|
98 ## Let us consider a trivial problem with a diagonal matrix (we exploit the |
|
|
99 ## sparsity of A) |
|
|
100 ## |
|
|
101 ## @example |
|
|
102 ## @group |
|
|
103 ## N = 10; |
|
|
104 ## A = diag([1:N]); A = sparse(A); |
|
|
105 ## b = rand(N,1); |
|
|
106 ## @end group |
|
|
107 ## @end example |
|
|
108 ## |
|
|
109 ## @sc{Example 1:} Simplest use of @code{pcr} |
|
|
110 ## |
|
|
111 ## @example |
|
|
112 ## x = pcr(A, b) |
|
|
113 ## @end example |
|
|
114 ## |
|
|
115 ## @sc{Example 2:} @code{pcr} with a function which computes |
|
5838
|
116 ## @code{@var{a} * @var{x}}. |
|
5837
|
117 ## |
|
|
118 ## @example |
|
|
119 ## @group |
|
|
120 ## function y = applyA(x) |
|
|
121 ## y = [1:10]'.*x; |
|
|
122 ## endfunction |
|
|
123 ## |
|
|
124 ## x = pcr('applyA',b) |
|
|
125 ## @end group |
|
|
126 ## @end example |
|
|
127 ## |
|
|
128 ## @sc{Example 3:} Preconditioned iteration, with full diagnostics. The |
|
|
129 ## preconditioner (quite strange, because even the original matrix |
|
5838
|
130 ## @var{a} is trivial) is defined as a function |
|
5837
|
131 ## |
|
|
132 ## @example |
|
|
133 ## @group |
|
|
134 ## function y = applyM(x) |
|
|
135 ## K = floor(length(x)-2); |
|
|
136 ## y = x; |
|
|
137 ## y(1:K) = x(1:K)./[1:K]'; |
|
|
138 ## endfunction |
|
|
139 ## |
|
|
140 ## [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM') |
|
|
141 ## semilogy([1:iter+1], resvec); |
|
|
142 ## @end group |
|
|
143 ## @end example |
|
|
144 ## |
|
|
145 ## @sc{Example 4:} Finally, a preconditioner which depends on a |
|
|
146 ## parameter @var{k}. |
|
|
147 ## |
|
|
148 ## @example |
|
|
149 ## @group |
|
|
150 ## function y = applyM(x, varargin) |
|
|
151 ## K = varargin@{1@}; |
|
|
152 ## y = x; y(1:K) = x(1:K)./[1:K]'; |
|
|
153 ## endfunction |
|
|
154 ## |
|
|
155 ## [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM',[],3) |
|
|
156 ## @end group |
|
|
157 ## @end example |
|
|
158 ## |
|
|
159 ## @sc{References} |
|
|
160 ## |
|
|
161 ## [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of |
|
|
162 ## Equations", section 9.5.4; Springer, 1994 |
|
|
163 ## |
|
|
164 ## @seealso{sparse, pcg} |
|
|
165 ## @end deftypefn |
|
|
166 |
|
5838
|
167 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> |
|
5837
|
168 |
|
5838
|
169 function [x, flag, relres, iter, resvec] = pcr (A, b, tol, maxit, M, x0, varargin) |
|
5837
|
170 |
|
|
171 breakdown = false; |
|
|
172 |
|
5838
|
173 if (nargin < 6 || isempty (x0)) |
|
|
174 x = zeros (size (b)); |
|
5837
|
175 else |
|
|
176 x = x0; |
|
|
177 endif |
|
|
178 |
|
|
179 if (nargin < 5) |
|
|
180 M = []; |
|
|
181 endif |
|
|
182 |
|
5838
|
183 if (nargin < 4 || isempty (maxit)) |
|
5837
|
184 maxit = 20; |
|
|
185 endif |
|
|
186 |
|
5838
|
187 maxit += 2; |
|
5837
|
188 |
|
5838
|
189 if (nargin < 3 || isempty (tol)) |
|
5837
|
190 tol = 1e-6; |
|
|
191 endif |
|
|
192 |
|
|
193 if (nargin < 2) |
|
5838
|
194 print_usage (); |
|
5837
|
195 endif |
|
|
196 |
|
|
197 ## init |
|
5838
|
198 if (isnumeric (A)) # is A a matrix? |
|
|
199 r = b - A*x; |
|
5837
|
200 else # then A should be a function! |
|
5838
|
201 r = b - feval (A, x, varargin{:}); |
|
5837
|
202 endif |
|
|
203 |
|
5838
|
204 if (isnumeric (M)) # is M a matrix? |
|
|
205 if (isempty (M)) # if M is empty, use no precond |
|
5837
|
206 p = r; |
|
|
207 else # otherwise, apply the precond |
|
5838
|
208 p = M \ r; |
|
5837
|
209 endif |
|
|
210 else # then M should be a function! |
|
5838
|
211 p = feval (M, r, varargin{:}); |
|
5837
|
212 endif |
|
|
213 |
|
|
214 iter = 2; |
|
|
215 |
|
|
216 b_bot_old = 1; |
|
5838
|
217 q_old = p_old = s_old = zeros (size (x)); |
|
5837
|
218 |
|
5838
|
219 if (isnumeric (A)) # is A a matrix? |
|
|
220 q = A * p; |
|
5837
|
221 else # then A should be a function! |
|
5838
|
222 q = feval (A, p, varargin{:}); |
|
5837
|
223 endif |
|
|
224 |
|
5838
|
225 resvec(1) = abs (norm (r)); |
|
5837
|
226 |
|
|
227 ## iteration |
|
5838
|
228 while (resvec(iter-1) > tol*resvec(1) && iter < maxit) |
|
|
229 |
|
|
230 if (isnumeric (M)) # is M a matrix? |
|
|
231 if (isempty (M)) # if M is empty, use no precond |
|
5837
|
232 s = q; |
|
|
233 else # otherwise, apply the precond |
|
5838
|
234 s = M \ q; |
|
5837
|
235 endif |
|
|
236 else # then M should be a function! |
|
5838
|
237 s = feval (M, q, varargin{:}); |
|
5837
|
238 endif |
|
5838
|
239 b_top = r' * s; |
|
|
240 b_bot = q' * s; |
|
5837
|
241 |
|
|
242 if (b_bot == 0.0) |
|
|
243 breakdown = true; |
|
|
244 break; |
|
|
245 endif |
|
5838
|
246 lambda = b_top / b_bot; |
|
5837
|
247 |
|
5838
|
248 x += lambda*p; |
|
|
249 r -= lambda*q; |
|
5837
|
250 |
|
|
251 if (isnumeric(A)) # is A a matrix? |
|
|
252 t = A*s; |
|
|
253 else # then A should be a function! |
|
5838
|
254 t = feval (A, s, varargin{:}); |
|
5837
|
255 endif |
|
|
256 |
|
5838
|
257 alpha0 = (t'*s) / b_bot; |
|
|
258 alpha1 = (t'*s_old) / b_bot_old; |
|
5837
|
259 |
|
5838
|
260 p_temp = p; |
|
|
261 q_temp = q; |
|
|
262 |
|
5837
|
263 p = s - alpha0*p - alpha1*p_old; |
|
|
264 q = t - alpha0*q - alpha1*q_old; |
|
|
265 |
|
|
266 s_old = s; |
|
|
267 p_old = p_temp; |
|
|
268 q_old = q_temp; |
|
|
269 b_bot_old = b_bot; |
|
|
270 |
|
5838
|
271 resvec(iter) = abs (norm (r)); |
|
|
272 iter++; |
|
5837
|
273 endwhile |
|
|
274 |
|
|
275 flag = 0; |
|
5838
|
276 relres = resvec(iter-1) ./ resvec(1); |
|
|
277 iter -= 2; |
|
|
278 if (iter >= maxit-2) |
|
5837
|
279 flag = 1; |
|
|
280 if (nargout < 2) |
|
6498
|
281 warning ("pcr: maximum number of iterations (%d) reached\n", iter); |
|
|
282 warning ("the initial residual norm was reduced %g times.\n", 1.0/relres); |
|
5837
|
283 endif |
|
5838
|
284 elseif (nargout < 2 && ! breakdown) |
|
6498
|
285 fprintf (stderr, "pcr: converged in %d iterations. \n", iter); |
|
|
286 fprintf (stderr, "the initial residual norm was reduced %g times.\n", |
|
5838
|
287 1.0/relres); |
|
5837
|
288 endif |
|
5838
|
289 |
|
5837
|
290 if (breakdown) |
|
|
291 flag = 3; |
|
|
292 if (nargout < 2) |
|
7001
|
293 warning ("pcr: breakdown occurred:\n"); |
|
6498
|
294 warning ("system matrix singular or preconditioner indefinite?\n"); |
|
5837
|
295 endif |
|
|
296 endif |
|
|
297 |
|
|
298 endfunction |
|
|
299 |
|
|
300 %!demo |
|
|
301 %! |
|
|
302 %! # Simplest usage of PCR (see also 'help pcr') |
|
|
303 %! |
|
|
304 %! N = 20; |
|
|
305 %! A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution |
|
|
306 %! x = pcr(A,b); |
|
|
307 %! printf('The solution relative error is %g\n', norm(x-y)/norm(y)); |
|
|
308 %! |
|
|
309 %! # You shouldn't be afraid if PCR issues some warning messages in this |
|
|
310 %! # example: watch out in the second example, why it takes N iterations |
|
|
311 %! # of PCR to converge to (a very accurate, by the way) solution |
|
|
312 %!demo |
|
|
313 %! |
|
|
314 %! # Full output from PCR |
|
|
315 %! # We use this output to plot the convergence history |
|
|
316 %! |
|
|
317 %! N = 20; |
|
|
318 %! A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution |
|
|
319 %! [x, flag, relres, iter, resvec] = pcr(A,b); |
|
|
320 %! printf('The solution relative error is %g\n', norm(x-X)/norm(X)); |
|
|
321 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); |
|
|
322 %! semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;'); |
|
|
323 %!demo |
|
|
324 %! |
|
|
325 %! # Full output from PCR |
|
|
326 %! # We use indefinite matrix based on the Hilbert matrix, with one |
|
|
327 %! # strongly negative eigenvalue |
|
|
328 %! # Hilbert matrix is extremely ill conditioned, so is ours, |
|
|
329 %! # and that's why PCR WILL have problems |
|
|
330 %! |
|
|
331 %! N = 10; |
|
|
332 %! A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution |
|
|
333 %! printf('Condition number of A is %g\n', cond(A)); |
|
|
334 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],200); |
|
|
335 %! if (flag == 3) |
|
|
336 %! printf('PCR breakdown. System matrix is [close to] singular\n'); |
|
|
337 %! end |
|
|
338 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); |
|
|
339 %! semilogy([0:iter],resvec,'o-g;absolute residual;'); |
|
|
340 %!demo |
|
|
341 %! |
|
|
342 %! # Full output from PCR |
|
|
343 %! # We use an indefinite matrix based on the 1-D Laplacian matrix for A, |
|
|
344 %! # and here we have cond(A) = O(N^2) |
|
|
345 %! # That's the reason we need some preconditioner; here we take |
|
|
346 %! # a very simple and not powerful Jacobi preconditioner, |
|
|
347 %! # which is the diagonal of A |
|
|
348 %! |
|
|
349 %! # Note that we use here indefinite preconditioners! |
|
|
350 %! |
|
|
351 %! N = 100; |
|
|
352 %! A = zeros(N,N); |
|
|
353 %! for i=1:N-1 # form 1-D Laplacian matrix |
|
|
354 %! A(i:i+1,i:i+1) = [2 -1; -1 2]; |
|
|
355 %! endfor |
|
|
356 %! A = [A, zeros(size(A)); zeros(size(A)), -A]; |
|
|
357 %! b = rand(2*N,1); X = A\b; #X is the true solution |
|
|
358 %! maxit = 80; |
|
|
359 %! printf('System condition number is %g\n',cond(A)); |
|
|
360 %! # No preconditioner: the convergence is very slow! |
|
|
361 %! |
|
|
362 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit); |
|
|
363 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); |
|
|
364 %! semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;'); |
|
|
365 %! |
|
|
366 %! pause(1); |
|
|
367 %! # Test Jacobi preconditioner: it will not help much!!! |
|
|
368 %! |
|
|
369 %! M = diag(diag(A)); # Jacobi preconditioner |
|
|
370 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M); |
|
|
371 %! hold on; |
|
|
372 %! semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;'); |
|
|
373 %! |
|
|
374 %! pause(1); |
|
|
375 %! # Test nonoverlapping block Jacobi preconditioner: this one should give |
|
|
376 %! # some convergence speedup! |
|
|
377 %! |
|
|
378 %! M = zeros(N,N);k=4; |
|
|
379 %! for i=1:k:N # get k x k diagonal blocks of A |
|
|
380 %! M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1); |
|
|
381 %! endfor |
|
|
382 %! M = [M, zeros(size(M)); zeros(size(M)), -M]; |
|
|
383 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M); |
|
|
384 %! semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;'); |
|
|
385 %! hold off; |
|
|
386 %!test |
|
|
387 %! |
|
|
388 %! #solve small indefinite diagonal system |
|
|
389 %! |
|
|
390 %! N = 10; |
|
|
391 %! A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution |
|
|
392 %! [x, flag] = pcr(A,b,[],N+1); |
|
|
393 %! assert(norm(x-X)/norm(X)<1e-10); |
|
|
394 %! assert(flag,0); |
|
|
395 %! |
|
|
396 %!test |
|
|
397 %! |
|
|
398 %! #solve tridiagonal system, do not converge in default 20 iterations |
|
|
399 %! #should perform max allowable default number of iterations |
|
|
400 %! |
|
|
401 %! N = 100; |
|
|
402 %! A = zeros(N,N); |
|
|
403 %! for i=1:N-1 # form 1-D Laplacian matrix |
|
|
404 %! A(i:i+1,i:i+1) = [2 -1; -1 2]; |
|
|
405 %! endfor |
|
|
406 %! b = ones(N,1); X = A\b; #X is the true solution |
|
|
407 %! [x, flag, relres, iter, resvec] = pcr(A,b,1e-12); |
|
|
408 %! assert(flag,1); |
|
|
409 %! assert(relres>0.6); |
|
|
410 %! assert(iter,20); |
|
|
411 %! |
|
|
412 %!test |
|
|
413 %! |
|
|
414 %! #solve tridiagonal system with 'prefect' preconditioner |
|
|
415 %! #converges in one iteration |
|
|
416 %! |
|
|
417 %! N = 100; |
|
|
418 %! A = zeros(N,N); |
|
|
419 %! for i=1:N-1 # form 1-D Laplacian matrix |
|
|
420 %! A(i:i+1,i:i+1) = [2 -1; -1 2]; |
|
|
421 %! endfor |
|
|
422 %! b = ones(N,1); X = A\b; #X is the true solution |
|
|
423 %! [x, flag, relres, iter] = pcr(A,b,[],[],A,b); |
|
|
424 %! assert(norm(x-X)/norm(X)<1e-6); |
|
|
425 %! assert(relres<1e-6); |
|
|
426 %! assert(flag,0); |
|
|
427 %! assert(iter,1); #should converge in one iteration |
|
|
428 %! |
