Mercurial > octave-nkf

comparison scripts/sparse/pcg.m @ 14237:11949c9795a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Add clf() to all demos using plot features to get reproducibility. Use 64 as input to all colormaps (jet (64)) to get reproducibility. * bicubic.m, cell2mat.m, celldisp.m, cplxpair.m, interp1.m, interp2.m, interpft.m, interpn.m, profile.m, profshow.m, convhull.m, delaunay.m, griddata.m, inpolygon.m, voronoi.m, autumn.m, bone.m, contrast.m, cool.m, copper.m, flag.m, gmap40.m, gray.m, hot.m, hsv.m, image.m, imshow.m, jet.m, ocean.m, pink.m, prism.m, rainbow.m, spring.m, summer.m, white.m, winter.m, condest.m, onenormest.m, axis.m, clabel.m, colorbar.m, comet.m, comet3.m, compass.m, contour.m, contour3.m, contourf.m, cylinder.m, daspect.m, ellipsoid.m, errorbar.m, ezcontour.m, ezcontourf.m, ezmesh.m, ezmeshc.m, ezplot.m, ezplot3.m, ezpolar.m, ezsurf.m, ezsurfc.m, feather.m, fill.m, fplot.m, grid.m, hold.m, isosurface.m, legend.m, loglog.m, loglogerr.m, pareto.m, patch.m, pbaspect.m, pcolor.m, pie.m, pie3.m, plot3.m, plotmatrix.m, plotyy.m, polar.m, quiver.m, quiver3.m, rectangle.m, refreshdata.m, ribbon.m, rose.m, scatter.m, scatter3.m, semilogx.m, semilogxerr.m, semilogy.m, semilogyerr.m, shading.m, slice.m, sombrero.m, stairs.m, stem.m, stem3.m, subplot.m, surf.m, surfc.m, surfl.m, surfnorm.m, text.m, title.m, trimesh.m, triplot.m, trisurf.m, uigetdir.m, uigetfile.m, uimenu.m, uiputfile.m, waitbar.m, xlim.m, ylim.m, zlim.m, mkpp.m, pchip.m, polyaffine.m, spline.m, bicgstab.m, cgs.m, gplot.m, pcg.m, pcr.m, treeplot.m, strtok.m, demo.m, example.m, rundemos.m, speed.m, test.m, calendar.m, datestr.m, datetick.m, weekday.m: Revamp %!demos to use Octave coding conventions on spacing, etc.
author Rik <octave@nomad.inbox5.com>
date Fri, 20 Jan 2012 12:59:53 -0800
parents 72c96de7a403
children ce2b59a6d0e5
comparison
equal deleted inserted replaced
14236:35903f035390 14237:11949c9795a0
387 warning ("pcg: matrix not positive definite?\n"); 387 warning ("pcg: matrix not positive definite?\n");
388 endif 388 endif
389 endif 389 endif
390 endfunction 390 endfunction
391 391
392
392 %!demo 393 %!demo
393 %! 394 %! # Simplest usage of pcg (see also 'help pcg')
394 %! # Simplest usage of pcg (see also 'help pcg') 395 %!
395 %! 396 %! N = 10;
396 %! N = 10; 397 %! A = diag ([1:N]); b = rand (N, 1);
397 %! A = diag ([1:N]); b = rand (N, 1); y = A \ b; #y is the true solution 398 %! y = A \ b; # y is the true solution
398 %! x = pcg (A, b); 399 %! x = pcg (A, b);
399 %! printf('The solution relative error is %g\n', norm (x - y) / norm (y)); 400 %! printf ("The solution relative error is %g\n", norm (x - y) / norm (y));
400 %! 401 %!
401 %! # You shouldn't be afraid if pcg issues some warning messages in this 402 %! # You shouldn't be afraid if pcg issues some warning messages in this
402 %! # example: watch out in the second example, why it takes N iterations 403 %! # example: watch out in the second example, why it takes N iterations
403 %! # of pcg to converge to (a very accurate, by the way) solution 404 %! # of pcg to converge to (a very accurate, by the way) solution
405
404 %!demo 406 %!demo
405 %! 407 %! # Full output from pcg, except for the eigenvalue estimates
406 %! # Full output from pcg, except for the eigenvalue estimates 408 %! # We use this output to plot the convergence history
407 %! # We use this output to plot the convergence history 409 %!
408 %! 410 %! N = 10;
409 %! N = 10; 411 %! A = diag ([1:N]); b = rand (N, 1);
410 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution 412 %! X = A \ b; # X is the true solution
411 %! [x, flag, relres, iter, resvec] = pcg (A, b); 413 %! [x, flag, relres, iter, resvec] = pcg (A, b);
412 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X)); 414 %! printf ("The solution relative error is %g\n", norm (x - X) / norm (X));
413 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); 415 %! title ("Convergence history");
414 %! semilogy([0:iter], resvec / resvec(1),'o-g'); 416 %! semilogy ([0:iter], resvec / resvec(1), "o-g");
415 %! legend('relative residual'); 417 %! xlabel ("Iteration"); ylabel ("log(||b-Ax||/||b||)");
418 %! legend ("relative residual");
419
416 %!demo 420 %!demo
417 %! 421 %! # Full output from pcg, including the eigenvalue estimates
418 %! # Full output from pcg, including the eigenvalue estimates 422 %! # Hilbert matrix is extremely ill-conditioned, so pcg WILL have problems
419 %! # Hilbert matrix is extremely ill conditioned, so pcg WILL have problems 423 %!
420 %! 424 %! N = 10;
421 %! N = 10; 425 %! A = hilb (N); b = rand (N, 1);
422 %! A = hilb (N); b = rand (N, 1); X = A \ b; #X is the true solution 426 %! X = A \ b; # X is the true solution
423 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200); 427 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200);
424 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X)); 428 %! printf ("The solution relative error is %g\n", norm (x - X) / norm (X));
425 %! printf('Condition number estimate is %g\n', eigest(2) / eigest (1)); 429 %! printf ("Condition number estimate is %g\n", eigest(2) / eigest(1));
426 %! printf('Actual condition number is %g\n', cond (A)); 430 %! printf ("Actual condition number is %g\n", cond (A));
427 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); 431 %! title ("Convergence history");
428 %! semilogy([0:iter], resvec,['o-g';'+-r']); 432 %! semilogy ([0:iter], resvec, ["o-g";"+-r"]);
429 %! legend('absolute residual','absolute preconditioned residual'); 433 %! xlabel ("Iteration"); ylabel ("log(||b-Ax||)");
434 %! legend ("absolute residual", "absolute preconditioned residual");
435
430 %!demo 436 %!demo
431 %! 437 %!
432 %! # Full output from pcg, including the eigenvalue estimates 438 %! # Full output from pcg, including the eigenvalue estimates
433 %! # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2) 439 %! # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2)
434 %! # and that's the reasone we need some preconditioner; here we take 440 %! # and that's the reason we need some preconditioner; here we take
435 %! # a very simple and not powerful Jacobi preconditioner, 441 %! # a very simple and not powerful Jacobi preconditioner,
436 %! # which is the diagonal of A 442 %! # which is the diagonal of A
437 %! 443 %!
438 %! N = 100; 444 %! N = 100;
439 %! A = zeros (N, N); 445 %! A = zeros (N, N);
440 %! for i=1 : N - 1 # form 1-D Laplacian matrix 446 %! for i = 1 : N - 1 # form 1-D Laplacian matrix
441 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; 447 %! A(i:i+1, i:i+1) = [2 -1; -1 2];
442 %! endfor 448 %! endfor
443 %! b = rand (N, 1); X = A \ b; #X is the true solution 449 %! b = rand (N, 1);
444 %! maxit = 80; 450 %! X = A \ b; # X is the true solution
445 %! printf('System condition number is %g\n', cond (A)); 451 %! maxit = 80;
446 %! # No preconditioner: the convergence is very slow! 452 %! printf ("System condition number is %g\n", cond (A));
447 %! 453 %! # No preconditioner: the convergence is very slow!
448 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit); 454 %!
449 %! printf('System condition number estimate is %g\n', eigest(2) / eigest(1)); 455 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit);
450 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); 456 %! printf ("System condition number estimate is %g\n", eigest(2) / eigest(1));
451 %! semilogy([0:iter], resvec(:,1), 'o-g'); 457 %! title ("Convergence history");
452 %! legend('NO preconditioning: absolute residual'); 458 %! semilogy ([0:iter], resvec(:,1), "o-g");
453 %! 459 %! xlabel ("Iteration"); ylabel ("log(||b-Ax||)");
454 %! pause(1); 460 %! legend ("NO preconditioning: absolute residual");
455 %! # Test Jacobi preconditioner: it will not help much!!! 461 %!
456 %! 462 %! pause (1);
457 %! M = diag (diag (A)); # Jacobi preconditioner 463 %! # Test Jacobi preconditioner: it will not help much!!!
458 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); 464 %!
459 %! printf('JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1)); 465 %! M = diag (diag (A)); # Jacobi preconditioner
460 %! hold on; 466 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
461 %! semilogy([0:iter], resvec(:,1), 'o-r'); 467 %! printf ("JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1));
462 %! legend('NO preconditioning: absolute residual', ... 468 %! hold on;
463 %! 'JACOBI preconditioner: absolute residual'); 469 %! semilogy ([0:iter], resvec(:,1), "o-r");
464 %! 470 %! legend ("NO preconditioning: absolute residual", ...
465 %! pause(1); 471 %! "JACOBI preconditioner: absolute residual");
466 %! # Test nonoverlapping block Jacobi preconditioner: it will help much! 472 %!
467 %! 473 %! pause (1);
468 %! M = zeros (N, N); k = 4; 474 %! # Test nonoverlapping block Jacobi preconditioner: it will help much!
469 %! for i = 1 : k : N # form 1-D Laplacian matrix 475 %!
470 %! M (i:i+k-1, i:i+k-1) = A (i:i+k-1, i:i+k-1); 476 %! M = zeros (N, N); k = 4;
471 %! endfor 477 %! for i = 1 : k : N # form 1-D Laplacian matrix
472 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); 478 %! M(i:i+k-1, i:i+k-1) = A(i:i+k-1, i:i+k-1);
473 %! printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1)); 479 %! endfor
474 %! semilogy ([0:iter], resvec(:,1),'o-b'); 480 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M);
475 %! legend('NO preconditioning: absolute residual', ... 481 %! printf ("BLOCK JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1));
476 %! 'JACOBI preconditioner: absolute residual', ... 482 %! semilogy ([0:iter], resvec(:,1), "o-b");
477 %! 'BLOCK JACOBI preconditioner: absolute residual'); 483 %! legend ("NO preconditioning: absolute residual", ...
478 %! hold off; 484 %! "JACOBI preconditioner: absolute residual", ...
485 %! "BLOCK JACOBI preconditioner: absolute residual");
486 %! hold off;
487
479 %!test 488 %!test
480 %! 489 %! # solve small diagonal system
481 %! #solve small diagonal system 490 %!
482 %! 491 %! N = 10;
483 %! N = 10; 492 %! A = diag ([1:N]); b = rand (N, 1);
484 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution 493 %! X = A \ b; # X is the true solution
485 %! [x, flag] = pcg (A, b, [], N+1); 494 %! [x, flag] = pcg (A, b, [], N+1);
486 %! assert(norm (x - X) / norm (X), 0, 1e-10); 495 %! assert (norm (x - X) / norm (X), 0, 1e-10);
487 %! assert(flag, 0); 496 %! assert (flag, 0);
488 %! 497
489 %!test 498 %!test
490 %! 499 %! # solve small indefinite diagonal system
491 %! #solve small indefinite diagonal system 500 %! # despite A is indefinite, the iteration continues and converges
492 %! #despite A is indefinite, the iteration continues and converges 501 %! # indefiniteness of A is detected
493 %! #indefiniteness of A is detected 502 %!
494 %! 503 %! N = 10;
495 %! N = 10; 504 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1);
496 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1); X = A \ b; #X is the true solution 505 %! X = A \ b; # X is the true solution
497 %! [x, flag] = pcg (A, b, [], N+1); 506 %! [x, flag] = pcg (A, b, [], N+1);
498 %! assert(norm (x - X) / norm (X), 0, 1e-10); 507 %! assert (norm (x - X) / norm (X), 0, 1e-10);
499 %! assert(flag, 3); 508 %! assert (flag, 3);
500 %! 509
501 %!test 510 %!test
502 %! 511 %! # solve tridiagonal system, do not converge in default 20 iterations
503 %! #solve tridiagonal system, do not converge in default 20 iterations 512 %!
504 %! 513 %! N = 100;
505 %! N = 100; 514 %! A = zeros (N, N);
506 %! A = zeros (N, N); 515 %! for i = 1 : N - 1 # form 1-D Laplacian matrix
507 %! for i = 1 : N - 1 # form 1-D Laplacian matrix 516 %! A(i:i+1, i:i+1) = [2 -1; -1 2];
508 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; 517 %! endfor
509 %! endfor 518 %! b = ones (N, 1);
510 %! b = ones (N, 1); X = A \ b; #X is the true solution 519 %! X = A \ b; # X is the true solution
511 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12); 520 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12);
512 %! assert(flag); 521 %! assert (flag);
513 %! assert(relres > 1.0); 522 %! assert (relres > 1.0);
514 %! assert(iter, 20); #should perform max allowable default number of iterations 523 %! assert (iter, 20); # should perform max allowable default number of iterations
515 %! 524
516 %!test 525 %!test
517 %! 526 %! # solve tridiagonal system with 'perfect' preconditioner
518 %! #solve tridiagonal system with 'prefect' preconditioner 527 %! # which converges in one iteration, so the eigest does not
519 %! #converges in one iteration, so the eigest does not work 528 %! # work and issues a warning
520 %! #and issues a warning 529 %!
521 %! 530 %! N = 100;
522 %! N = 100; 531 %! A = zeros (N, N);
523 %! A = zeros (N, N); 532 %! for i = 1 : N - 1 # form 1-D Laplacian matrix
524 %! for i = 1 : N - 1 # form 1-D Laplacian matrix 533 %! A (i:i+1, i:i+1) = [2 -1; -1 2];
525 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; 534 %! endfor
526 %! endfor 535 %! b = ones (N, 1);
527 %! b = ones (N, 1); X = A \ b; #X is the true solution 536 %! X = A \ b; # X is the true solution
528 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b); 537 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
529 %! assert(norm (x - X) / norm (X), 0, 1e-6); 538 %! assert (norm (x - X) / norm (X), 0, 1e-6);
530 %! assert(flag, 0); 539 %! assert (flag, 0);
531 %! assert(iter, 1); #should converge in one iteration 540 %! assert (iter, 1); # should converge in one iteration
532 %! assert(isnan (eigest), isnan ([NaN, NaN])); 541 %! assert (isnan (eigest), isnan ([NaN, NaN]));
533 %! 542