Mercurial > octave-antonio
annotate scripts/sparse/pcg.m @ 20164:df437a52bcaf stable
doc: Update more docstrings to have one sentence summary as first line.
Reviewed miscellaneous, sparse, strings in scripts directory.
* scripts/miscellaneous/bzip2.m, scripts/miscellaneous/citation.m,
scripts/miscellaneous/compare_versions.m, scripts/miscellaneous/computer.m,
scripts/miscellaneous/debug.m, scripts/miscellaneous/dir.m,
scripts/miscellaneous/edit.m, scripts/miscellaneous/error_ids.m,
scripts/miscellaneous/fileattrib.m, scripts/miscellaneous/fullfile.m,
scripts/miscellaneous/genvarname.m, scripts/miscellaneous/gzip.m,
scripts/miscellaneous/mkoctfile.m, scripts/miscellaneous/news.m,
scripts/miscellaneous/open.m, scripts/miscellaneous/parseparams.m,
scripts/miscellaneous/recycle.m, scripts/miscellaneous/run.m,
scripts/miscellaneous/swapbytes.m, scripts/miscellaneous/tar.m,
scripts/miscellaneous/tmpnam.m, scripts/miscellaneous/unpack.m,
scripts/miscellaneous/what.m, scripts/sparse/bicg.m, scripts/sparse/bicgstab.m,
scripts/sparse/cgs.m, scripts/sparse/colperm.m, scripts/sparse/eigs.m,
scripts/sparse/etreeplot.m, scripts/sparse/gmres.m, scripts/sparse/gplot.m,
scripts/sparse/ichol.m, scripts/sparse/ilu.m, scripts/sparse/pcg.m,
scripts/sparse/pcr.m, scripts/sparse/qmr.m, scripts/sparse/spaugment.m,
scripts/sparse/spconvert.m, scripts/sparse/spdiags.m, scripts/sparse/spfun.m,
scripts/sparse/spones.m, scripts/sparse/sprandsym.m, scripts/sparse/spstats.m,
scripts/sparse/spy.m, scripts/sparse/svds.m, scripts/sparse/treelayout.m,
scripts/sparse/treeplot.m, scripts/strings/base2dec.m,
scripts/strings/bin2dec.m, scripts/strings/blanks.m, scripts/strings/cstrcat.m,
scripts/strings/deblank.m, scripts/strings/dec2base.m,
scripts/strings/dec2bin.m, scripts/strings/dec2hex.m,
scripts/strings/findstr.m, scripts/strings/hex2dec.m, scripts/strings/index.m,
scripts/strings/isletter.m, scripts/strings/isstrprop.m,
scripts/strings/mat2str.m, scripts/strings/ostrsplit.m,
scripts/strings/regexptranslate.m, scripts/strings/rindex.m,
scripts/strings/str2num.m, scripts/strings/strcat.m, scripts/strings/strchr.m,
scripts/strings/strjoin.m, scripts/strings/strjust.m,
scripts/strings/strmatch.m, scripts/strings/strsplit.m,
scripts/strings/strtok.m, scripts/strings/strtrim.m,
scripts/strings/strtrunc.m, scripts/strings/substr.m,
scripts/strings/untabify.m, scripts/time/datenum.m:
Update more docstrings to have one sentence summary as first line.
author | Rik <rik@octave.org> |
---|---|
date | Mon, 04 May 2015 14:22:02 -0700 |
parents | 9fc020886ae9 |
children |
rev | line source |
---|---|
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19596
diff
changeset
|
1 ## Copyright (C) 2004-2015 Piotr Krzyzanowski |
5837 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
5837 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
5837 | 18 |
19 ## -*- texinfo -*- | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
20 ## @deftypefn {Function File} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{}) |
5837 | 21 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{}) |
22 ## | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
23 ## 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
|
24 ## 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
|
25 ## |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
26 ## The input arguments are |
5837 | 27 ## |
28 ## @itemize | |
29 ## @item | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
30 ## @var{A} can be either a square (preferably sparse) matrix or a function |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
31 ## handle, inline function or string containing the name of a function which |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
32 ## computes @w{@code{@var{A} * @var{x}}}. In principle, @var{A} should be |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
33 ## symmetric and positive definite; if @code{pcg} finds @var{A} not to be |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
34 ## positive definite, a warning is printed and the @var{flag} output will be |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
35 ## set. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
36 ## |
5837 | 37 ## @item |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
38 ## @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
|
39 ## |
5837 | 40 ## @item |
41 ## @var{tol} is the required relative tolerance for the residual error, | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
42 ## @w{@code{@var{b} - @var{A} * @var{x}}}. The iteration stops if |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
43 ## @w{@code{norm (@var{b} - @var{A} * @var{x})} @leq{} |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
44 ## @w{@var{tol} * norm (@var{b})}}. |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
45 ## 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
|
46 ## |
5837 | 47 ## @item |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
48 ## @var{maxit} is the maximum allowable number of iterations; if @var{maxit} |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
49 ## 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
|
50 ## |
5837 | 51 ## @item |
10846
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
52 ## @var{m} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that |
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
53 ## the iteration is (theoretically) equivalent to solving by @code{pcg} |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
54 ## @w{@code{@var{P} * @var{x} = @var{m} \ @var{b}}}, with |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
55 ## @w{@code{@var{P} = @var{m} \ @var{A}}}. |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
56 ## Note that a proper choice of the preconditioner may dramatically improve |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
57 ## the overall performance of the method. Instead of matrices @var{m1} and |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
58 ## @var{m2}, the user may pass two functions which return the results of |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
59 ## applying the inverse of @var{m1} and @var{m2} to a vector (usually this is |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
60 ## the preferred way of using the preconditioner). If @var{m1} is omitted or |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
61 ## empty @code{[]} then no preconditioning is applied. If @var{m2} is |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
62 ## omitted, @var{m} = @var{m1} will be used as a preconditioner. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
63 ## |
5837 | 64 ## @item |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
65 ## @var{x0} is the initial guess. If @var{x0} is omitted or empty then the |
5837 | 66 ## function sets @var{x0} to a zero vector by default. |
67 ## @end itemize | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
68 ## |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
69 ## The arguments which follow @var{x0} are treated as parameters, and passed in |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
70 ## a proper way to any of the functions (@var{A} or @var{m}) which are passed |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
71 ## to @code{pcg}. See the examples below for further details. The output |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
72 ## arguments are |
5837 | 73 ## |
74 ## @itemize | |
75 ## @item | |
76 ## @var{x} is the computed approximation to the solution of | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
77 ## @w{@code{@var{A} * @var{x} = @var{b}}}. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
78 ## |
5837 | 79 ## @item |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
80 ## @var{flag} reports on the convergence. A value of 0 means the solution |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
81 ## converged and the tolerance criterion given by @var{tol} is satisfied. |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
82 ## A value of 1 means that the @var{maxit} limit for the iteration count was |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
83 ## reached. A value of 3 indicates that the (preconditioned) matrix was found |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
84 ## not to be positive definite. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
85 ## |
5837 | 86 ## @item |
87 ## @var{relres} is the ratio of the final residual to its initial value, | |
88 ## 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
|
89 ## |
5837 | 90 ## @item |
91 ## @var{iter} is the actual number of iterations performed. | |
92 ## | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
93 ## @item |
5837 | 94 ## @var{resvec} describes the convergence history of the method. |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
95 ## @code{@var{resvec}(i,1)} is the Euclidean norm of the residual, and |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
96 ## @code{@var{resvec}(i,2)} is the preconditioned residual norm, after the |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
97 ## (@var{i}-1)-th iteration, @code{@var{i} = 1, 2, @dots{}, @var{iter}+1}. |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
98 ## The preconditioned residual norm is defined as |
5838 | 99 ## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
100 ## @code{@var{r} = @var{b} - @var{A} * @var{x}}, see also the |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
101 ## description of @var{m}. If @var{eigest} is not required, only |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
102 ## @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
|
103 ## |
5837 | 104 ## @item |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
105 ## @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
|
106 ## and largest @code{@var{eigest}(2)} eigenvalues of the preconditioned matrix |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
107 ## @w{@code{@var{P} = @var{m} \ @var{A}}}. In particular, if no |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
108 ## 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
|
109 ## @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
|
110 ## @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
|
111 ## @code{@var{eigest}(2) / @var{eigest}(1)} is a lower bound for |
19593
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
17744
diff
changeset
|
112 ## @code{cond (@var{P}, 2)}, which nevertheless in the limit should |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
113 ## theoretically be equal to the actual value of the condition number. |
5837 | 114 ## The method which computes @var{eigest} works only for symmetric positive |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
115 ## definite @var{A} and @var{m}, and the user is responsible for verifying this |
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
116 ## assumption. |
5837 | 117 ## @end itemize |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
118 ## |
5837 | 119 ## Let us consider a trivial problem with a diagonal matrix (we exploit the |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
120 ## sparsity of A) |
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
121 ## |
5837 | 122 ## @example |
123 ## @group | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
124 ## n = 10; |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
125 ## A = diag (sparse (1:n)); |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
126 ## b = rand (n, 1); |
19148 | 127 ## [l, u, p] = ilu (A, struct ("droptol", 1.e-3)); |
5837 | 128 ## @end group |
129 ## @end example | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
130 ## |
5837 | 131 ## @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
|
132 ## |
5837 | 133 ## @example |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
134 ## x = pcg (A, b) |
5837 | 135 ## @end example |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
136 ## |
5837 | 137 ## @sc{Example 2:} @code{pcg} with a function which computes |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
138 ## @code{@var{A} * @var{x}} |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
139 ## |
5837 | 140 ## @example |
141 ## @group | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
142 ## function y = apply_a (x) |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
143 ## y = [1:N]' .* x; |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
144 ## endfunction |
5837 | 145 ## |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
146 ## x = pcg ("apply_a", b) |
5837 | 147 ## @end group |
148 ## @end example | |
6747 | 149 ## |
150 ## @sc{Example 3:} @code{pcg} with a preconditioner: @var{l} * @var{u} | |
151 ## | |
152 ## @example | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
153 ## x = pcg (A, b, 1.e-6, 500, l*u) |
6747 | 154 ## @end example |
155 ## | |
156 ## @sc{Example 4:} @code{pcg} with a preconditioner: @var{l} * @var{u}. | |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
157 ## Faster than @sc{Example 3} since lower and upper triangular matrices are |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
158 ## easier to invert |
6747 | 159 ## |
160 ## @example | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
161 ## x = pcg (A, b, 1.e-6, 500, l, u) |
6747 | 162 ## @end example |
163 ## | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
164 ## @sc{Example 5:} Preconditioned iteration, with full diagnostics. The |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
165 ## preconditioner (quite strange, because even the original matrix @var{A} is |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
166 ## trivial) is defined as a function |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
167 ## |
5837 | 168 ## @example |
169 ## @group | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
170 ## function y = apply_m (x) |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
171 ## k = floor (length (x) - 2); |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
172 ## y = x; |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
173 ## y(1:k) = x(1:k) ./ [1:k]'; |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
174 ## endfunction |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
175 ## |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
176 ## [x, flag, relres, iter, resvec, eigest] = ... |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
177 ## pcg (A, b, [], [], "apply_m"); |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
178 ## semilogy (1:iter+1, resvec); |
5837 | 179 ## @end group |
180 ## @end example | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
181 ## |
20164
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
182 ## @sc{Example 6:} Finally, a preconditioner which depends on a parameter |
df437a52bcaf
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
183 ## @var{k}. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
184 ## |
5837 | 185 ## @example |
186 ## @group | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
187 ## function y = apply_M (x, varargin) |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
188 ## K = varargin@{1@}; |
6547 | 189 ## y = x; |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
190 ## y(1:K) = x(1:K) ./ [1:K]'; |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
191 ## endfunction |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
192 ## |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
193 ## [x, flag, relres, iter, resvec, eigest] = ... |
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
194 ## pcg (A, b, [], [], "apply_m", [], [], 3) |
5837 | 195 ## @end group |
196 ## @end example | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
197 ## |
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
198 ## References: |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
199 ## |
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
200 ## @enumerate |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
201 ## @item |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
202 ## 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
|
203 ## 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
|
204 ## |
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
205 ## @item |
19596
0e1f5a750d00
maint: Periodic merge of gui-release to default.
John W. Eaton <jwe@octave.org>
diff
changeset
|
206 ## @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
|
207 ## @nospell{PWS} 1996. (condition number estimate from PCG) |
a4969508008e
doc: Periodic spellcheck of the documentation.
Rik <rik@octave.org>
parents:
14872
diff
changeset
|
208 ## Revised version of this book is available online at |
a4969508008e
doc: Periodic spellcheck of the documentation.
Rik <rik@octave.org>
parents:
14872
diff
changeset
|
209 ## @url{http://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
|
210 ## @end enumerate |
5837 | 211 ## |
212 ## @seealso{sparse, pcr} | |
213 ## @end deftypefn | |
214 | |
5838 | 215 ## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl> |
6747 | 216 ## Modified by: Vittoria Rezzonico <vittoria.rezzonico@epfl.ch> |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
217 ## - Add the ability to provide the pre-conditioner as two separate matrices |
5837 | 218 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
219 function [x, flag, relres, iter, resvec, eigest] = pcg (A, b, tol, maxit, m1, m2, x0, varargin) |
5837 | 220 |
8507 | 221 ## M = M1*M2 |
6747 | 222 |
223 if (nargin < 7 || isempty (x0)) | |
5838 | 224 x = zeros (size (b)); |
5837 | 225 else |
226 x = x0; | |
227 endif | |
228 | |
8507 | 229 if (nargin < 5 || isempty (m1)) |
230 exist_m1 = 0; | |
231 else | |
232 exist_m1 = 1; | |
233 endif | |
6747 | 234 |
8507 | 235 if (nargin < 6 || isempty (m2)) |
236 exist_m2 = 0; | |
237 else | |
238 exist_m2 = 1; | |
239 endif | |
5837 | 240 |
5838 | 241 if (nargin < 4 || isempty (maxit)) |
14872
c2dbdeaa25df
maint: use rows() and columns() to clarify m-files.
Rik <octave@nomad.inbox5.com>
parents:
14868
diff
changeset
|
242 maxit = min (rows (b), 20); |
5837 | 243 endif |
244 | |
5838 | 245 maxit += 2; |
5837 | 246 |
5838 | 247 if (nargin < 3 || isempty (tol)) |
5837 | 248 tol = 1e-6; |
249 endif | |
250 | |
251 preconditioned_residual_out = false; | |
252 if (nargout > 5) | |
5838 | 253 T = zeros (maxit, maxit); |
5837 | 254 preconditioned_residual_out = true; |
255 endif | |
256 | |
8506 | 257 ## Assume A is positive definite. |
258 matrix_positive_definite = true; | |
5837 | 259 |
5838 | 260 p = zeros (size (b)); |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
261 oldtau = 1; |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
262 if (isnumeric (A)) |
8506 | 263 ## A is a matrix. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
264 r = b - A*x; |
8506 | 265 else |
266 ## A should be a function. | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
267 r = b - feval (A, x, varargin{:}); |
5837 | 268 endif |
269 | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
270 b_norm = norm (b); |
5838 | 271 resvec(1,1) = norm (r); |
5837 | 272 alpha = 1; |
273 iter = 2; | |
274 | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
275 while (resvec(iter-1,1) > tol * b_norm && iter < maxit) |
8507 | 276 if (exist_m1) |
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
277 if (isnumeric (m1)) |
10549 | 278 y = m1 \ r; |
6747 | 279 else |
10549 | 280 y = feval (m1, r, varargin{:}); |
5837 | 281 endif |
6747 | 282 else |
283 y = r; | |
284 endif | |
8507 | 285 if (exist_m2) |
286 if (isnumeric (m2)) | |
10549 | 287 z = m2 \ y; |
6747 | 288 else |
10549 | 289 z = feval (m2, y, varargin{:}); |
6747 | 290 endif |
291 else | |
292 z = y; | |
5837 | 293 endif |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11563
diff
changeset
|
294 tau = z' * r; |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
295 resvec(iter-1,2) = sqrt (tau); |
5838 | 296 beta = tau / oldtau; |
5837 | 297 oldtau = tau; |
6747 | 298 p = z + beta * p; |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
299 if (isnumeric (A)) |
8506 | 300 ## A is a matrix. |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
301 w = A * p; |
8506 | 302 else |
303 ## A should be a function. | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10846
diff
changeset
|
304 w = feval (A, p, varargin{:}); |
5837 | 305 endif |
8506 | 306 ## Needed only for eigest. |
307 oldalpha = alpha; | |
5838 | 308 alpha = tau / (p'*w); |
8506 | 309 if (alpha <= 0.0) |
310 ## Negative matrix. | |
5837 | 311 matrix_positive_definite = false; |
312 endif | |
6747 | 313 x += alpha * p; |
314 r -= alpha * w; | |
5838 | 315 if (nargout > 5 && iter > 2) |
5837 | 316 T(iter-1:iter, iter-1:iter) = T(iter-1:iter, iter-1:iter) + ... |
10549 | 317 [1 sqrt(beta); sqrt(beta) beta]./oldalpha; |
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
318 ## EVS = eig (T(2:iter-1,2:iter-1)); |
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
319 ## fprintf (stderr,"PCG condest: %g (iteration: %d)\n", max (EVS)/min (EVS),iter); |
5837 | 320 endif |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
321 resvec(iter,1) = norm (r); |
5838 | 322 iter++; |
5837 | 323 endwhile |
324 | |
325 if (nargout > 5) | |
5838 | 326 if (matrix_positive_definite) |
5837 | 327 if (iter > 3) |
10549 | 328 T = T(2:iter-2,2:iter-2); |
329 l = eig (T); | |
330 eigest = [min(l), max(l)]; | |
331 ## fprintf (stderr, "pcg condest: %g\n", eigest(2)/eigest(1)); | |
5837 | 332 else |
10549 | 333 eigest = [NaN, NaN]; |
11588
d5bd2766c640
style fixes for warning and error messages in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
334 warning ("pcg: eigenvalue estimate failed: iteration converged too fast"); |
5837 | 335 endif |
336 else | |
5838 | 337 eigest = [NaN, NaN]; |
5837 | 338 endif |
339 | |
8506 | 340 ## Apply the preconditioner once more and finish with the precond |
341 ## residual. | |
8507 | 342 if (exist_m1) |
343 if (isnumeric (m1)) | |
10549 | 344 y = m1 \ r; |
6747 | 345 else |
10549 | 346 y = feval (m1, r, varargin{:}); |
5837 | 347 endif |
6747 | 348 else |
349 y = r; | |
5837 | 350 endif |
8507 | 351 if (exist_m2) |
352 if (isnumeric (m2)) | |
10549 | 353 z = m2 \ y; |
6747 | 354 else |
10549 | 355 z = feval (m2, y, varargin{:}); |
6747 | 356 endif |
357 else | |
358 z = y; | |
359 endif | |
360 | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
361 resvec(iter-1,2) = sqrt (r' * z); |
5837 | 362 else |
5838 | 363 resvec = resvec(:,1); |
5837 | 364 endif |
365 | |
366 flag = 0; | |
17681
f6fded839513
pcg.m: Change success criterion to match reference text (bug #40057)
Rik <rik@octave.org>
parents:
17336
diff
changeset
|
367 relres = resvec(iter-1,1) ./ resvec(1,1); |
5838 | 368 iter -= 2; |
6747 | 369 if (iter >= maxit - 2) |
5837 | 370 flag = 1; |
371 if (nargout < 2) | |
6498 | 372 warning ("pcg: maximum number of iterations (%d) reached\n", iter); |
6747 | 373 warning ("the initial residual norm was reduced %g times.\n", ... |
10549 | 374 1.0 / relres); |
5837 | 375 endif |
5838 | 376 elseif (nargout < 2) |
6498 | 377 fprintf (stderr, "pcg: converged in %d iterations. ", iter); |
378 fprintf (stderr, "the initial residual norm was reduced %g times.\n",... | |
10549 | 379 1.0/relres); |
5837 | 380 endif |
381 | |
5838 | 382 if (! matrix_positive_definite) |
5837 | 383 flag = 3; |
384 if (nargout < 2) | |
6498 | 385 warning ("pcg: matrix not positive definite?\n"); |
5837 | 386 endif |
387 endif | |
388 endfunction | |
389 | |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
390 |
5837 | 391 %!demo |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
392 %! ## Simplest usage of pcg (see also 'help pcg') |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
393 %! |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
394 %! N = 10; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
395 %! A = diag ([1:N]); b = rand (N, 1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
396 %! y = A \ b; # y is the true solution |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
397 %! x = pcg (A, b); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
398 %! printf ("The solution relative error is %g\n", norm (x - y) / norm (y)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
399 %! |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
400 %! ## You shouldn't be afraid if pcg issues some warning messages in this |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
401 %! ## example: watch out in the second example, why it takes N iterations |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
402 %! ## of pcg to converge to (a very accurate, by the way) solution |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
403 |
5837 | 404 %!demo |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
405 %! ## Full output from pcg, except for the eigenvalue estimates |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
406 %! ## We use this output to plot the convergence history |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
407 %! |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
408 %! N = 10; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
409 %! A = diag ([1:N]); b = rand (N, 1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
410 %! X = A \ b; # X is the true solution |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
411 %! [x, flag, relres, iter, resvec] = pcg (A, b); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
412 %! printf ("The solution relative error is %g\n", norm (x - X) / norm (X)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
413 %! title ("Convergence history"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
414 %! semilogy ([0:iter], resvec / resvec(1), "o-g"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
415 %! xlabel ("Iteration"); ylabel ("log(||b-Ax||/||b||)"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
416 %! legend ("relative residual"); |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
417 |
5837 | 418 %!demo |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
419 %! ## Full output from pcg, including the eigenvalue estimates |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
420 %! ## Hilbert matrix is extremely ill-conditioned, so pcg WILL have problems |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
421 %! |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
422 %! N = 10; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
423 %! A = hilb (N); b = rand (N, 1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
424 %! X = A \ b; # X is the true solution |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
425 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
426 %! printf ("The solution relative error is %g\n", norm (x - X) / norm (X)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
427 %! printf ("Condition number estimate is %g\n", eigest(2) / eigest(1)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
428 %! printf ("Actual condition number is %g\n", cond (A)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
429 %! title ("Convergence history"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
430 %! semilogy ([0:iter], resvec, ["o-g";"+-r"]); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
431 %! xlabel ("Iteration"); ylabel ("log(||b-Ax||)"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
432 %! legend ("absolute residual", "absolute preconditioned residual"); |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
433 |
5837 | 434 %!demo |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
435 %! ## Full output from pcg, including the eigenvalue estimates |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
436 %! ## We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2) |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
437 %! ## and that's the reason we need some preconditioner; here we take |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
438 %! ## a very simple and not powerful Jacobi preconditioner, |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
439 %! ## which is the diagonal of A. |
5837 | 440 %! |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
441 %! N = 100; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
442 %! A = zeros (N, N); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
443 %! for i = 1 : N - 1 # form 1-D Laplacian matrix |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
444 %! A(i:i+1, i:i+1) = [2 -1; -1 2]; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
445 %! endfor |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
446 %! b = rand (N, 1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
447 %! X = A \ b; # X is the true solution |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
448 %! maxit = 80; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
449 %! printf ("System condition number is %g\n", cond (A)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
450 %! ## No preconditioner: the convergence is very slow! |
5837 | 451 %! |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
452 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
453 %! printf ("System condition number estimate is %g\n", eigest(2) / eigest(1)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
454 %! title ("Convergence history"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
455 %! semilogy ([0:iter], resvec(:,1), "o-g"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
456 %! xlabel ("Iteration"); ylabel ("log(||b-Ax||)"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
457 %! legend ("NO preconditioning: absolute residual"); |
5837 | 458 %! |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
459 %! pause (1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
460 %! ## Test Jacobi preconditioner: it will not help much!!! |
5837 | 461 %! |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
462 %! M = diag (diag (A)); # Jacobi preconditioner |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
463 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
464 %! printf ("JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
465 %! hold on; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
466 %! semilogy ([0:iter], resvec(:,1), "o-r"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
467 %! legend ("NO preconditioning: absolute residual", ... |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
468 %! "JACOBI preconditioner: absolute residual"); |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
469 %! |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
470 %! pause (1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
471 %! ## Test nonoverlapping block Jacobi preconditioner: it will help much! |
5837 | 472 %! |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
473 %! M = zeros (N, N); k = 4; |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
474 %! for i = 1 : k : N # form 1-D Laplacian matrix |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
475 %! M(i:i+k-1, i:i+k-1) = A(i:i+k-1, i:i+k-1); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
476 %! endfor |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
477 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
478 %! printf ("BLOCK JACOBI preconditioned system condition number estimate is %g\n", eigest(2) / eigest(1)); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
479 %! semilogy ([0:iter], resvec(:,1), "o-b"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
480 %! legend ("NO preconditioning: absolute residual", ... |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
481 %! "JACOBI preconditioner: absolute residual", ... |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
482 %! "BLOCK JACOBI preconditioner: absolute residual"); |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
483 %! hold off; |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
484 |
5837 | 485 %!test |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
486 %! ## solve small diagonal system |
5837 | 487 %! |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
488 %! N = 10; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
489 %! A = diag ([1:N]); b = rand (N, 1); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
490 %! X = A \ b; # X is the true solution |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
491 %! [x, flag] = pcg (A, b, [], N+1); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
492 %! assert (norm (x - X) / norm (X), 0, 1e-10); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
493 %! assert (flag, 0); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
494 |
5837 | 495 %!test |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
496 %! ## solve small indefinite diagonal system |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
497 %! ## despite A is indefinite, the iteration continues and converges |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
498 %! ## indefiniteness of A is detected |
5837 | 499 %! |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
500 %! N = 10; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
501 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
502 %! X = A \ b; # X is the true solution |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
503 %! [x, flag] = pcg (A, b, [], N+1); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
504 %! assert (norm (x - X) / norm (X), 0, 1e-10); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
505 %! assert (flag, 3); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
506 |
5837 | 507 %!test |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
508 %! ## solve tridiagonal system, do not converge in default 20 iterations |
5837 | 509 %! |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
510 %! N = 100; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
511 %! A = zeros (N, N); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
512 %! for i = 1 : N - 1 # form 1-D Laplacian matrix |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
513 %! A(i:i+1, i:i+1) = [2 -1; -1 2]; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
514 %! endfor |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
515 %! b = ones (N, 1); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
516 %! X = A \ b; # X is the true solution |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
517 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
518 %! assert (flag); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
519 %! assert (relres > 1.0); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
520 %! assert (iter, 20); # should perform max allowable default number of iterations |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
521 |
5837 | 522 %!test |
17336
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
523 %! ## solve tridiagonal system with 'perfect' preconditioner |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
524 %! ## which converges in one iteration, so the eigest does not |
b81b9d079515
Use '##' for comments which stand alone on a line.
Rik <rik@octave.org>
parents:
16826
diff
changeset
|
525 %! ## work and issues a warning |
5837 | 526 %! |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
527 %! N = 100; |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
528 %! A = zeros (N, N); |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
529 %! for i = 1 : N - 1 # form 1-D Laplacian matrix |
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
530 %! A(i:i+1, i:i+1) = [2 -1; -1 2]; |
14237
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
531 %! endfor |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
532 %! b = ones (N, 1); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
533 %! X = A \ b; # X is the true solution |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
534 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
535 %! assert (norm (x - X) / norm (X), 0, 1e-6); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
536 %! assert (flag, 0); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
537 %! assert (iter, 1); # should converge in one iteration |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
538 %! assert (isnan (eigest), isnan ([NaN, NaN])); |
11949c9795a0
Revamp %!demos in m-files to use Octave coding conventions on spacing, etc.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
539 |