Mercurial > hg > octave-nkf

comparison scripts/sparse/pcr.m @ 10549:95c3e38098bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Untabify .m scripts
author Rik <code@nomad.inbox5.com>
date Fri, 23 Apr 2010 11:28:50 -0700
parents 1bf0ce0930be
children 3140cb7a05a1
comparison
equal deleted inserted replaced
10548:479536c5bb10 10549:95c3e38098bf
98 ## Let us consider a trivial problem with a diagonal matrix (we exploit the 98 ## Let us consider a trivial problem with a diagonal matrix (we exploit the
99 ## sparsity of A) 99 ## sparsity of A)
100 ## 100 ##
101 ## @example 101 ## @example
102 ## @group 102 ## @group
103 ## n = 10; 103 ## n = 10;
104 ## a = sparse (diag (1:n)); 104 ## a = sparse (diag (1:n));
105 ## b = rand (N, 1); 105 ## b = rand (N, 1);
106 ## @end group 106 ## @end group
107 ## @end example 107 ## @end example
108 ## 108 ##
109 ## @sc{Example 1:} Simplest use of @code{pcr} 109 ## @sc{Example 1:} Simplest use of @code{pcr}
110 ## 110 ##
129 ## preconditioner (quite strange, because even the original matrix 129 ## preconditioner (quite strange, because even the original matrix
130 ## @var{a} is trivial) is defined as a function 130 ## @var{a} is trivial) is defined as a function
131 ## 131 ##
132 ## @example 132 ## @example
133 ## @group 133 ## @group
134 ## function y = apply_m (x) 134 ## function y = apply_m (x)
135 ## k = floor (length(x)-2); 135 ## k = floor (length(x)-2);
136 ## y = x; 136 ## y = x;
137 ## y(1:k) = x(1:k)./[1:k]'; 137 ## y(1:k) = x(1:k)./[1:k]';
138 ## endfunction 138 ## endfunction
139 ## 139 ##
140 ## [x, flag, relres, iter, resvec] = ... 140 ## [x, flag, relres, iter, resvec] = ...
141 ## pcr (a, b, [], [], "apply_m") 141 ## pcr (a, b, [], [], "apply_m")
142 ## semilogy([1:iter+1], resvec); 142 ## semilogy([1:iter+1], resvec);
148 ## 148 ##
149 ## @example 149 ## @example
150 ## @group 150 ## @group
151 ## function y = apply_m (x, varargin) 151 ## function y = apply_m (x, varargin)
152 ## k = varargin@{1@}; 152 ## k = varargin@{1@};
153 ## y = x; y(1:k) = x(1:k)./[1:k]'; 153 ## y = x; y(1:k) = x(1:k)./[1:k]';
154 ## endfunction 154 ## endfunction
155 ## 155 ##
156 ## [x, flag, relres, iter, resvec] = ... 156 ## [x, flag, relres, iter, resvec] = ...
157 ## pcr (a, b, [], [], "apply_m"', [], 3) 157 ## pcr (a, b, [], [], "apply_m"', [], 3)
158 ## @end group 158 ## @end group
159 ## @end example 159 ## @end example
160 ## 160 ##
161 ## @sc{References} 161 ## @sc{References}
162 ## 162 ##
163 ## [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of 163 ## [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of
164 ## Equations", section 9.5.4; Springer, 1994 164 ## Equations", section 9.5.4; Springer, 1994
165 ## 165 ##
166 ## @seealso{sparse, pcg} 166 ## @seealso{sparse, pcg}
167 ## @end deftypefn 167 ## @end deftypefn
168 168
169 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> 169 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
195 if (nargin < 2) 195 if (nargin < 2)
196 print_usage (); 196 print_usage ();
197 endif 197 endif
198 198
199 ## init 199 ## init
200 if (isnumeric (a)) # is A a matrix? 200 if (isnumeric (a)) # is A a matrix?
201 r = b - a*x; 201 r = b - a*x;
202 else # then A should be a function! 202 else # then A should be a function!
203 r = b - feval (a, x, varargin{:}); 203 r = b - feval (a, x, varargin{:});
204 endif 204 endif
205 205
206 if (isnumeric (m)) # is M a matrix? 206 if (isnumeric (m)) # is M a matrix?
207 if (isempty (m)) # if M is empty, use no precond 207 if (isempty (m)) # if M is empty, use no precond
208 p = r; 208 p = r;
209 else # otherwise, apply the precond 209 else # otherwise, apply the precond
210 p = m \ r; 210 p = m \ r;
211 endif 211 endif
212 else # then M should be a function! 212 else # then M should be a function!
213 p = feval (m, r, varargin{:}); 213 p = feval (m, r, varargin{:});
214 endif 214 endif
215 215
216 iter = 2; 216 iter = 2;
217 217
218 b_bot_old = 1; 218 b_bot_old = 1;
219 q_old = p_old = s_old = zeros (size (x)); 219 q_old = p_old = s_old = zeros (size (x));
220 220
221 if (isnumeric (a)) # is A a matrix? 221 if (isnumeric (a)) # is A a matrix?
222 q = a * p; 222 q = a * p;
223 else # then A should be a function! 223 else # then A should be a function!
224 q = feval (a, p, varargin{:}); 224 q = feval (a, p, varargin{:});
225 endif 225 endif
226 226
227 resvec(1) = abs (norm (r)); 227 resvec(1) = abs (norm (r));
228 228
229 ## iteration 229 ## iteration
230 while (resvec(iter-1) > tol*resvec(1) && iter < maxit) 230 while (resvec(iter-1) > tol*resvec(1) && iter < maxit)
231 231
232 if (isnumeric (m)) # is M a matrix? 232 if (isnumeric (m)) # is M a matrix?
233 if (isempty (m)) # if M is empty, use no precond 233 if (isempty (m)) # if M is empty, use no precond
234 s = q; 234 s = q;
235 else # otherwise, apply the precond 235 else # otherwise, apply the precond
236 s = m \ q; 236 s = m \ q;
237 endif 237 endif
238 else # then M should be a function! 238 else # then M should be a function!
239 s = feval (m, q, varargin{:}); 239 s = feval (m, q, varargin{:});
240 endif 240 endif
241 b_top = r' * s; 241 b_top = r' * s;
242 b_bot = q' * s; 242 b_bot = q' * s;
243 243
244 if (b_bot == 0.0) 244 if (b_bot == 0.0)
245 breakdown = true; 245 breakdown = true;
246 break; 246 break;
247 endif 247 endif
248 lambda = b_top / b_bot; 248 lambda = b_top / b_bot;
249 249
250 x += lambda*p; 250 x += lambda*p;
251 r -= lambda*q; 251 r -= lambda*q;
252 252
253 if (isnumeric(a)) # is A a matrix? 253 if (isnumeric(a)) # is A a matrix?
254 t = a*s; 254 t = a*s;
255 else # then A should be a function! 255 else # then A should be a function!
256 t = feval (a, s, varargin{:}); 256 t = feval (a, s, varargin{:});
257 endif 257 endif
258 258
259 alpha0 = (t'*s) / b_bot; 259 alpha0 = (t'*s) / b_bot;
260 alpha1 = (t'*s_old) / b_bot_old; 260 alpha1 = (t'*s_old) / b_bot_old;
261 261
262 p_temp = p; 262 p_temp = p;
263 q_temp = q; 263 q_temp = q;
264 264
265 p = s - alpha0*p - alpha1*p_old; 265 p = s - alpha0*p - alpha1*p_old;
266 q = t - alpha0*q - alpha1*q_old; 266 q = t - alpha0*q - alpha1*q_old;
267 267
268 s_old = s; 268 s_old = s;
269 p_old = p_temp; 269 p_old = p_temp;
270 q_old = q_temp; 270 q_old = q_temp;
271 b_bot_old = b_bot; 271 b_bot_old = b_bot;
272 272
273 resvec(iter) = abs (norm (r)); 273 resvec(iter) = abs (norm (r));
274 iter++; 274 iter++;
275 endwhile 275 endwhile
276 276
277 flag = 0; 277 flag = 0;
284 warning ("the initial residual norm was reduced %g times.\n", 1.0/relres); 284 warning ("the initial residual norm was reduced %g times.\n", 1.0/relres);
285 endif 285 endif
286 elseif (nargout < 2 && ! breakdown) 286 elseif (nargout < 2 && ! breakdown)
287 fprintf (stderr, "pcr: converged in %d iterations. \n", iter); 287 fprintf (stderr, "pcr: converged in %d iterations. \n", iter);
288 fprintf (stderr, "the initial residual norm was reduced %g times.\n", 288 fprintf (stderr, "the initial residual norm was reduced %g times.\n",
289 1.0 / relres); 289 1.0 / relres);
290 endif 290 endif
291 291
292 if (breakdown) 292 if (breakdown)
293 flag = 3; 293 flag = 3;
294 if (nargout < 2) 294 if (nargout < 2)
299 299
300 endfunction 300 endfunction
301 301
302 %!demo 302 %!demo
303 %! 303 %!
304 %! # Simplest usage of PCR (see also 'help pcr') 304 %! # Simplest usage of PCR (see also 'help pcr')
305 %! 305 %!
306 %! N = 20; 306 %! N = 20;
307 %! A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution 307 %! A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
308 %! x = pcr(A,b); 308 %! x = pcr(A,b);
309 %! printf('The solution relative error is %g\n', norm(x-y)/norm(y)); 309 %! printf('The solution relative error is %g\n', norm(x-y)/norm(y));
310 %! 310 %!
311 %! # You shouldn't be afraid if PCR issues some warning messages in this 311 %! # You shouldn't be afraid if PCR issues some warning messages in this
312 %! # example: watch out in the second example, why it takes N iterations 312 %! # example: watch out in the second example, why it takes N iterations
313 %! # of PCR to converge to (a very accurate, by the way) solution 313 %! # of PCR to converge to (a very accurate, by the way) solution
314 %!demo 314 %!demo
315 %! 315 %!
316 %! # Full output from PCR 316 %! # Full output from PCR
317 %! # We use this output to plot the convergence history 317 %! # We use this output to plot the convergence history
318 %! 318 %!
319 %! N = 20; 319 %! N = 20;
320 %! A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution 320 %! A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
321 %! [x, flag, relres, iter, resvec] = pcr(A,b); 321 %! [x, flag, relres, iter, resvec] = pcr(A,b);
322 %! printf('The solution relative error is %g\n', norm(x-X)/norm(X)); 322 %! printf('The solution relative error is %g\n', norm(x-X)/norm(X));
323 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); 323 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
324 %! semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;'); 324 %! semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');
325 %!demo 325 %!demo
326 %! 326 %!
327 %! # Full output from PCR 327 %! # Full output from PCR
328 %! # We use indefinite matrix based on the Hilbert matrix, with one 328 %! # We use indefinite matrix based on the Hilbert matrix, with one
329 %! # strongly negative eigenvalue 329 %! # strongly negative eigenvalue
330 %! # Hilbert matrix is extremely ill conditioned, so is ours, 330 %! # Hilbert matrix is extremely ill conditioned, so is ours,
331 %! # and that's why PCR WILL have problems 331 %! # and that's why PCR WILL have problems
332 %! 332 %!
333 %! N = 10; 333 %! N = 10;
334 %! A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution 334 %! A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
335 %! printf('Condition number of A is %g\n', cond(A)); 335 %! printf('Condition number of A is %g\n', cond(A));
336 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],200); 336 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],200);
337 %! if (flag == 3) 337 %! if (flag == 3)
338 %! printf('PCR breakdown. System matrix is [close to] singular\n'); 338 %! printf('PCR breakdown. System matrix is [close to] singular\n');
339 %! end 339 %! end
340 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); 340 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
341 %! semilogy([0:iter],resvec,'o-g;absolute residual;'); 341 %! semilogy([0:iter],resvec,'o-g;absolute residual;');
342 %!demo 342 %!demo
343 %! 343 %!
344 %! # Full output from PCR 344 %! # Full output from PCR
345 %! # We use an indefinite matrix based on the 1-D Laplacian matrix for A, 345 %! # We use an indefinite matrix based on the 1-D Laplacian matrix for A,
346 %! # and here we have cond(A) = O(N^2) 346 %! # and here we have cond(A) = O(N^2)
347 %! # That's the reason we need some preconditioner; here we take 347 %! # That's the reason we need some preconditioner; here we take
348 %! # a very simple and not powerful Jacobi preconditioner, 348 %! # a very simple and not powerful Jacobi preconditioner,
349 %! # which is the diagonal of A 349 %! # which is the diagonal of A
350 %! 350 %!
351 %! # Note that we use here indefinite preconditioners! 351 %! # Note that we use here indefinite preconditioners!
352 %! 352 %!
353 %! N = 100; 353 %! N = 100;
354 %! A = zeros(N,N); 354 %! A = zeros(N,N);
355 %! for i=1:N-1 # form 1-D Laplacian matrix 355 %! for i=1:N-1 # form 1-D Laplacian matrix
356 %! A(i:i+1,i:i+1) = [2 -1; -1 2]; 356 %! A(i:i+1,i:i+1) = [2 -1; -1 2];
357 %! endfor 357 %! endfor
358 %! A = [A, zeros(size(A)); zeros(size(A)), -A]; 358 %! A = [A, zeros(size(A)); zeros(size(A)), -A];
359 %! b = rand(2*N,1); X = A\b; #X is the true solution 359 %! b = rand(2*N,1); X = A\b; #X is the true solution
360 %! maxit = 80; 360 %! maxit = 80;
361 %! printf('System condition number is %g\n',cond(A)); 361 %! printf('System condition number is %g\n',cond(A));
362 %! # No preconditioner: the convergence is very slow! 362 %! # No preconditioner: the convergence is very slow!
363 %! 363 %!
364 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit); 364 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit);
365 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); 365 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
366 %! semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;'); 366 %! semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;');
367 %! 367 %!
368 %! pause(1); 368 %! pause(1);
369 %! # Test Jacobi preconditioner: it will not help much!!! 369 %! # Test Jacobi preconditioner: it will not help much!!!
370 %! 370 %!
371 %! M = diag(diag(A)); # Jacobi preconditioner 371 %! M = diag(diag(A)); # Jacobi preconditioner
372 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M); 372 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
373 %! hold on; 373 %! hold on;
374 %! semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;'); 374 %! semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');
375 %! 375 %!
376 %! pause(1); 376 %! pause(1);
377 %! # Test nonoverlapping block Jacobi preconditioner: this one should give 377 %! # Test nonoverlapping block Jacobi preconditioner: this one should give
378 %! # some convergence speedup! 378 %! # some convergence speedup!
379 %! 379 %!
380 %! M = zeros(N,N);k=4; 380 %! M = zeros(N,N);k=4;
381 %! for i=1:k:N # get k x k diagonal blocks of A 381 %! for i=1:k:N # get k x k diagonal blocks of A
382 %! M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1); 382 %! M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
383 %! endfor 383 %! endfor
384 %! M = [M, zeros(size(M)); zeros(size(M)), -M]; 384 %! M = [M, zeros(size(M)); zeros(size(M)), -M];
385 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M); 385 %! [x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
386 %! semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;'); 386 %! semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
387 %! hold off; 387 %! hold off;
388 %!test 388 %!test
389 %! 389 %!
390 %! #solve small indefinite diagonal system 390 %! #solve small indefinite diagonal system
391 %! 391 %!
392 %! N = 10; 392 %! N = 10;
393 %! A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution 393 %! A = diag(linspace(-10.1,10,N)); b = ones(N,1); X = A\b; #X is the true solution
394 %! [x, flag] = pcr(A,b,[],N+1); 394 %! [x, flag] = pcr(A,b,[],N+1);
395 %! assert(norm(x-X)/norm(X)<1e-10); 395 %! assert(norm(x-X)/norm(X)<1e-10);
396 %! assert(flag,0); 396 %! assert(flag,0);
397 %! 397 %!
398 %!test 398 %!test
399 %! 399 %!
400 %! #solve tridiagonal system, do not converge in default 20 iterations 400 %! #solve tridiagonal system, do not converge in default 20 iterations
401 %! #should perform max allowable default number of iterations 401 %! #should perform max allowable default number of iterations
402 %! 402 %!
403 %! N = 100; 403 %! N = 100;
404 %! A = zeros(N,N); 404 %! A = zeros(N,N);
405 %! for i=1:N-1 # form 1-D Laplacian matrix 405 %! for i=1:N-1 # form 1-D Laplacian matrix
406 %! A(i:i+1,i:i+1) = [2 -1; -1 2]; 406 %! A(i:i+1,i:i+1) = [2 -1; -1 2];
407 %! endfor 407 %! endfor
408 %! b = ones(N,1); X = A\b; #X is the true solution 408 %! b = ones(N,1); X = A\b; #X is the true solution
409 %! [x, flag, relres, iter, resvec] = pcr(A,b,1e-12); 409 %! [x, flag, relres, iter, resvec] = pcr(A,b,1e-12);
410 %! assert(flag,1); 410 %! assert(flag,1);
411 %! assert(relres>0.6); 411 %! assert(relres>0.6);
412 %! assert(iter,20); 412 %! assert(iter,20);
413 %! 413 %!
414 %!test 414 %!test
415 %! 415 %!
416 %! #solve tridiagonal system with 'prefect' preconditioner 416 %! #solve tridiagonal system with 'prefect' preconditioner
417 %! #converges in one iteration 417 %! #converges in one iteration
418 %! 418 %!
419 %! N = 100; 419 %! N = 100;
420 %! A = zeros(N,N); 420 %! A = zeros(N,N);
421 %! for i=1:N-1 # form 1-D Laplacian matrix 421 %! for i=1:N-1 # form 1-D Laplacian matrix
422 %! A(i:i+1,i:i+1) = [2 -1; -1 2]; 422 %! A(i:i+1,i:i+1) = [2 -1; -1 2];
423 %! endfor 423 %! endfor
424 %! b = ones(N,1); X = A\b; #X is the true solution 424 %! b = ones(N,1); X = A\b; #X is the true solution
425 %! [x, flag, relres, iter] = pcr(A,b,[],[],A,b); 425 %! [x, flag, relres, iter] = pcr(A,b,[],[],A,b);
426 %! assert(norm(x-X)/norm(X)<1e-6); 426 %! assert(norm(x-X)/norm(X)<1e-6);
427 %! assert(relres<1e-6); 427 %! assert(relres<1e-6);
428 %! assert(flag,0); 428 %! assert(flag,0);
429 %! assert(iter,1); #should converge in one iteration 429 %! assert(iter,1); #should converge in one iteration
430 %! 430 %!