octave: scripts/sparse/pcg.m annotate

annotate scripts/sparse/pcg.m @ 26376:00f796120a6d stable

maint: Update copyright dates in all source files.

author	John W. Eaton <jwe@octave.org>
date	Wed, 02 Jan 2019 16:32:43 -0500
parents	ddf8a03beffb
children	4e632e942e84

rev	line source
26376 00f796120a6d maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 25755 diff changeset	1 ## Copyright (C) 2004-2019 Piotr Krzyzanowski
00f796120a6d maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 25755 diff changeset	2 ## Copyright (C) 2016-2019 Cristiano Dorigo
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	3 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	4 ## This file is part of Octave.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	5 ##
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	6 ## Octave is free software: you can redistribute it and/or modify it
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	7 ## under the terms of the GNU General Public License as published by
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	8 ## the Free Software Foundation, either version 3 of the License, or
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	9 ## (at your option) any later version.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	10 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	11 ## 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	12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	14 ## GNU General Public License for more details.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	15 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	16 ## 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	17 ## along with Octave; see the file COPYING. If not, see
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	18 ## <https://www.gnu.org/licenses/>.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	19
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	20 ## -- texinfo --
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	21 ## @deftypefn {} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{})
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	22 ## @deftypefnx {} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, [], @var{x0}, @dots{})
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	23 ## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@var{A}, @var{b}, @dots{})
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	24 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	25 ## Solve the linear system of equations @w{@code{@var{A} * @var{x} = @var{b}}}
20164 df437a52bcaf doc: Update more docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19833 diff changeset	26 ## by means of the Preconditioned Conjugate Gradient iterative method.
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 19833 diff changeset	27 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	28 ## The input arguments are:
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	29 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	30 ## @itemize
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	31 ## @item @var{A} is the matrix of the linear system and it must be square.
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	32 ## @var{A} can be passed as a matrix, function handle, or inline function
40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	33 ## @code{Afun} such that @code{Afun(x) = A * x}. Additional parameters to
40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	34 ## @code{Afun} may be passed after @var{x0}.
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	35 ##
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	36 ## @var{A} has to be Hermitian and Positive Definite (@nospell{HPD})@. If
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	37 ## @code{pcg} detects @var{A} not to be positive definite, a warning is printed
40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	38 ## and the @var{flag} output is set.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	39 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	40 ## @item
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	41 ## @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	42 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	43 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	44 ## @var{tol} is the required relative tolerance for the residual error,
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	45 ## @w{@code{@var{b} - @var{A} * @var{x}}}. The iteration stops if
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	46 ## @w{@code{norm (@var{b} - @var{A} * @var{x})} @leq{}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	47 ## @w{@code{@var{tol} * norm (@var{b})}}}.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	48 ## If @var{tol} is omitted or empty, then a tolerance of 1e-6 is used.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	49 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	50 ## @item
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	51 ## @var{maxit} is the maximum allowed number of iterations; if @var{maxit}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	52 ## is omitted or empty then a value of 20 is used.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	53 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	54 ## @item
25003 2365c2661b3c doc: Spellcheck documentation ahead of 4.4 release. Rik <rik@octave.org> parents: 24992 diff changeset	55 ## @var{m} is a @nospell{HPD} preconditioning matrix. For any decomposition
24992 3d5a39077708 avoid some warnings from old versions of makeinfo (bug #53479) Guillaume Flandin parents: 24985 diff changeset	56 ## @code{@var{m} = @var{p1} * @var{p2}} such that
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	57 ## @w{@code{inv (@var{p1}) * @var{A} * inv (@var{p2})}} is @nospell{HPD}, the
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	58 ## conjugate gradient method is formally applied to the linear system
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	59 ## @w{@code{inv (@var{p1}) * @var{A} * inv (@var{p2}) * @var{y} = inv
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	60 ## (@var{p1}) * @var{b}}},
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	61 ## with @code{@var{x} = inv (@var{p2}) * @var{y}} (split preconditioning).
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	62 ## In practice, at each iteration of the conjugate gradient method a
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	63 ## linear system with matrix @var{m} is solved with @code{mldivide}.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	64 ## If a particular factorization
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	65 ## @code{@var{m} = @var{m1} * @var{m2}} is available (for instance, an
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	66 ## incomplete Cholesky factorization of @var{a}), the two matrices
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	67 ## @var{m1} and @var{m2} can be passed and the relative linear systems
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	68 ## are solved with the @code{mldivide} operator.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	69 ## Note that a proper choice of the preconditioner may dramatically improve
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	70 ## the overall performance of the method. Instead of matrices @var{m1} and
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	71 ## @var{m2}, the user may pass two functions which return the results of
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	72 ## applying the inverse of @var{m1} and @var{m2} to a vector.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	73 ## If @var{m1} is omitted or empty @code{[]}, then no preconditioning
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	74 ## is applied. If no factorization of @var{m} is available, @var{m2}
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	75 ## can be omitted or left [], and the input variable @var{m1} can be
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	76 ## used to pass the preconditioner @var{m}.
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
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	79 ## @var{x0} is the initial guess. If @var{x0} is omitted or empty then the
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	80 ## function sets @var{x0} to a zero vector by default.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	81 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	82 ##
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	83 ## The arguments which follow @var{x0} are treated as parameters, and passed in
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	84 ## an appropriate manner to any of the functions (@var{A} or @var{m1} or
40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	85 ## @var{m2}) that have been given to @code{pcg}.
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	86 ## See the examples below for further details.
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	87 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	88 ## The output arguments are:
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	89 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	90 ## @itemize
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	91 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	92 ## @var{x} is the computed approximation to the solution of
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	93 ## @w{@code{@var{A} * @var{x} = @var{b}}}. If the algorithm did not converge,
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	94 ## then @var{x} is the iteration which has the minimum residual.
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	95 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	96 ## @item
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	97 ## @var{flag} reports on the convergence:
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	98 ##
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	99 ## @itemize
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	100 ## @item 0: The algorithm converged to within the prescribed tolerance.
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	101 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	102 ## @item 1: The algorithm did not converge and it reached the maximum
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	103 ## number of iterations.
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	104 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	105 ## @item 2: The preconditioner matrix is singular.
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	106 ##
d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	107 ## @item 3: The algorithm stagnated, i.e., the absolute value of the
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	108 ## difference between the current iteration @var{x} and the previous is less
40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	109 ## than @code{@var{eps} * norm (@var{x},2)}.
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	110 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	111 ## @item 4: The algorithm detects that the input (preconditioned) matrix is not
25003 2365c2661b3c doc: Spellcheck documentation ahead of 4.4 release. Rik <rik@octave.org> parents: 24992 diff changeset	112 ## @nospell{HPD}.
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	113 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	114 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	115 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	116 ## @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	117 ## 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	118 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	119 ## @item
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	120 ## @var{iter} indicates the iteration of @var{x} which it was
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	121 ## computed. Since the output @var{x} corresponds to the minimal
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	122 ## residual solution, the total number of iterations that
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	123 ## the method performed is given by @code{length(resvec) - 1}.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	124 ##
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	125 ## @item
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	126 ## @var{resvec} describes the convergence history of the method.
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	127 ## @code{@var{resvec} (@var{i}, 1)} is the Euclidean norm of the residual, and
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	128 ## @code{@var{resvec} (@var{i}, 2)} is the preconditioned residual
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	129 ## norm, after the
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	130 ## (@var{i}-1)-th iteration, @code{@var{i} = 1, 2, @dots{}, @var{iter}+1}.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	131 ## The preconditioned residual norm is defined as
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	132 ## @code{@var{r}' * (@var{m} \ @var{r})} where
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	133 ## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	134 ## description of @var{m}. If @var{eigest} is not required, only
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	135 ## @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	136 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	137 ## @item
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	138 ## @var{eigest} returns the estimate for the smallest @code{@var{eigest}(1)}
f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	139 ## and largest @code{@var{eigest}(2)} eigenvalues of the preconditioned matrix
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	140 ## @w{@code{@var{P} = @var{m} \ @var{A}}}. In particular, if no
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	141 ## preconditioning is used, the estimates for the extreme eigenvalues of
f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	142 ## @var{A} are returned. @code{@var{eigest}(1)} is an overestimate and
f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	143 ## @code{@var{eigest}(2)} is an underestimate, so that
f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	144 ## @code{@var{eigest}(2) / @var{eigest}(1)} is a lower bound for
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	145 ## @code{cond (@var{P}, 2)}, which nevertheless in the limit should
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	146 ## 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	147 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	148 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	149 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	150 ## Let us consider a trivial problem with a tridiagonal matrix
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	151 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	152 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	153 ## @group
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	154 ## n = 10;
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	155 ## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	156 ## b = A * ones (n, 1);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	157 ## M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	158 ## M2 = M1';
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	159 ## M = M1 * M2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	160 ## Afun = @@(x) A * x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	161 ## Mfun = @@(x) M \ x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	162 ## M1fun = @@(x) M1 \ x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	163 ## M2fun = @@(x) M2 \ x;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	164 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	165 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	166 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	167 ## @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	168 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	169 ## @example
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	170 ## x = pcg (A, b)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	171 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	172 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	173 ## @sc{Example 2:} @code{pcg} with a function which computes
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	174 ## @code{@var{A} * @var{x}}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	175 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	176 ## @example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	177 ## x = pcg (Afun, b)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	178 ## @end example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	179 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	180 ## @sc{Example 3:} @code{pcg} with a preconditioner matrix @var{M}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	181 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	182 ## @example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	183 ## x = pcg (A, b, 1e-06, 100, M)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	184 ## @end example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	185 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	186 ## @sc{Example 4:} @code{pcg} with a function as preconditioner
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	187 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	188 ## @example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	189 ## x = pcg (Afun, b, 1e-6, 100, Mfun)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	190 ## @end example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	191 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	192 ## @sc{Example 5:} @code{pcg} with preconditioner matrices @var{M1}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	193 ## and @var{M2}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	194 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	195 ## @example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	196 ## x = pcg (A, b, 1e-6, 100, M1, M2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	197 ## @end example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	198 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	199 ## @sc{Example 6:} @code{pcg} with functions as preconditioners
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	200 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	201 ## @example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	202 ## x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	203 ## @end example
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	204 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	205 ## @sc{Example 7:} @code{pcg} with as input a function requiring an argument
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	206 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	207 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	208 ## @group
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	209 ## function y = Ap (A, x, p) # compute A^p * x
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	210 ## y = x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	211 ## for i = 1:p
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	212 ## y = A * y;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	213 ## endfor
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	214 ## endfunction
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	215 ## Apfun = @@(x, p) Ap (A, x, p);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	216 ## x = pcg (Apfun, b, [], [], [], [], [], 2);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	217 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	218 ## @end example
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	219 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	220 ## @sc{Example 8:} explicit example to show that @code{pcg} uses a
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	221 ## split preconditioner
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	222 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	223 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	224 ## @group
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	225 ## M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	226 ## M2 = M1';
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	227 ## M = M1 * M2;
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	228 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	229 ## ## reference solution computed by pcg after two iterations
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	230 ## [x_ref, fl] = pcg (A, b, [], 2, M)
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	231 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	232 ## ## split preconditioning
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	233 ## [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	234 ## x = M2 \ y # compare x and x_ref
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	235 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	236 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	237 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	238 ##
10791 3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	239 ## References:
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	240 ##
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	241 ## @enumerate
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	242 ## @item
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	243 ## 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	244 ## 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	245 ##
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	246 ## @item
19596 0e1f5a750d00 maint: Periodic merge of gui-release to default. John W. Eaton <jwe@octave.org> parents: 19148 19593 diff changeset	247 ## @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear Systems},
16826 a4969508008e doc: Periodic spellcheck of the documentation. Rik <rik@octave.org> parents: 14872 diff changeset	248 ## @nospell{PWS} 1996. (condition number estimate from PCG)
a4969508008e doc: Periodic spellcheck of the documentation. Rik <rik@octave.org> parents: 14872 diff changeset	249 ## Revised version of this book is available online at
25143 13fd0610480f doc: Use https whenever possible in @url entries. Rik <rik@octave.org> parents: 25054 diff changeset	250 ## @url{https://www-users.cs.umn.edu/~saad/books.html}
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	251 ## @end enumerate
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	252 ##
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	253 ## @seealso{sparse, pcr, gmres, bicg, bicgstab, cgs}
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	254 ## @end deftypefn
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	255
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	256 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	257 ## 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	258 ## - 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	259
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	260
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	261 function [x_min, flag, relres, iter_min, resvec, eigest] =...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	262 pcg (A, b, tol = [], maxit = [], M1 = [], M2 = [], x0 = [], varargin)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	263
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	264 ## Insert the default input (if necessary)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	265 [tol, maxit, x0] = __default__input__ ({1e-6, min(rows (b), 20),...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	266 zeros(size (b))}, tol, maxit, x0);
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	267
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	268 if (tol >= 1)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	269 warning ("Input tol is bigger than 1. \n Try to use a smaller tolerance.");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	270 elseif (tol <= eps / 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	271 warning ("Input tol may not be achievable by pcg. \n Try to use a bigger tolerance");
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8506 diff changeset	272 endif
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	273
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	274 ## Check if the input data A,b,m1,m2 are consistent (i.e. if they are
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	275 ## matrix or function handle)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	276
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	277 [Afun, M1fun, M2fun] = __alltohandles__ (A, b, M1, M2, "pcg");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	278
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	279 maxit += 2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	280 n_arg_out = nargout;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	281
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	282 ## Set Initial data
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	283 b_norm = norm (b);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	284 if (b_norm == 0)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	285 if (n_arg_out < 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	286 printf("The right hand side vector is all zero so pcg \n");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	287 printf ("returned an all zero solution without iterating.\n");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	288 endif
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	289 x_min = b;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	290 flag = 0;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	291 relres = 0;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	292 resvec = 0;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	293 iter_min = 0;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	294 eigest = [NaN, NaN];
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	295 return
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	296 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	297
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	298 x = x_pr = x_min = x0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	299
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	300 ## x_pr (x previous) needs to check the stagnation
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	301 ## x_min needs to save the iterated with minimum residual
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	302
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	303 r = b - feval (Afun, x, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	304 iter = 2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	305 iter_min = 0;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	306 flag = 1;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	307 resvec = zeros (maxit + 1, 2);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	308 resvec(1, 1) = norm (r);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	309 p = zeros (size (b));
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	310 alpha = old_tau = 1;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	311
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	312 if (n_arg_out > 5)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	313 T = zeros (maxit, maxit);
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	314 else
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	315 T = [];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	316 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	317
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	318 while (resvec(iter-1,1) > tol * b_norm && iter < maxit)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	319 if (iter == 2) # Check whether M1 or M2 are singular
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	320 try
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	321 warning ("error","Octave:singular-matrix","local")
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	322 z = feval (M1fun, r, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	323 z = feval (M2fun, z, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	324 catch
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	325 flag = 2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	326 break;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	327 end_try_catch
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	328 else
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	329 z = feval (M1fun, r, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	330 z = feval (M2fun, z, varargin{:});
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	331 endif
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	332
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	333 tau = z' * r;
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	334 resvec(iter - 1, 2) = sqrt (tau);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	335 beta = tau / old_tau;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	336 old_tau = tau;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	337 p = z + beta * p;
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	338 w = feval (Afun, p, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	339
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	340 ## Needed only for eigest.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	341
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	342 old_alpha = alpha;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	343 den = p' * w;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	344 alpha = tau / den;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	345
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	346 ## Check if alpha is negative and/or if it has a consistent
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	347 ## imaginary part: if yes then A probably is not positive definite
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	348 if ((abs (imag (tau)) >= abs (real (tau)) * tol) \|\| ...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	349 real (tau) <= 0 \|\| ...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	350 (abs (imag (den)) >= abs (real (den)) * tol) \|\| ...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	351 (real (den) <= 0))
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	352 flag = 4;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	353 break;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	354 endif
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	355
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	356 x += alpha * p;
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	357 r -= alpha * w;
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	358 resvec(iter, 1) = norm (r);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	359 ## Chek if the iterated has minimum residual
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	360 if (resvec (iter,1) <= resvec (iter_min + 1,1))
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	361 x_min = x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	362 iter_min = iter - 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	363 endif
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	364 if (n_arg_out > 5 && iter > 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	365 T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	366 [1, sqrt(beta); sqrt(beta), beta] ./ ...
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	367 old_alpha;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	368 endif
20735 418ae0cb752f Replace ++,-- with in-place operators for performance. Rik <rik@octave.org> parents: 20727 diff changeset	369 iter += 1;
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	370 if (norm (x - x_pr) <= eps * norm (x)) # Check the stagnation
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	371 flag = 3;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	372 break;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	373 endif
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	374 x_pr = x;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	375 endwhile
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	376
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	377 if (n_arg_out > 5)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	378 ## Apply the preconditioner once more and finish with the precond
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	379 ## residual.
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	380 z = feval (M1fun, r, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	381 z = feval (M2fun, z, varargin{:});
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	382 endif
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	383
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	384 ## (Eventually) computes the eigenvalue of inv(m2)inv(m1)A
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	385 if (n_arg_out > 5)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	386 if (flag != 4)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	387 if (iter > 3)
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	388 T = T(2:iter-2,2:iter-2);
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	389 l = eig (T);
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	390 eigest = [min(l), max(l)];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	391 else
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	392 eigest = [NaN, NaN];
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	393 warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	394 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	395 else
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	396 eigest = [NaN, NaN];
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	397 warning ('pcg: eigenvalue estimate failed: matrix not positive definite?')
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	398 endif
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	399 resvec(iter - 1, 2) = sqrt (r' * z);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	400 resvec = resvec (1:(iter-1), :);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	401 else
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	402 eigest = [NaN, NaN];
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	403 resvec = resvec(1:(iter-1),1);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	404 endif
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	405
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	406 ## Set the last variables
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	407
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	408 if (flag == 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	409 relres = 1;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	410 elseif (resvec (1, 1) == 0)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	411 relres = 0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	412 else
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	413 relres = resvec(iter_min+1, 1) ./ resvec(1, 1);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	414 endif
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	415
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	416 iter -= 2; # compatibility
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	417
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	418 ## Set the flag in the proper way if flag not 3, 4 or 2
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	419 if (flag == 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	420 flag = 2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	421 elseif (flag == 1) && (relres <= tol)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	422 flag = 0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	423 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	424
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	425 if (n_arg_out < 2)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	426 switch (flag)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	427 case {0}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	428 printf ("pcg converged at iteration %d ", iter_min);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	429 printf ("with relative residual %d\n", relres);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	430 case {1}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	431 printf ("pcg stopped at iteration %d ", iter+1);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	432 printf ("without converging to the desired tolerance %d ", tol);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	433 printf ("because the maximum number of iteration was reached, \n");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	434 printf ("The iterated returned (number %d) ",iter_min);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	435 printf ("has relative residual %d \n", relres);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	436 case {2}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	437 printf ("pcg stopped at iteration %d ", iter+1)
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	438 printf ("without converging to the desired tolerance %d ", tol);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	439 printf ("because the preconditioned matrix is singular.\n");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	440 printf ("The iterated returned (number %d) ", iter_min);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	441 printf ("has relative residual %d \n", relres);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	442 case {3}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	443 printf ("pcg stopped at iteration %d ", iter+1);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	444 printf ("without converging to the desired tolerance %d ", tol);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	445 printf ("because of stagnation. \n");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	446 printf ("The iterated returned (number %d) ", iter_min);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	447 printf ("has relative residual %d.\n", relres);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	448 case {4}
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	449 printf ("pcg stopped at iteration %d ", iter + 1);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	450 printf ("without converging to the desired tolerance %d ",tol);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	451 printf ("because the (preconditioned) matrix is not positive definite. \n");
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	452 printf ("The iterate returned (number %d) ", iter_min);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	453 printf ("has relative residual %d \n", relres);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	454 endswitch
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	455 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	456 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	457
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	458
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	459 %!demo # simplest use
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	460 %! n = 10;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	461 %! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	462 %! b = A * ones (n, 1);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	463 %! M1 = ichol (A); # for this tridiagonal case it corresponds to chol (A)'
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	464 %! M2 = M1';
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	465 %! M = M1 * M2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	466 %! x = pcg (A, b);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	467 %! Afun = @(x) A * x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	468 %! x = pcg (Afun, b);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	469 %! x = pcg (A, b, 1e-6, 100, M);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	470 %! x = pcg (A, b, 1e-6, 100, M1, M2);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	471 %! Mfun = @(x) M \ x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	472 %! x = pcg (Afun, b, 1e-6, 100, Mfun);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	473 %! M1fun = @(x) M1 \ x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	474 %! M2fun = @(x) M2 \ x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	475 %! x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	476 %! function y = Ap (A, x, p) # compute A^p * x
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	477 %! y = x;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	478 %! for i = 1:p
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	479 %! y = A * y;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	480 %! endfor
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	481 %! endfunction
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	482 %! Afun = @(x, p) Ap (A, x, p);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	483 %! ## solution of A^2 * x = b
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	484 %! x = pcg (Afun, b, [], [], [], [], [], 2);
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	485
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	486 %!demo
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	487 %! n = 10;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	488 %! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	489 %! b = A * ones (n, 1);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	490 %! M1 = ichol (A + 0.1 * eye (n)); # Perturb the factorization of A
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	491 %! M2 = M1';
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	492 %! M = M1 * M2;
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	493 %!
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	494 %! ## Reference solution computed by pcg after two iterations
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	495 %! [x_ref, fl] = pcg (A, b, [], 2, M);
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	496 %! x_ref
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	497 %!
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	498 %! ## Split preconditioning
24860 6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	499 %! [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	500 %! x = M2 \ y # compare x and x_ref
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	501 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	502 %! ## Check that all type of inputs work
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	503 %! A = toeplitz (sparse ([2, 1 ,0, 0, 0]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	504 %! b = A * ones (5, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	505 %! M1 = diag (sqrt (diag (A)));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	506 %! M2 = M1; # M1 * M2 is the Jacobi preconditioner
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	507 %! Afun = @(z) A*z;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	508 %! M1_fun = @(z) M1 \ z;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	509 %! M2_fun = @(z) M2 \ z;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	510 %! [x, flag, ~, iter] = pcg (A,b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	511 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	512 %! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	513 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	514 %! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	515 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	516 %! [x, flag] = pcg (A, b, [], [], M1_fun, M2_fun);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	517 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	518 %! [x, flag] = pcg (A, b,[],[], M1_fun, M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	519 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	520 %! [x, flag] = pcg (A, b,[],[], M1, M2_fun);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	521 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	522 %! [x, flag] = pcg (Afun, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	523 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	524 %! [x, flag] = pcg (Afun, b,[],[], M1 * M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	525 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	526 %! [x, flag] = pcg (Afun, b,[],[], M1, M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	527 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	528 %! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	529 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	530 %! [x, flag] = pcg (Afun, b,[],[], M1, M2_fun);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	531 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	532 %! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2_fun);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	533 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	534
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	535 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	536 %! ## solve a small diagonal system
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	537 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	538 %! A = diag ([1:N]); b = rand (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	539 %! [x, flag] = pcg (A, b, [], N+1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	540 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	541 %! assert (norm (b - A*x) / norm (b), 0, 1e-6);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	542
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	543 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	544 %! ## A is not positive definite
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	545 %! ## The indefiniteness of A is detected.
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	546 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	547 %! A = -diag ([1:N]); b = sum (A, 2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	548 %! [x, flag] = pcg (A, b, [], N + 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	549 %! assert (flag, 4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	550
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	551 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	552 %! ## solve tridiagonal system, do not converge in default 20 iterations
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	553 %! N = 100;
25755 ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	554 %! ## Form 1-D Laplacian matrix
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	555 %! A = 2 * eye (N,N);
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	556 %! A(2:(N+1):end) = -1;
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	557 %! A((N+1):(N+1):end) = -1;
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	558 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	559 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	560 %! assert (flag);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	561 %! assert (relres >= 1.0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	562
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	563 %!warning <iteration converged too fast>
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	564 %! ## solve tridiagonal system with "perfect" preconditioner which converges
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	565 %! ## in one iteration, so the eigest does not work and issues a warning.
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	566 %! N = 100;
25755 ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	567 %! ## Form 1-D Laplacian matrix
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	568 %! A = 2 * eye (N,N);
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	569 %! A(2:(N+1):end) = -1;
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	570 %! A((N+1):(N+1):end) = -1;
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	571 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	572 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	573 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	574 %! assert (norm (b - A*x) / norm (b), 0, 1e-6);
25436 996d78102a71 maint: Strip trailing whitespace from source files. John W. Eaton <jwe@octave.org> parents: 25151 diff changeset	575 %!
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	576 %! assert (isnan (eigest), isnan ([NaN, NaN]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	577
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	578 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	579 %! ## pcg detect a non-Hermitian matrix, with a considerable imaginary part.
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	580 %! ## In this example, Matlab does not recognize the wrong type of matrix and
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	581 %! ## makes iterations until it reaches maxit.
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	582 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	583 %! A = diag (1:N) + 1e-4*i;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	584 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	585 %! [x, flag] = pcg (A, b, []);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	586 %! assert (flag, 4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	587
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	588 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	589 %! ## The imaginary part is not influent (it is too small), so pcg doesn't stop
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	590 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	591 %! A = diag (1:N) + 1e-10*i;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	592 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	593 %! [x, flag] = pcg (A, b, [], N+1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	594 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	595 %! assert (norm (b - A*x) / norm (b), 0, 1e-6);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	596
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	597 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	598 %! ## pcg solves linear system with A Hermitian positive definite
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	599 %! N = 20;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	600 %! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ...
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	601 %! i * sparse (toeplitz ([0, 1, zeros(1, 18)], [0, -1, zeros(1,18)]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	602 %! b = A * ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	603 %! Hermitian_A = ishermitian (A);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	604 %! [x, flag] = pcg (A, b, [], 2*N);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	605 %! assert (Hermitian_A, true);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	606 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	607 %! assert (x, ones (N, 1), -1e-4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	608
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	609 %!testif HAVE_CHOLMOD
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	610 %! ## pcg solves preconditioned linear system with A HPD
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	611 %! N = 20;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	612 %! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ...
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	613 %! i * sparse (toeplitz ([0, 1, zeros(1, 18)], [0, -1, zeros(1,18)]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	614 %! b = A * ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	615 %! M2 = chol (A + 0.1 * eye (N)); # Factor of a perturbed matrix
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	616 %! M = M2' * M2;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	617 %! Hermitian_A = ishermitian (A);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	618 %! Hermitian_M = ishermitian (M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	619 %! [x, flag] = pcg (A, b, [], 2*N, M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	620 %! assert (Hermitian_A, true);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	621 %! assert (Hermitian_M, true);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	622 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	623 %! assert (x, ones (N, 1), -1e-4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	624
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	625 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	626 %! ## pcg recognizes that the preconditioner matrix is singular
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	627 %! N = 3;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	628 %! A = toeplitz ([2, 1, 0]);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	629 %! M = [1 0 0; 0 1 0; 0 0 0]; # the last row is zero
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	630 %! [x, flag] = pcg (A, ones (3, 1), [], [], M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	631 %! assert (flag, 2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	632
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	633 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	634 %! A = rand (4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	635 %! A = A' * A;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	636 %! [x, flag] = pcg (A, zeros (4, 1), [], [], [], [], ones (4, 1));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	637 %! assert (x, zeros (4, 1));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	638
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	639 ## Test return types
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	640 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	641 %! A = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	642 %! b = 1;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	643 %! [x, flag] = pcg (A, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	644 %! assert (class (x), "single");
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	645
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	646 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	647 %! A = 1;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	648 %! b = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	649 %! [x, flag] = pcg (A, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	650 %! assert (class (x), "single");
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	651
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	652 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	653 %! A = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	654 %! b = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	655 %! [x, flag] = pcg (A, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	656 %! assert (class (x), "single");
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	657
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	658 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	659 %!function y = Afun (x)
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	660 %! A = toeplitz ([2, 1, 0, 0]);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	661 %! y = A * x;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	662 %!endfunction
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	663 %! [x, flag] = pcg ("Afun", [3; 4; 4; 3]);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	664 %! assert (x, ones (4, 1), 1e-6);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	665
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	666 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	667 %! ## unpreconditioned residual
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	668 %! A = toeplitz (sparse ([4, 1, 0, 0, 0]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	669 %! b = sum (A, 2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	670 %! M = toeplitz (sparse ([2, 1, 0, 0, 0]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	671 %! [x, flag, relres] = pcg (A, b, [], 2, M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	672 %! assert (norm (b - A * x) / norm (b), relres, 8 * eps);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	673