octave-lyh: scripts/sparse/pcg.m annotate

annotate scripts/sparse/pcg.m @ 16508:f19e24c97b20

move common warndlg, errordlg, helpdlg, and msgbox code to private function * message_dialog.m: New file. * scripts/ui/module.mk: Include it in the list of functions. * errordlg.m, helpdlg.m, warndlg.m, msgbox.m: Call message_dialog to do most of the work.

author	John W. Eaton <jwe@octave.org>
date	Fri, 12 Apr 2013 18:17:26 -0400
parents	c2dbdeaa25df
children	a4969508008e

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: 7007 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: 7007 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: 7007 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: 7007 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: 10846 diff changeset	20 ## @deftypefn {Function File} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @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}, @var{eigest}] =} pcg (@dots{})
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	22 ##
14093 050bc580cb60 doc: Various docstring improvements before 3.6.0 release. Rik <octave@nomad.inbox5.com> parents: 11588 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 Gradient 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: 10846 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: 10846 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: 10846 diff changeset	32 ## @var{A} should be symmetric and positive definite; if @code{pcg}
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	33 ## finds @var{A} to not be positive definite, 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,
11563 3c6e8aaa9555 Grammarcheck m-files before 3.4 release. Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	41 ## @code{@var{b} - @var{A} * @var{x}}. The iteration stops if
11587 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{pcg} 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.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	51 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	52 ## @item
10846 a4f482e66b65 Grammarcheck more of the documentation. Rik <octave@nomad.inbox5.com> parents: 10793 diff changeset	53 ## @var{m} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that
a4f482e66b65 Grammarcheck more of the documentation. Rik <octave@nomad.inbox5.com> parents: 10793 diff changeset	54 ## the iteration is (theoretically) equivalent to solving by @code{pcg}
a4f482e66b65 Grammarcheck more of the documentation. Rik <octave@nomad.inbox5.com> parents: 10793 diff changeset	55 ## @code{@var{P} *
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	56 ## @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	57 ## 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	58 ## improve the overall performance of the method. Instead of matrices
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	59 ## @var{m1} and @var{m2}, the user may pass two functions which return
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	60 ## the results of applying the inverse of @var{m1} and @var{m2} to
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	61 ## a vector (usually this is the preferred way of using the preconditioner).
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	62 ## If @code{[]} is supplied for @var{m1}, or @var{m1} is omitted, no
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	63 ## preconditioning is applied. If @var{m2} is omitted, @var{m} = @var{m1}
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	64 ## will be used as preconditioner.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	65 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	66 ## @item
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	67 ## @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	68 ## function sets @var{x0} to a zero vector by default.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	69 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	70 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	71 ## 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: 10846 diff changeset	72 ## 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	73 ## which are passed to @code{pcg}. See the examples below for further
1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	74 ## details. The output arguments are
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	75 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	76 ## @itemize
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	77 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	78 ## @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: 10846 diff changeset	79 ## @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	80 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	81 ## @item
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	82 ## @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	83 ## 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	84 ## 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	85 ## for the iteration count was reached. @code{@var{flag} = 3} reports that
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	86 ## the (preconditioned) matrix was found not positive definite.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	87 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	88 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	89 ## @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	90 ## 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	91 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	92 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	93 ## @var{iter} is the actual number of iterations performed.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	94 ##
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	95 ## @item
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	96 ## @var{resvec} describes the convergence history of the method.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	97 ## @code{@var{resvec} (i,1)} is the Euclidean norm of the residual, and
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	98 ## @code{@var{resvec} (i,2)} is the preconditioned residual norm,
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	99 ## after the (@var{i}-1)-th iteration, @code{@var{i} =
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	100 ## 1, 2, @dots{}, @var{iter}+1}. The preconditioned residual norm
6547 4fb053f24fd6 [project @ 2007-04-19 21:47:40 by jwe] jwe parents: 6498 diff changeset	101 ## is defined as
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	102 ## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	103 ## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	104 ## description of @var{m}. If @var{eigest} is not required, only
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	105 ## @code{@var{resvec} (:,1)} is returned.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	106 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	107 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	108 ## @var{eigest} returns the estimate for the smallest @code{@var{eigest}
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	109 ## (1)} and largest @code{@var{eigest} (2)} eigenvalues of the
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	110 ## preconditioned matrix @code{@var{P} = @var{m} \ @var{A}}. In
7007 6304d9ea0a30 [project @ 2007-10-11 16:26:36 by jwe] jwe parents: 6747 diff changeset	111 ## particular, if no preconditioning is used, the estimates for the
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	112 ## extreme eigenvalues of @var{A} are returned. @code{@var{eigest} (1)}
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	113 ## is an overestimate and @code{@var{eigest} (2)} is an underestimate,
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	114 ## so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	115 ## for @code{cond (@var{P}, 2)}, which nevertheless in the limit should
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	116 ## theoretically be equal to the actual value of the condition number.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	117 ## The method which computes @var{eigest} works only for symmetric positive
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	118 ## definite @var{A} and @var{m}, and the user is responsible for
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	119 ## verifying this assumption.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	120 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	121 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	122 ## 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	123 ## sparsity of A)
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	124 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	125 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	126 ## @group
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	127 ## n = 10;
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	128 ## A = diag (sparse (1:n));
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	129 ## b = rand (n, 1);
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	130 ## [l, u, p, q] = luinc (A, 1.e-3);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	131 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	132 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	133 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	134 ## @sc{Example 1:} Simplest use of @code{pcg}
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	135 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	136 ## @example
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	137 ## x = pcg (A,b)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	138 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	139 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	140 ## @sc{Example 2:} @code{pcg} with a function which computes
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	141 ## @code{@var{A} * @var{x}}
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	142 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	143 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	144 ## @group
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	145 ## function y = apply_a (x)
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	146 ## y = [1:N]' .* x;
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	147 ## endfunction
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	148 ##
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	149 ## x = pcg ("apply_a", b)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	150 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	151 ## @end example
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	152 ##
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	153 ## @sc{Example 3:} @code{pcg} with a preconditioner: @var{l} * @var{u}
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	154 ##
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	155 ## @example
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	156 ## x = pcg (A, b, 1.e-6, 500, l*u)
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	157 ## @end example
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	158 ##
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	159 ## @sc{Example 4:} @code{pcg} with a preconditioner: @var{l} * @var{u}.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	160 ## Faster than @sc{Example 3} since lower and upper triangular matrices
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	161 ## are easier to invert
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	162 ##
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	163 ## @example
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	164 ## x = pcg (A, b, 1.e-6, 500, l, u)
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	165 ## @end example
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	166 ##
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	167 ## @sc{Example 5:} Preconditioned iteration, with full diagnostics. The
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	168 ## 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: 10846 diff changeset	169 ## @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	170 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	171 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	172 ## @group
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	173 ## function y = apply_m (x)
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	174 ## k = floor (length (x) - 2);
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	175 ## y = x;
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	176 ## y(1:k) = x(1:k) ./ [1:k]';
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	177 ## endfunction
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	178 ##
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	179 ## [x, flag, relres, iter, resvec, eigest] = ...
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	180 ## pcg (A, b, [], [], "apply_m");
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	181 ## semilogy (1:iter+1, resvec);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	182 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	183 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	184 ##
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	185 ## @sc{Example 6:} Finally, a preconditioner which depends on a
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	186 ## parameter @var{k}.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	187 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	188 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	189 ## @group
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	190 ## 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	191 ## K = varargin@{1@};
6547 4fb053f24fd6 [project @ 2007-04-19 21:47:40 by jwe] jwe parents: 6498 diff changeset	192 ## y = x;
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	193 ## y(1:K) = x(1:K) ./ [1:K]';
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	194 ## endfunction
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	195 ##
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	196 ## [x, flag, relres, iter, resvec, eigest] = ...
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	197 ## pcg (A, b, [], [], "apply_m", [], [], 3)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	198 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	199 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	200 ##
10791 3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	201 ## References:
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	202 ##
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	203 ## @enumerate
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	204 ## @item
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	205 ## C.T. Kelley, @cite{Iterative Methods for Linear and Nonlinear Equations},
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	206 ## SIAM, 1995. (the base PCG algorithm)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	207 ##
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	208 ## @item
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	209 ## Y. Saad, @cite{Iterative Methods for Sparse Linear Systems}, PWS 1996.
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	210 ## (condition number estimate from PCG) Revised version of this book is
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	211 ## available online at @url{http://www-users.cs.umn.edu/~saad/books.html}
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	212 ## @end enumerate
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	213 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	214 ## @seealso{sparse, pcr}
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	215 ## @end deftypefn
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	216
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	217 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	218 ## Modified by: Vittoria Rezzonico <vittoria.rezzonico@epfl.ch>
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	219 ## - Add the ability to provide the pre-conditioner as two separate matrices
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	220
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	221 function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, tol, maxit, m1, m2, x0, varargin)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	222
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	223 ## M = M1*M2
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	224
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	225 if (nargin < 7 \|\| isempty (x0))
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	226 x = zeros (size (b));
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	227 else
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	228 x = x0;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	229 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	230
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	231 if (nargin < 5 \|\| isempty (m1))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	232 exist_m1 = 0;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	233 else
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	234 exist_m1 = 1;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	235 endif
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	236
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	237 if (nargin < 6 \|\| isempty (m2))
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	238 exist_m2 = 0;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	239 else
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	240 exist_m2 = 1;
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	241 endif
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	242
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	243 if (nargin < 4 \|\| isempty (maxit))
14872 c2dbdeaa25df maint: use rows() and columns() to clarify m-files. Rik <octave@nomad.inbox5.com> parents: 14868 diff changeset	244 maxit = min (rows (b), 20);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	245 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	246
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	247 maxit += 2;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	248
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	249 if (nargin < 3 \|\| isempty (tol))
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	250 tol = 1e-6;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	251 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	252
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	253 preconditioned_residual_out = false;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	254 if (nargout > 5)
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	255 T = zeros (maxit, maxit);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	256 preconditioned_residual_out = true;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	257 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	258
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	259 ## Assume A is positive definite.
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	260 matrix_positive_definite = true;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	261
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	262 p = zeros (size (b));
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	263 oldtau = 1;
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	264 if (isnumeric (A))
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	265 ## A is a matrix.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	266 r = b - A*x;
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	267 else
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	268 ## A should be a function.
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	269 r = b - feval (A, x, varargin{:});
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	270 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	271
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	272 resvec(1,1) = norm (r);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	273 alpha = 1;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	274 iter = 2;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	275
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	276 while (resvec (iter-1,1) > tol * resvec (1,1) && iter < maxit)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	277 if (exist_m1)
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	278 if (isnumeric (m1))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	279 y = m1 \ r;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	280 else
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	281 y = feval (m1, r, varargin{:});
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	282 endif
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	283 else
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	284 y = r;
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	285 endif
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	286 if (exist_m2)
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	287 if (isnumeric (m2))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	288 z = m2 \ y;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	289 else
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	290 z = feval (m2, y, varargin{:});
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	291 endif
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	292 else
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	293 z = y;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	294 endif
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	295 tau = z' * r;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	296 resvec (iter-1,2) = sqrt (tau);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	297 beta = tau / oldtau;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	298 oldtau = tau;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	299 p = z + beta * p;
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	300 if (isnumeric (A))
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	301 ## A is a matrix.
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	302 w = A * p;
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	303 else
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	304 ## A should be a function.
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 10846 diff changeset	305 w = feval (A, p, varargin{:});
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	306 endif
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	307 ## Needed only for eigest.
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	308 oldalpha = alpha;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	309 alpha = tau / (p'*w);
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	310 if (alpha <= 0.0)
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	311 ## Negative matrix.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	312 matrix_positive_definite = false;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	313 endif
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	314 x += alpha * p;
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	315 r -= alpha * w;
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	316 if (nargout > 5 && iter > 2)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	317 T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	318 [1 sqrt(beta); sqrt(beta) beta]./oldalpha;
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	319 ## EVS = eig (T(2:iter-1,2:iter-1));
5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	320 ## fprintf (stderr,"PCG condest: %g (iteration: %d)\n", max (EVS)/min (EVS),iter);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	321 endif
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	322 resvec (iter,1) = norm (r);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	323 iter++;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	324 endwhile
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	325
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	326 if (nargout > 5)
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	327 if (matrix_positive_definite)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	328 if (iter > 3)
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	329 T = T(2:iter-2,2:iter-2);
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	330 l = eig (T);
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	331 eigest = [min(l), max(l)];
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	332 ## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1));
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	333 else
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	334 eigest = [NaN, NaN];
11588 d5bd2766c640 style fixes for warning and error messages in script files John W. Eaton <jwe@octave.org> parents: 11587 diff changeset	335 warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	336 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	337 else
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	338 eigest = [NaN, NaN];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	339 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	340
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	341 ## Apply the preconditioner once more and finish with the precond
bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	342 ## residual.
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	343 if (exist_m1)
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	344 if (isnumeric (m1))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	345 y = m1 \ r;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	346 else
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	347 y = feval (m1, r, varargin{:});
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	348 endif
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	349 else
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	350 y = r;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	351 endif
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	352 if (exist_m2)
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	353 if (isnumeric (m2))
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	354 z = m2 \ y;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	355 else
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	356 z = feval (m2, y, varargin{:});
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	357 endif
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	358 else
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	359 z = y;
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	360 endif
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	361
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	362 resvec (iter-1,2) = sqrt (r' * z);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	363 else
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	364 resvec = resvec(:,1);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	365 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	366
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	367 flag = 0;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	368 relres = resvec (iter-1,1) ./ resvec(1,1);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	369 iter -= 2;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	370 if (iter >= maxit - 2)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	371 flag = 1;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	372 if (nargout < 2)
6498 2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe] jwe parents: 5838 diff changeset	373 warning ("pcg: maximum number of iterations (%d) reached\n", iter);
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	374 warning ("the initial residual norm was reduced %g times.\n", ...
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	375 1.0 / relres);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	376 endif
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	377 elseif (nargout < 2)
6498 2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe] jwe parents: 5838 diff changeset	378 fprintf (stderr, "pcg: converged in %d iterations. ", iter);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe] jwe parents: 5838 diff changeset	379 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	380 1.0/relres);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	381 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	382
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	383 if (! matrix_positive_definite)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	384 flag = 3;
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	385 if (nargout < 2)
6498 2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe] jwe parents: 5838 diff changeset	386 warning ("pcg: matrix not positive definite?\n");
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	387 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	388 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	389 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	390
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	391
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	392 %!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	393 %! # Simplest usage of pcg (see also 'help pcg')
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	394 %!
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	395 %! 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	396 %! A = diag ([1: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	397 %! 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	398 %! x = pcg (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	399 %! printf ("The solution relative error is %g\n", norm (x - y) / norm (y));
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	400 %!
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	401 %! # You shouldn't be afraid if pcg 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	402 %! # 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	403 %! # of pcg 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	404
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	405 %!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	406 %! # Full output from pcg, except for the eigenvalue estimates
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	407 %! # We use this output to plot the 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	408 %!
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	409 %! 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	410 %! A = diag ([1: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	411 %! 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	412 %! [x, flag, relres, iter, resvec] = pcg (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	413 %! 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	414 %! 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	415 %! semilogy ([0:iter], resvec / resvec(1), "o-g");
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	416 %! 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	417 %! legend ("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	418
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	419 %!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	420 %! # Full output from pcg, including the eigenvalue estimates
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	421 %! # Hilbert matrix is extremely ill-conditioned, so pcg WILL have problems
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	422 %!
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	423 %! 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	424 %! A = hilb (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	425 %! 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	426 %! [x, flag, relres, iter, resvec, eigest] = pcg (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	427 %! 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	428 %! printf ("Condition number estimate is %g\n", eigest(2) / eigest(1));
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	429 %! printf ("Actual 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	430 %! 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	431 %! semilogy ([0:iter], resvec, ["o-g";"+-r"]);
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	432 %! 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	433 %! legend ("absolute residual", "absolute preconditioned residual");
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	434
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	435 %!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	436 %! # Full output from pcg, including the eigenvalue estimates
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	437 %! # We use the 1-D Laplacian matrix for A, and 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	438 %! # and 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	439 %! # 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	440 %! # which is the diagonal of A
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	441 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	442 %! 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	443 %! 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	444 %! 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	445 %! 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	446 %! endfor
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	447 %! 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	448 %! 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	449 %! 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	450 %! 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	451 %! # No preconditioner: the convergence is very slow!
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	452 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	453 %! [x, flag, relres, iter, resvec, eigest] = pcg (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	454 %! printf ("System condition number estimate is %g\n", eigest(2) / eigest(1));
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	455 %! 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	456 %! semilogy ([0:iter], resvec(:,1), "o-g");
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	457 %! 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	458 %! legend ("NO preconditioning: absolute residual");
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	459 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	460 %! 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	461 %! # Test Jacobi preconditioner: it will not help much!!!
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	462 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	463 %! 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	464 %! [x, flag, relres, iter, resvec, eigest] = pcg (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	465 %! printf ("JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1));
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	466 %! 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	467 %! semilogy ([0:iter], resvec(:,1), "o-r");
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	468 %! legend ("NO preconditioning: 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	469 %! "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	470 %!
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	471 %! 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	472 %! # Test nonoverlapping block Jacobi preconditioner: it will help much!
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	473 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	474 %! 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	475 %! for i = 1 : k : N # 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	476 %! 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	477 %! endfor
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	478 %! [x, flag, relres, iter, resvec, eigest] = pcg (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	479 %! printf ("BLOCK JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1));
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	480 %! semilogy ([0:iter], resvec(:,1), "o-b");
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	481 %! legend ("NO preconditioning: 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	482 %! "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	483 %! "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	484 %! 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	485
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	486 %!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	487 %! # solve small diagonal system
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	488 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	489 %! 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	490 %! A = diag ([1: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	491 %! 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	492 %! [x, flag] = pcg (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	493 %! assert (norm (x - X) / norm (X), 0, 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	494 %! 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	495
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	496 %!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	497 %! # solve small indefinite diagonal system
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	498 %! # despite A is indefinite, the iteration continues and converges
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	499 %! # indefiniteness of A is detected
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	500 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	501 %! 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	502 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); 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	503 %! 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	504 %! [x, flag] = pcg (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	505 %! assert (norm (x - X) / norm (X), 0, 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	506 %! assert (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	507
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	508 %!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	509 %! # solve tridiagonal system, do not converge in default 20 iterations
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	510 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	511 %! 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	512 %! 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	513 %! 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	514 %! 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	515 %! endfor
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	516 %! 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	517 %! 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	518 %! [x, flag, relres, iter, resvec, eigest] = pcg (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	519 %! assert (flag);
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	520 %! assert (relres > 1.0);
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	521 %! assert (iter, 20); # should perform max allowable default number of iterations
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	522
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	523 %!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	524 %! # 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	525 %! # which converges in one iteration, so the eigest does not
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	526 %! # work and issues a warning
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	527 %!
14237 11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	528 %! 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	529 %! 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	530 %! 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	531 %! 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	532 %! endfor
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	533 %! 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	534 %! 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	535 %! [x, flag, relres, iter, resvec, eigest] = pcg (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	536 %! assert (norm (x - X) / norm (X), 0, 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	537 %! 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	538 %! 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	539 %! assert (isnan (eigest), isnan ([NaN, NaN]));
11949c9795a0 Revamp %!demos in m-files to use Octave coding conventions on spacing, etc. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	540

Mercurial > hg > octave-lyh

annotate scripts/sparse/pcg.m @ 16508:f19e24c97b20