octave-antonio: scripts/sparse/pcr.m annotate

annotate scripts/sparse/pcr.m @ 14237:11949c9795a0

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

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

Mercurial > octave-antonio

annotate scripts/sparse/pcr.m @ 14237:11949c9795a0