Mercurial > octave-nkf

comparison scripts/sparse/pcg.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
120 ## Let us consider a trivial problem with a diagonal matrix (we exploit the 120 ## Let us consider a trivial problem with a diagonal matrix (we exploit the
121 ## sparsity of A) 121 ## sparsity of A)
122 ## 122 ##
123 ## @example 123 ## @example
124 ## @group 124 ## @group
125 ## n = 10; 125 ## n = 10;
126 ## a = diag (sparse (1:n)); 126 ## a = diag (sparse (1:n));
127 ## b = rand (n, 1); 127 ## b = rand (n, 1);
128 ## [l, u, p, q] = luinc (a, 1.e-3); 128 ## [l, u, p, q] = luinc (a, 1.e-3);
129 ## @end group 129 ## @end group
130 ## @end example 130 ## @end example
131 ## 131 ##
132 ## @sc{Example 1:} Simplest use of @code{pcg} 132 ## @sc{Example 1:} Simplest use of @code{pcg}
196 ## @end group 196 ## @end group
197 ## @end example 197 ## @end example
198 ## 198 ##
199 ## @sc{References} 199 ## @sc{References}
200 ## 200 ##
201 ## [1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations', 201 ## [1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations',
202 ## SIAM, 1995 (the base PCG algorithm) 202 ## SIAM, 1995 (the base PCG algorithm)
203 ## 203 ##
204 ## [2] Y.Saad, 'Iterative methods for sparse linear systems', PWS 1996 204 ## [2] Y.Saad, 'Iterative methods for sparse linear systems', PWS 1996
205 ## (condition number estimate from PCG) Revised version of this book is 205 ## (condition number estimate from PCG) Revised version of this book is
206 ## available online at http://www-users.cs.umn.edu/~saad/books.html 206 ## available online at http://www-users.cs.umn.edu/~saad/books.html
207 ## 207 ##
208 ## 208 ##
209 ## @seealso{sparse, pcr} 209 ## @seealso{sparse, pcr}
210 ## @end deftypefn 210 ## @end deftypefn
211 211
270 iter = 2; 270 iter = 2;
271 271
272 while (resvec (iter-1,1) > tol * resvec (1,1) && iter < maxit) 272 while (resvec (iter-1,1) > tol * resvec (1,1) && iter < maxit)
273 if (exist_m1) 273 if (exist_m1)
274 if(isnumeric (m1)) 274 if(isnumeric (m1))
275 y = m1 \ r; 275 y = m1 \ r;
276 else 276 else
277 y = feval (m1, r, varargin{:}); 277 y = feval (m1, r, varargin{:});
278 endif 278 endif
279 else 279 else
280 y = r; 280 y = r;
281 endif 281 endif
282 if (exist_m2) 282 if (exist_m2)
283 if (isnumeric (m2)) 283 if (isnumeric (m2))
284 z = m2 \ y; 284 z = m2 \ y;
285 else 285 else
286 z = feval (m2, y, varargin{:}); 286 z = feval (m2, y, varargin{:});
287 endif 287 endif
288 else 288 else
289 z = y; 289 z = y;
290 endif 290 endif
291 tau = z' * r; 291 tau = z' * r;
309 endif 309 endif
310 x += alpha * p; 310 x += alpha * p;
311 r -= alpha * w; 311 r -= alpha * w;
312 if (nargout > 5 && iter > 2) 312 if (nargout > 5 && iter > 2)
313 T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ... 313 T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
314 [1 sqrt(beta); sqrt(beta) beta]./oldalpha; 314 [1 sqrt(beta); sqrt(beta) beta]./oldalpha;
315 ## EVS = eig(T(2:iter-1,2:iter-1)); 315 ## EVS = eig(T(2:iter-1,2:iter-1));
316 ## fprintf(stderr,"PCG condest: %g (iteration: %d)\n", max(EVS)/min(EVS),iter); 316 ## fprintf(stderr,"PCG condest: %g (iteration: %d)\n", max(EVS)/min(EVS),iter);
317 endif 317 endif
318 resvec (iter,1) = norm (r); 318 resvec (iter,1) = norm (r);
319 iter++; 319 iter++;
320 endwhile 320 endwhile
321 321
322 if (nargout > 5) 322 if (nargout > 5)
323 if (matrix_positive_definite) 323 if (matrix_positive_definite)
324 if (iter > 3) 324 if (iter > 3)
325 T = T(2:iter-2,2:iter-2); 325 T = T(2:iter-2,2:iter-2);
326 l = eig (T); 326 l = eig (T);
327 eigest = [min(l), max(l)]; 327 eigest = [min(l), max(l)];
328 ## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1)); 328 ## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1));
329 else 329 else
330 eigest = [NaN, NaN]; 330 eigest = [NaN, NaN];
331 warning ("pcg: eigenvalue estimate failed: iteration converged too fast."); 331 warning ("pcg: eigenvalue estimate failed: iteration converged too fast.");
332 endif 332 endif
333 else 333 else
334 eigest = [NaN, NaN]; 334 eigest = [NaN, NaN];
335 endif 335 endif
336 336
337 ## Apply the preconditioner once more and finish with the precond 337 ## Apply the preconditioner once more and finish with the precond
338 ## residual. 338 ## residual.
339 if (exist_m1) 339 if (exist_m1)
340 if (isnumeric (m1)) 340 if (isnumeric (m1))
341 y = m1 \ r; 341 y = m1 \ r;
342 else 342 else
343 y = feval (m1, r, varargin{:}); 343 y = feval (m1, r, varargin{:});
344 endif 344 endif
345 else 345 else
346 y = r; 346 y = r;
347 endif 347 endif
348 if (exist_m2) 348 if (exist_m2)
349 if (isnumeric (m2)) 349 if (isnumeric (m2))
350 z = m2 \ y; 350 z = m2 \ y;
351 else 351 else
352 z = feval (m2, y, varargin{:}); 352 z = feval (m2, y, varargin{:});
353 endif 353 endif
354 else 354 else
355 z = y; 355 z = y;
356 endif 356 endif
357 357
366 if (iter >= maxit - 2) 366 if (iter >= maxit - 2)
367 flag = 1; 367 flag = 1;
368 if (nargout < 2) 368 if (nargout < 2)
369 warning ("pcg: maximum number of iterations (%d) reached\n", iter); 369 warning ("pcg: maximum number of iterations (%d) reached\n", iter);
370 warning ("the initial residual norm was reduced %g times.\n", ... 370 warning ("the initial residual norm was reduced %g times.\n", ...
371 1.0 / relres); 371 1.0 / relres);
372 endif 372 endif
373 elseif (nargout < 2) 373 elseif (nargout < 2)
374 fprintf (stderr, "pcg: converged in %d iterations. ", iter); 374 fprintf (stderr, "pcg: converged in %d iterations. ", iter);
375 fprintf (stderr, "the initial residual norm was reduced %g times.\n",... 375 fprintf (stderr, "the initial residual norm was reduced %g times.\n",...
376 1.0/relres); 376 1.0/relres);
377 endif 377 endif
378 378
379 if (! matrix_positive_definite) 379 if (! matrix_positive_definite)
380 flag = 3; 380 flag = 3;
381 if (nargout < 2) 381 if (nargout < 2)
384 endif 384 endif
385 endfunction 385 endfunction
386 386
387 %!demo 387 %!demo
388 %! 388 %!
389 %! # Simplest usage of pcg (see also 'help pcg') 389 %! # Simplest usage of pcg (see also 'help pcg')
390 %! 390 %!
391 %! N = 10; 391 %! N = 10;
392 %! A = diag ([1:N]); b = rand (N, 1); y = A \ b; #y is the true solution 392 %! A = diag ([1:N]); b = rand (N, 1); y = A \ b; #y is the true solution
393 %! x = pcg (A, b); 393 %! x = pcg (A, b);
394 %! printf('The solution relative error is %g\n', norm (x - y) / norm (y)); 394 %! printf('The solution relative error is %g\n', norm (x - y) / norm (y));
395 %! 395 %!
396 %! # You shouldn't be afraid if pcg issues some warning messages in this 396 %! # You shouldn't be afraid if pcg issues some warning messages in this
397 %! # example: watch out in the second example, why it takes N iterations 397 %! # example: watch out in the second example, why it takes N iterations
398 %! # of pcg to converge to (a very accurate, by the way) solution 398 %! # of pcg to converge to (a very accurate, by the way) solution
399 %!demo 399 %!demo
400 %! 400 %!
401 %! # Full output from pcg, except for the eigenvalue estimates 401 %! # Full output from pcg, except for the eigenvalue estimates
402 %! # We use this output to plot the convergence history 402 %! # We use this output to plot the convergence history
403 %! 403 %!
404 %! N = 10; 404 %! N = 10;
405 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution 405 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution
406 %! [x, flag, relres, iter, resvec] = pcg (A, b); 406 %! [x, flag, relres, iter, resvec] = pcg (A, b);
407 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X)); 407 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X));
408 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); 408 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
409 %! semilogy([0:iter], resvec / resvec(1),'o-g'); 409 %! semilogy([0:iter], resvec / resvec(1),'o-g');
410 %! legend('relative residual'); 410 %! legend('relative residual');
411 %!demo 411 %!demo
412 %! 412 %!
413 %! # Full output from pcg, including the eigenvalue estimates 413 %! # Full output from pcg, including the eigenvalue estimates
414 %! # Hilbert matrix is extremely ill conditioned, so pcg WILL have problems 414 %! # Hilbert matrix is extremely ill conditioned, so pcg WILL have problems
415 %! 415 %!
416 %! N = 10; 416 %! N = 10;
417 %! A = hilb (N); b = rand (N, 1); X = A \ b; #X is the true solution 417 %! A = hilb (N); b = rand (N, 1); X = A \ b; #X is the true solution
418 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200); 418 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200);
419 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X)); 419 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X));
420 %! printf('Condition number estimate is %g\n', eigest(2) / eigest (1)); 420 %! printf('Condition number estimate is %g\n', eigest(2) / eigest (1));
421 %! printf('Actual condition number is %g\n', cond (A)); 421 %! printf('Actual condition number is %g\n', cond (A));
422 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); 422 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
423 %! semilogy([0:iter], resvec,['o-g';'+-r']); 423 %! semilogy([0:iter], resvec,['o-g';'+-r']);
424 %! legend('absolute residual','absolute preconditioned residual'); 424 %! legend('absolute residual','absolute preconditioned residual');
425 %!demo 425 %!demo
426 %! 426 %!
427 %! # Full output from pcg, including the eigenvalue estimates 427 %! # Full output from pcg, including the eigenvalue estimates
428 %! # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2) 428 %! # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2)
429 %! # and that's the reasone we need some preconditioner; here we take 429 %! # and that's the reasone we need some preconditioner; here we take
430 %! # a very simple and not powerful Jacobi preconditioner, 430 %! # a very simple and not powerful Jacobi preconditioner,
431 %! # which is the diagonal of A 431 %! # which is the diagonal of A
432 %! 432 %!
433 %! N = 100; 433 %! N = 100;
434 %! A = zeros (N, N); 434 %! A = zeros (N, N);
435 %! for i=1 : N - 1 # form 1-D Laplacian matrix 435 %! for i=1 : N - 1 # form 1-D Laplacian matrix
436 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; 436 %! A (i:i+1, i:i+1) = [2 -1; -1 2];
437 %! endfor 437 %! endfor
438 %! b = rand (N, 1); X = A \ b; #X is the true solution 438 %! b = rand (N, 1); X = A \ b; #X is the true solution
439 %! maxit = 80; 439 %! maxit = 80;
440 %! printf('System condition number is %g\n', cond (A)); 440 %! printf('System condition number is %g\n', cond (A));
441 %! # No preconditioner: the convergence is very slow! 441 %! # No preconditioner: the convergence is very slow!
442 %! 442 %!
443 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit); 443 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit);
444 %! printf('System condition number estimate is %g\n', eigest(2) / eigest(1)); 444 %! printf('System condition number estimate is %g\n', eigest(2) / eigest(1));
445 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); 445 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
446 %! semilogy([0:iter], resvec(:,1), 'o-g'); 446 %! semilogy([0:iter], resvec(:,1), 'o-g');
447 %! legend('NO preconditioning: absolute residual'); 447 %! legend('NO preconditioning: absolute residual');
448 %! 448 %!
449 %! pause(1); 449 %! pause(1);
450 %! # Test Jacobi preconditioner: it will not help much!!! 450 %! # Test Jacobi preconditioner: it will not help much!!!
451 %! 451 %!
452 %! M = diag (diag (A)); # Jacobi preconditioner 452 %! M = diag (diag (A)); # Jacobi preconditioner
453 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); 453 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
454 %! printf('JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1)); 454 %! printf('JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1));
455 %! hold on; 455 %! hold on;
456 %! semilogy([0:iter], resvec(:,1), 'o-r'); 456 %! semilogy([0:iter], resvec(:,1), 'o-r');
457 %! legend('NO preconditioning: absolute residual', ... 457 %! legend('NO preconditioning: absolute residual', ...
458 %! 'JACOBI preconditioner: absolute residual'); 458 %! 'JACOBI preconditioner: absolute residual');
459 %! 459 %!
460 %! pause(1); 460 %! pause(1);
461 %! # Test nonoverlapping block Jacobi preconditioner: it will help much! 461 %! # Test nonoverlapping block Jacobi preconditioner: it will help much!
462 %! 462 %!
463 %! M = zeros (N, N); k = 4; 463 %! M = zeros (N, N); k = 4;
464 %! for i = 1 : k : N # form 1-D Laplacian matrix 464 %! for i = 1 : k : N # form 1-D Laplacian matrix
465 %! M (i:i+k-1, i:i+k-1) = A (i:i+k-1, i:i+k-1); 465 %! M (i:i+k-1, i:i+k-1) = A (i:i+k-1, i:i+k-1);
466 %! endfor 466 %! endfor
467 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); 467 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
468 %! printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1)); 468 %! printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1));
469 %! semilogy ([0:iter], resvec(:,1),'o-b'); 469 %! semilogy ([0:iter], resvec(:,1),'o-b');
470 %! legend('NO preconditioning: absolute residual', ... 470 %! legend('NO preconditioning: absolute residual', ...
471 %! 'JACOBI preconditioner: absolute residual', ... 471 %! 'JACOBI preconditioner: absolute residual', ...
472 %! 'BLOCK JACOBI preconditioner: absolute residual'); 472 %! 'BLOCK JACOBI preconditioner: absolute residual');
473 %! hold off; 473 %! hold off;
474 %!test 474 %!test
475 %! 475 %!
476 %! #solve small diagonal system 476 %! #solve small diagonal system
477 %! 477 %!
478 %! N = 10; 478 %! N = 10;
479 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution 479 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution
480 %! [x, flag] = pcg (A, b, [], N+1); 480 %! [x, flag] = pcg (A, b, [], N+1);
481 %! assert(norm (x - X) / norm (X), 0, 1e-10); 481 %! assert(norm (x - X) / norm (X), 0, 1e-10);
482 %! assert(flag, 0); 482 %! assert(flag, 0);
483 %! 483 %!
484 %!test 484 %!test
485 %! 485 %!
486 %! #solve small indefinite diagonal system 486 %! #solve small indefinite diagonal system
487 %! #despite A is indefinite, the iteration continues and converges 487 %! #despite A is indefinite, the iteration continues and converges
488 %! #indefiniteness of A is detected 488 %! #indefiniteness of A is detected
489 %! 489 %!
490 %! N = 10; 490 %! N = 10;
491 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1); X = A \ b; #X is the true solution 491 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1); X = A \ b; #X is the true solution
492 %! [x, flag] = pcg (A, b, [], N+1); 492 %! [x, flag] = pcg (A, b, [], N+1);
493 %! assert(norm (x - X) / norm (X), 0, 1e-10); 493 %! assert(norm (x - X) / norm (X), 0, 1e-10);
494 %! assert(flag, 3); 494 %! assert(flag, 3);
495 %! 495 %!
496 %!test 496 %!test
497 %! 497 %!
498 %! #solve tridiagonal system, do not converge in default 20 iterations 498 %! #solve tridiagonal system, do not converge in default 20 iterations
499 %! 499 %!
500 %! N = 100; 500 %! N = 100;
501 %! A = zeros (N, N); 501 %! A = zeros (N, N);
502 %! for i = 1 : N - 1 # form 1-D Laplacian matrix 502 %! for i = 1 : N - 1 # form 1-D Laplacian matrix
503 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; 503 %! A (i:i+1, i:i+1) = [2 -1; -1 2];
504 %! endfor 504 %! endfor
505 %! b = ones (N, 1); X = A \ b; #X is the true solution 505 %! b = ones (N, 1); X = A \ b; #X is the true solution
506 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12); 506 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12);
507 %! assert(flag); 507 %! assert(flag);
508 %! assert(relres > 1.0); 508 %! assert(relres > 1.0);
509 %! assert(iter, 20); #should perform max allowable default number of iterations 509 %! assert(iter, 20); #should perform max allowable default number of iterations
510 %! 510 %!
511 %!test 511 %!test
512 %! 512 %!
513 %! #solve tridiagonal system with 'prefect' preconditioner 513 %! #solve tridiagonal system with 'prefect' preconditioner
514 %! #converges in one iteration, so the eigest does not work 514 %! #converges in one iteration, so the eigest does not work
515 %! #and issues a warning 515 %! #and issues a warning
516 %! 516 %!
517 %! N = 100; 517 %! N = 100;
518 %! A = zeros (N, N); 518 %! A = zeros (N, N);
519 %! for i = 1 : N - 1 # form 1-D Laplacian matrix 519 %! for i = 1 : N - 1 # form 1-D Laplacian matrix
520 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; 520 %! A (i:i+1, i:i+1) = [2 -1; -1 2];
521 %! endfor 521 %! endfor
522 %! b = ones (N, 1); X = A \ b; #X is the true solution 522 %! b = ones (N, 1); X = A \ b; #X is the true solution
523 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b); 523 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
524 %! assert(norm (x - X) / norm (X), 0, 1e-6); 524 %! assert(norm (x - X) / norm (X), 0, 1e-6);
525 %! assert(flag, 0); 525 %! assert(flag, 0);
526 %! assert(iter, 1); #should converge in one iteration 526 %! assert(iter, 1); #should converge in one iteration
527 %! assert(isnan (eigest), isnan ([NaN, NaN])); 527 %! assert(isnan (eigest), isnan ([NaN, NaN]));
528 %! 528 %!