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