octave: scripts/sparse/pcg.m annotate

annotate scripts/sparse/pcg.m @ 31706:597f3ee61a48 stable

update Octave Project Developers copyright for the new year

author	John W. Eaton <jwe@octave.org>
date	Fri, 06 Jan 2023 13:11:27 -0500
parents	e1788b1a315f
children	fe178c793ea0

rev	line source
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	1 ########################################################################
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	2 ##
31706 597f3ee61a48 update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 30893 diff changeset	3 ## Copyright (C) 2004-2023 The Octave Project Developers
27918 b442ec6dda5c use centralized file for copyright info for individual contributors John W. Eaton <jwe@octave.org> parents: 27356 diff changeset	4 ##
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	6 ## distribution or <https://octave.org/copyright/>.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	7 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	8 ## This file is part of Octave.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	9 ##
24534 194eb4bd202b maint: Update punctuation for GPL v3 license text. Rik <rik@octave.org> parents: 23220 diff changeset	10 ## 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	11 ## 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	12 ## 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	13 ## (at your option) any later version.
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	14 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	15 ## 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	16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755 3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting. Rik <rik@octave.org> parents: 22323 diff changeset	18 ## GNU General Public License for more details.
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 ## 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	21 ## 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	22 ## <https://www.gnu.org/licenses/>.
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	23 ##
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	24 ########################################################################
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	25
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	26 ## -- texinfo --
24985 d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	27 ## @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	28 ## @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	29 ## @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	30 ##
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 ## 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	32 ## 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	33 ##
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	34 ## The input arguments are:
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	35 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	36 ## @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	37 ## @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	38 ## @var{A} can be passed as a matrix, function handle, or inline function
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	39 ## @code{Afcn} such that @code{Afcn(x) = A * x}. Additional parameters to
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	40 ## @code{Afcn} 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	41 ##
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	42 ## @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	43 ## @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	44 ## 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	45 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	46 ## @item
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	47 ## @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	48 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	49 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	50 ## @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	51 ## @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	52 ## @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	53 ## @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	54 ## 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	55 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	56 ## @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	57 ## @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	58 ## 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	59 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	60 ## @item
25003 2365c2661b3c doc: Spellcheck documentation ahead of 4.4 release. Rik <rik@octave.org> parents: 24992 diff changeset	61 ## @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	62 ## @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	63 ## @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	64 ## 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	65 ## @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	66 ## (@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	67 ## 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	68 ## 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	69 ## 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	70 ## 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	71 ## @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	72 ## 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	73 ## @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	74 ## 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	75 ## 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	76 ## 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	77 ## @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	78 ## 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	79 ## 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	80 ## 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	81 ## 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	82 ## 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	83 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	84 ## @item
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	85 ## @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	86 ## function sets @var{x0} to a zero vector by default.
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	87 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	88 ##
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	89 ## 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	90 ## 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	91 ## @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	92 ## 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	93 ##
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	94 ## The output arguments are:
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	95 ##
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	96 ## @itemize
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	97 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	98 ## @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	99 ## @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	100 ## 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	101 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	102 ## @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	103 ## @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	104 ##
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	105 ## @itemize
25151 40466945a451 pcg.m: Tweaks to documentation for clarity. Rik <rik@octave.org> parents: 25148 diff changeset	106 ## @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	107 ##
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	108 ## @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	109 ## number of iterations.
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 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	112 ##
d85b2485af9e doc: grammarcheck m-files ahead of 4.4 release. Rik <rik@octave.org> parents: 24863 diff changeset	113 ## @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	114 ## 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	115 ## 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	116 ##
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	117 ## @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	118 ## @nospell{HPD}.
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	119 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	120 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	121 ## @item
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	122 ## @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	123 ## 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	124 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	125 ## @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	126 ## @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	127 ## 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	128 ## 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	129 ## 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	130 ##
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	131 ## @item
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	132 ## @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	133 ## @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	134 ## @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	135 ## 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	136 ## (@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	137 ## 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	138 ## @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	139 ## @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	140 ## 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	141 ## @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	142 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	143 ## @item
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	144 ## @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	145 ## 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	146 ## @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	147 ## 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	148 ## @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	149 ## @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	150 ## @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	151 ## @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	152 ## 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	153 ## @end itemize
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	154 ##
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 ##
6266e321ef22 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 ## 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	157 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	158 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	159 ## @group
14327 4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	160 ## 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	161 ## 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	162 ## 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	163 ## 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	164 ## 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	165 ## M = M1 * M2;
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	166 ## Afcn = @@(x) A * x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	167 ## Mfcn = @@(x) M \ x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	168 ## M1fcn = @@(x) M1 \ x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	169 ## M2fcn = @@(x) M2 \ x;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	170 ## @end group
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 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	174 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	175 ## @example
17681 f6fded839513 pcg.m: Change success criterion to match reference text (bug #40057) Rik <rik@octave.org> parents: 17336 diff changeset	176 ## x = pcg (A, b)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	177 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	178 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	179 ## @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	180 ## @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	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
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	183 ## x = pcg (Afcn, 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	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 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	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 (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	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 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	193 ##
6266e321ef22 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 ## @example
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	195 ## x = pcg (Afcn, b, 1e-6, 100, Mfcn)
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	196 ## @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	197 ##
6266e321ef22 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 ## @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	199 ## 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	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 (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	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 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	206 ##
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	207 ## @example
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	208 ## x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn)
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 ## @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	210 ##
6266e321ef22 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 ## @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	212 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	213 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	214 ## @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	215 ## 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	216 ## 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	217 ## 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	218 ## 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	219 ## 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	220 ## endfunction
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	221 ## Apfcn = @@(x, p) Ap (A, x, p);
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	222 ## x = pcg (Apfcn, b, [], [], [], [], [], 2);
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	223 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	224 ## @end example
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	225 ##
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	226 ## @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	227 ## split preconditioner
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	228 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	229 ## @example
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	230 ## @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	231 ## 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	232 ## 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	233 ## M = M1 * M2;
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	234 ##
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	235 ## ## 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	236 ## [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	237 ##
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	238 ## ## 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	239 ## [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	240 ## 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	241 ##
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	242 ## @end group
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	243 ## @end example
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	244 ##
10791 3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	245 ## References:
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	246 ##
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	247 ## @enumerate
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	248 ## @item
be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	249 ## 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	250 ## 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	251 ##
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10791 diff changeset	252 ## @item
19596 0e1f5a750d00 maint: Periodic merge of gui-release to default. John W. Eaton <jwe@octave.org> parents: 19148 19593 diff changeset	253 ## @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	254 ## @nospell{PWS} 1996. (condition number estimate from PCG)
a4969508008e doc: Periodic spellcheck of the documentation. Rik <rik@octave.org> parents: 14872 diff changeset	255 ## 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	256 ## @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	257 ## @end enumerate
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	258 ##
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	259 ## @seealso{sparse, pcr, gmres, bicg, bicgstab, cgs}
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	260 ## @end deftypefn
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	261
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	262 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	263 pcg (A, b, tol = [], maxit = [], M1 = [], M2 = [], x0 = [], varargin)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	264
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	265 ## 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	266 [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	267 zeros(size (b))}, tol, maxit, x0);
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	268
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	269 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	270 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	271 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	272 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	273 endif
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	274
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	275 ## 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	276 ## 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	277
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	278 [Afcn, M1fcn, M2fcn] = __alltohandles__ (A, b, M1, M2, "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	279
6266e321ef22 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 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	281 n_arg_out = nargout;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	282
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	283 ## 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	284 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	285 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	286 if (n_arg_out < 2)
28912 0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir. Rik <rik@octave.org> parents: 27985 diff changeset	287 printf ("The right hand side vector is all zero so pcg \n");
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	288 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	289 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	290 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	291 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	292 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	293 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	294 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	295 eigest = [NaN, NaN];
28917 72d57dbcc305 maint: Add semicolon after break and return keywords. Rik <rik@octave.org> parents: 28912 diff changeset	296 return;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	297 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	298
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	299 x = x_pr = x_min = x0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	300
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	301 ## 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	302 ## 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	303
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	304 r = b - feval (Afcn, x, varargin{:});
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	305 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	306 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	307 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	308 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	309 resvec(1, 1) = norm (r);
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	310 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	311 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	312
6266e321ef22 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 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	314 T = zeros (maxit, maxit);
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7515 diff changeset	315 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	316 T = [];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	317 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	318
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	319 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	320 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	321 try
28928 ae7ce8358953 maint: Add semicolon to end of all warning() and error() invocations. Rik <rik@octave.org> parents: 28917 diff changeset	322 warning ("error","Octave:singular-matrix","local");
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	323 z = feval (M1fcn, r, varargin{:});
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	324 z = feval (M2fcn, z, varargin{:});
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	325 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	326 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	327 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	328 end_try_catch
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	329 else
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	330 z = feval (M1fcn, r, varargin{:});
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	331 z = feval (M2fcn, z, varargin{:});
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	332 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	333
11587 c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11563 diff changeset	334 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	335 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	336 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	337 old_tau = tau;
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	338 p = z + beta * p;
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	339 w = feval (Afcn, p, varargin{:});
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	340
6266e321ef22 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 ## 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	342
6266e321ef22 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 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	344 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	345 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	346
6266e321ef22 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 ## 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	348 ## 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	349 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	350 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	351 (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	352 (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	353 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	354 break;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	355 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	356
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	357 x += alpha * p;
392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	358 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	359 resvec(iter, 1) = norm (r);
27356 4e632e942e84 doc: Improve documentation of sparse functions. Kai T. Ohlhus <k.ohlhus@gmail.com> parents: 26376 diff changeset	360 ## Check if the iterated has minimum residual
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	361 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	362 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	363 iter_min = iter - 1;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	364 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	365 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	366 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	367 [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	368 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	369 endif
20735 418ae0cb752f Replace ++,-- with in-place operators for performance. Rik <rik@octave.org> parents: 20727 diff changeset	370 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	371 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	372 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	373 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	374 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	375 x_pr = x;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	376 endwhile
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	377
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	378 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	379 ## 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	380 ## residual.
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	381 z = feval (M1fcn, r, varargin{:});
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	382 z = feval (M2fcn, z, varargin{:});
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	383 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	384
6266e321ef22 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 ## (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	386 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	387 if (flag != 4)
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	388 if (iter > 3)
10549 95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	389 T = T(2:iter-2,2:iter-2);
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	390 l = eig (T);
95c3e38098bf Untabify .m scripts Rik <code@nomad.inbox5.com> parents: 9051 diff changeset	391 eigest = [min(l), max(l)];
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	392 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	393 eigest = [NaN, NaN];
24849 96c74e33d17a pcg.m: Overhaul function and BIST tests (bug #53258). Rik <rik@octave.org> parents: 24534 diff changeset	394 warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	395 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	396 else
5838 376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe] jwe parents: 5837 diff changeset	397 eigest = [NaN, NaN];
28928 ae7ce8358953 maint: Add semicolon to end of all warning() and error() invocations. Rik <rik@octave.org> parents: 28917 diff changeset	398 warning ('pcg: eigenvalue estimate failed: matrix not positive definite?');
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	399 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	400 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	401 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	402 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	403 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	404 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	405 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	406
6266e321ef22 Merge code from SoCiS 2016 "Improve iterative methods" project (http://socis16octave-improveiterativemethods.blogspot.com) Cristiano Dorigo <cristiano.dorigo@hotmail.it> parents: 24849 diff changeset	407 ## Set the last variables
6747 392d61107f11 [project @ 2007-06-19 20:30:30 by dbateman] dbateman parents: 6547 diff changeset	408
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	409 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	410 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	411 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	412 relres = 0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	413 else
29221 56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	414 relres = resvec(iter_min+1, 1) ./ b_norm;
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	415 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	416
6266e321ef22 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 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	418
6266e321ef22 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 ## 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	420 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	421 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	422 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	423 flag = 0;
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	424 endif
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	425
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	426 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	427 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	428 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	429 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	430 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	431 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	432 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	433 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	434 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	435 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	436 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	437 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	438 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	439 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	440 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	441 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	442 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	443 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	444 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	445 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	446 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	447 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	448 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	449 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	450 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	451 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	452 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	453 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	454 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	455 endswitch
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	456 endif
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28928 diff changeset	457
5837 55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	458 endfunction
55404f3b0da1 [project @ 2006-06-01 19:05:31 by jwe] jwe parents: diff changeset	459
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	460
6266e321ef22 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 %!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	462 %! 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	463 %! 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	464 %! b = A * ones (n, 1);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	465 %! 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	466 %! 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	467 %! 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	468 %! x = pcg (A, b);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	469 %! Afcn = @(x) A * x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	470 %! x = pcg (Afcn, 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	471 %! 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	472 %! x = pcg (A, b, 1e-6, 100, M1, M2);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	473 %! Mfcn = @(x) M \ x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	474 %! x = pcg (Afcn, b, 1e-6, 100, Mfcn);
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	475 %! M1fcn = @(x) M1 \ x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	476 %! M2fcn = @(x) M2 \ x;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	477 %! x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	478 %! 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	479 %! 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	480 %! 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	481 %! 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	482 %! 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	483 %! endfunction
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	484 %! Afcn = @(x, p) Ap (A, x, p);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	485 %! ## solution of A^2 * x = b
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	486 %! x = pcg (Afcn, 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	487
6266e321ef22 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 %!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	489 %! 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	490 %! 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	491 %! b = A * ones (n, 1);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	492 %! 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	493 %! 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	494 %! 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	495 %!
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	496 %! ## 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	497 %! [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	498 %! 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	499 %!
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	500 %! ## 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	501 %! [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	502 %! x = M2 \ y # compare x and x_ref
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	503 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	504 %! ## Check that all type of inputs work
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	505 %! A = toeplitz (sparse ([2, 1 ,0, 0, 0]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	506 %! b = A * ones (5, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	507 %! M1 = diag (sqrt (diag (A)));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	508 %! M2 = M1; # M1 * M2 is the Jacobi preconditioner
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	509 %! Afcn = @(z) A*z;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	510 %! M1_fcn = @(z) M1 \ z;
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	511 %! M2_fcn = @(z) M2 \ z;
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	512 %! [x, flag, ~, iter] = pcg (A,b);
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, ~ , iter] = pcg (A, b, [], [], M1, M2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	517 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	518 %! [x, flag] = pcg (A, b, [], [], M1_fcn, M2_fcn);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	519 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	520 %! [x, flag] = pcg (A, b,[],[], M1_fcn, M2);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	521 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	522 %! [x, flag] = pcg (A, b,[],[], M1, M2_fcn);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	523 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	524 %! [x, flag] = pcg (Afcn, b);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	525 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	526 %! [x, flag] = pcg (Afcn, b,[],[], M1 * M2);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	527 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	528 %! [x, flag] = pcg (Afcn, b,[],[], M1, M2);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	529 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	530 %! [x, flag] = pcg (Afcn, b,[],[], M1_fcn, M2);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	531 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	532 %! [x, flag] = pcg (Afcn, b,[],[], M1, M2_fcn);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	533 %! assert (flag, 0);
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	534 %! [x, flag] = pcg (Afcn, b,[],[], M1_fcn, M2_fcn);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	535 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	536
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	537 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	538 %! ## solve a small diagonal system
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	539 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	540 %! A = diag ([1:N]); b = rand (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	541 %! [x, flag] = pcg (A, b, [], N+1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	542 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	543 %! 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	544
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	545 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	546 %! ## A is not positive definite
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	547 %! ## The indefiniteness of A is detected.
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	548 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	549 %! A = -diag ([1:N]); b = sum (A, 2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	550 %! [x, flag] = pcg (A, b, [], N + 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	551 %! assert (flag, 4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	552
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	553 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	554 %! ## 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	555 %! N = 100;
25755 ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	556 %! ## Form 1-D Laplacian matrix
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	557 %! A = 2 * eye (N,N);
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	558 %! A(2:(N+1):end) = -1;
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	559 %! A((N+1):(N+1):end) = -1;
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	560 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	561 %! [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	562 %! assert (flag);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	563 %! assert (relres >= 1.0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	564
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	565 %!warning <iteration converged too fast>
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	566 %! ## solve tridiagonal system with "perfect" preconditioner which converges
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	567 %! ## 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	568 %! N = 100;
25755 ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	569 %! ## Form 1-D Laplacian matrix
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	570 %! A = 2 * eye (N,N);
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	571 %! A(2:(N+1):end) = -1;
ddf8a03beffb Calculate Laplacian matrix in BIST tests without for loop. Rik <rik@octave.org> parents: 25436 diff changeset	572 %! A((N+1):(N+1):end) = -1;
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	573 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	574 %! [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	575 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	576 %! 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	577 %!
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	578 %! assert (isnan (eigest), isnan ([NaN, NaN]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	579
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	580 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	581 %! ## 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	582 %! ## 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	583 %! ## makes iterations until it reaches maxit.
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	584 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	585 %! A = diag (1:N) + 1e-4*i;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	586 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	587 %! [x, flag] = pcg (A, b, []);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	588 %! assert (flag, 4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	589
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	590 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	591 %! ## 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	592 %! N = 10;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	593 %! A = diag (1:N) + 1e-10*i;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	594 %! b = ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	595 %! [x, flag] = pcg (A, b, [], N+1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	596 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	597 %! 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	598
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	599 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	600 %! ## pcg solves linear system with A Hermitian positive definite
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	601 %! N = 20;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	602 %! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ...
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	603 %! 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	604 %! b = A * ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	605 %! Hermitian_A = ishermitian (A);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	606 %! [x, flag] = pcg (A, b, [], 2*N);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	607 %! assert (Hermitian_A, true);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	608 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	609 %! assert (x, ones (N, 1), -1e-4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	610
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	611 %!testif HAVE_CHOLMOD
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	612 %! ## pcg solves preconditioned linear system with A HPD
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	613 %! N = 20;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	614 %! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ...
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	615 %! 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	616 %! b = A * ones (N, 1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	617 %! 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	618 %! M = M2' * M2;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	619 %! Hermitian_A = ishermitian (A);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	620 %! Hermitian_M = ishermitian (M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	621 %! [x, flag] = pcg (A, b, [], 2*N, M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	622 %! assert (Hermitian_A, true);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	623 %! assert (Hermitian_M, true);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	624 %! assert (flag, 0);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	625 %! assert (x, ones (N, 1), -1e-4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	626
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	627 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	628 %! ## pcg recognizes that the preconditioner matrix is singular
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	629 %! N = 3;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	630 %! A = toeplitz ([2, 1, 0]);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	631 %! 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	632 %! [x, flag] = pcg (A, ones (3, 1), [], [], M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	633 %! assert (flag, 2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	634
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	635 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	636 %! A = rand (4);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	637 %! A = A' * A;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	638 %! [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	639 %! assert (x, zeros (4, 1));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	640
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	641 ## Test return types
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	642 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	643 %! A = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	644 %! b = 1;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	645 %! [x, flag] = pcg (A, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	646 %! assert (class (x), "single");
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	647
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	648 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	649 %! A = 1;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	650 %! b = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	651 %! [x, flag] = pcg (A, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	652 %! assert (class (x), "single");
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	653
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	654 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	655 %! A = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	656 %! b = single (1);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	657 %! [x, flag] = pcg (A, b);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	658 %! assert (class (x), "single");
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	659
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	660 %!test
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	661 %!function y = Afcn (x)
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	662 %! A = toeplitz ([2, 1, 0, 0]);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	663 %! y = A * x;
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	664 %!endfunction
30893 e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files. Rik <rik@octave.org> parents: 30564 diff changeset	665 %! [x, flag] = pcg ("Afcn", [3; 4; 4; 3]);
25148 dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	666 %! assert (x, ones (4, 1), 1e-6);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	667
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	668 %!test
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	669 %! ## unpreconditioned residual
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	670 %! A = toeplitz (sparse ([4, 1, 0, 0, 0]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	671 %! b = sum (A, 2);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	672 %! M = toeplitz (sparse ([2, 1, 0, 0, 0]));
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	673 %! [x, flag, relres] = pcg (A, b, [], 2, M);
dedc0128645a pcg.m: Clean up BIST tests. Rik <rik@octave.org> parents: 25143 diff changeset	674 %! 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	675
29221 56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	676 %!test <*59776>
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	677 %! A = [ 1.00000000 -0.00054274 -0.00066848;
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	678 %! -0.00054274 1.00000000 -0.00060330;
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	679 %! -0.00066848 -0.00060330 1.00000000];
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	680 %! b = [1 1 1]';
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	681 %! [x, flag, relres, iter, resvec] = pcg (A, b, 1e-6, 4, [], [], [1; 1; 1]);
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	682 %! assert (flag, 0);
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	683 %! assert (relres, resvec(2) / norm (b));
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	684 %! assert (iter, 1);
56978d4b266c pcg.m: Return correct FLAG and correct RELRES output (bug #59776) Rik <rik@octave.org> parents: 27985 diff changeset	685