Mercurial > hg > octave-lyh
annotate scripts/sparse/pcg.m @ 10793:be55736a0783
Grammarcheck the documentation from m-files.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Sun, 18 Jul 2010 20:35:16 -0700 |
| parents | 3140cb7a05a1 |
| children | a4f482e66b65 |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2004, 2006, 2007, 2008, 2009 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 -*- | |
|
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
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 ## | |
| 5838 | 23 ## Solves the linear system of equations @code{@var{a} * @var{x} = |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## @var{b}} by means of the Preconditioned Conjugate Gradient iterative |
|
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 | |
| 5838 | 29 ## @var{a} can be either a square (preferably sparse) matrix or a |
| 5837 | 30 ## function handle, inline function or string containing the name |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
31 ## of a function which computes @code{@var{a} * @var{x}}. In principle |
| 5838 | 32 ## @var{a} should be symmetric and positive definite; if @code{pcg} |
| 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. |
| 35 ## | |
| 36 ## @item | |
| 37 ## @var{b} is the right hand side vector. | |
| 38 ## | |
| 39 ## @item | |
| 40 ## @var{tol} is the required relative tolerance for the residual error, | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
41 ## @code{@var{b} - @var{a} * @var{x}}. The iteration stops if @code{norm |
| 5838 | 42 ## (@var{b} - @var{a} * @var{x}) <= @var{tol} * norm (@var{b} - @var{a} * |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
43 ## @var{x0})}. If @var{tol} is empty or is omitted, the function sets |
| 5837 | 44 ## @code{@var{tol} = 1e-6} by default. |
| 45 ## | |
| 46 ## @item | |
| 47 ## @var{maxit} is the maximum allowable number of iterations; if | |
| 48 ## @code{[]} is supplied for @code{maxit}, or @code{pcg} has less | |
| 49 ## arguments, a default value equal to 20 is used. | |
| 50 ## | |
| 51 ## @item | |
| 6747 | 52 ## @var{m} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that the iteration is |
| 5837 | 53 ## (theoretically) equivalent to solving by @code{pcg} @code{@var{P} * |
| 5838 | 54 ## @var{x} = @var{m} \ @var{b}}, with @code{@var{P} = @var{m} \ @var{a}}. |
| 5837 | 55 ## 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
|
56 ## improve the overall performance of the method. Instead of matrices |
| 6747 | 57 ## @var{m1} and @var{m2}, the user may pass two functions which return |
| 58 ## the results of applying the inverse of @var{m1} and @var{m2} to | |
| 59 ## a vector (usually this is the preferred way of using the preconditioner). | |
| 60 ## 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
|
61 ## preconditioning is applied. If @var{m2} is omitted, @var{m} = @var{m1} |
| 6747 | 62 ## will be used as preconditioner. |
| 5837 | 63 ## |
| 64 ## @item | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
65 ## @var{x0} is the initial guess. If @var{x0} is empty or omitted, the |
| 5837 | 66 ## function sets @var{x0} to a zero vector by default. |
| 67 ## @end itemize | |
| 68 ## | |
| 69 ## The arguments which follow @var{x0} are treated as parameters, and | |
| 5838 | 70 ## 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
|
71 ## 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
|
72 ## details. The output arguments are |
| 5837 | 73 ## |
| 74 ## @itemize | |
| 75 ## @item | |
| 76 ## @var{x} is the computed approximation to the solution of | |
| 5838 | 77 ## @code{@var{a} * @var{x} = @var{b}}. |
| 5837 | 78 ## |
| 79 ## @item | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
80 ## @var{flag} reports on the convergence. @code{@var{flag} = 0} means |
| 5837 | 81 ## 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
|
82 ## 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
|
83 ## for the iteration count was reached. @code{@var{flag} = 3} reports that |
| 5837 | 84 ## the (preconditioned) matrix was found not positive definite. |
| 85 ## | |
| 86 ## @item | |
| 87 ## @var{relres} is the ratio of the final residual to its initial value, | |
| 88 ## measured in the Euclidean norm. | |
| 89 ## | |
| 90 ## @item | |
| 91 ## @var{iter} is the actual number of iterations performed. | |
| 92 ## | |
| 93 ## @item | |
| 94 ## @var{resvec} describes the convergence history of the method. | |
| 95 ## @code{@var{resvec} (i,1)} is the Euclidean norm of the residual, and | |
| 96 ## @code{@var{resvec} (i,2)} is the preconditioned residual norm, | |
| 97 ## 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
|
98 ## 1, 2, @dots{}, @var{iter}+1}. The preconditioned residual norm |
| 6547 | 99 ## is defined as |
| 5838 | 100 ## @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where |
| 101 ## @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
|
102 ## description of @var{m}. If @var{eigest} is not required, only |
| 5837 | 103 ## @code{@var{resvec} (:,1)} is returned. |
| 104 ## | |
| 105 ## @item | |
| 106 ## @var{eigest} returns the estimate for the smallest @code{@var{eigest} | |
| 107 ## (1)} and largest @code{@var{eigest} (2)} eigenvalues of the | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
108 ## preconditioned matrix @code{@var{P} = @var{m} \ @var{a}}. In |
| 7007 | 109 ## particular, if no preconditioning is used, the estimates for the |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
110 ## extreme eigenvalues of @var{a} are returned. @code{@var{eigest} (1)} |
| 5837 | 111 ## is an overestimate and @code{@var{eigest} (2)} is an underestimate, |
| 112 ## so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound | |
| 113 ## for @code{cond (@var{P}, 2)}, which nevertheless in the limit should | |
| 114 ## theoretically be equal to the actual value of the condition number. | |
| 115 ## The method which computes @var{eigest} works only for symmetric positive | |
| 5838 | 116 ## definite @var{a} and @var{m}, and the user is responsible for |
| 5837 | 117 ## verifying this assumption. |
| 118 ## @end itemize | |
| 119 ## | |
| 120 ## Let us consider a trivial problem with a diagonal matrix (we exploit the | |
| 121 ## sparsity of A) | |
| 122 ## | |
| 123 ## @example | |
| 124 ## @group | |
| 10549 | 125 ## n = 10; |
| 126 ## a = diag (sparse (1:n)); | |
| 127 ## b = rand (n, 1); | |
| 8507 | 128 ## [l, u, p, q] = luinc (a, 1.e-3); |
| 5837 | 129 ## @end group |
| 130 ## @end example | |
| 131 ## | |
| 132 ## @sc{Example 1:} Simplest use of @code{pcg} | |
| 133 ## | |
| 134 ## @example | |
| 135 ## x = pcg(A,b) | |
| 136 ## @end example | |
| 137 ## | |
| 138 ## @sc{Example 2:} @code{pcg} with a function which computes | |
| 5838 | 139 ## @code{@var{a} * @var{x}} |
| 5837 | 140 ## |
| 141 ## @example | |
| 142 ## @group | |
| 8507 | 143 ## function y = apply_a (x) |
| 5837 | 144 ## y = [1:N]'.*x; |
| 145 ## endfunction | |
| 146 ## | |
| 8507 | 147 ## x = pcg ("apply_a", b) |
| 5837 | 148 ## @end group |
| 149 ## @end example | |
| 6747 | 150 ## |
| 151 ## @sc{Example 3:} @code{pcg} with a preconditioner: @var{l} * @var{u} | |
| 152 ## | |
| 153 ## @example | |
| 8507 | 154 ## x = pcg (a, b, 1.e-6, 500, l*u); |
| 6747 | 155 ## @end example |
| 156 ## | |
| 157 ## @sc{Example 4:} @code{pcg} with a preconditioner: @var{l} * @var{u}. | |
| 158 ## Faster than @sc{Example 3} since lower and upper triangular matrices | |
| 159 ## are easier to invert | |
| 160 ## | |
| 161 ## @example | |
| 8507 | 162 ## x = pcg (a, b, 1.e-6, 500, l, u); |
| 6747 | 163 ## @end example |
| 164 ## | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
165 ## @sc{Example 5:} Preconditioned iteration, with full diagnostics. The |
| 5837 | 166 ## preconditioner (quite strange, because even the original matrix |
| 5838 | 167 ## @var{a} is trivial) is defined as a function |
| 5837 | 168 ## |
| 169 ## @example | |
| 170 ## @group | |
| 8507 | 171 ## function y = apply_m (x) |
| 172 ## k = floor (length (x) - 2); | |
| 6547 | 173 ## y = x; |
| 8507 | 174 ## y(1:k) = x(1:k)./[1:k]'; |
| 5837 | 175 ## endfunction |
| 176 ## | |
| 7031 | 177 ## [x, flag, relres, iter, resvec, eigest] = ... |
| 8507 | 178 ## pcg (a, b, [], [], "apply_m"); |
| 6547 | 179 ## semilogy (1:iter+1, resvec); |
| 5837 | 180 ## @end group |
| 181 ## @end example | |
| 182 ## | |
| 6747 | 183 ## @sc{Example 6:} Finally, a preconditioner which depends on a |
| 5837 | 184 ## parameter @var{k}. |
| 185 ## | |
| 186 ## @example | |
| 187 ## @group | |
| 8507 | 188 ## function y = apply_M (x, varargin) |
| 5837 | 189 ## K = varargin@{1@}; |
| 6547 | 190 ## y = x; |
| 191 ## y(1:K) = x(1:K)./[1:K]'; | |
| 7007 | 192 ## endfunction |
| 5837 | 193 ## |
| 194 ## [x, flag, relres, iter, resvec, eigest] = ... | |
| 8507 | 195 ## pcg (A, b, [], [], "apply_m", [], [], 3) |
| 5837 | 196 ## @end group |
| 197 ## @end example | |
| 198 ## | |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
199 ## References: |
| 5837 | 200 ## |
|
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
201 ## @enumerate |
|
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
202 ## @item |
|
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
203 ## 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
|
204 ## SIAM, 1995. (the base PCG algorithm) |
| 10549 | 205 ## |
|
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
206 ## @item |
|
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
207 ## 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
|
208 ## (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
|
209 ## 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
|
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> |
| 217 ## - Add the ability to provide the pre-conditioner as two separate | |
| 218 ## matrices | |
| 5837 | 219 |
| 8507 | 220 function [x, flag, relres, iter, resvec, eigest] = pcg (a, b, tol, maxit, m1, m2, x0, varargin) |
| 5837 | 221 |
| 8507 | 222 ## M = M1*M2 |
| 6747 | 223 |
| 224 if (nargin < 7 || isempty (x0)) | |
| 5838 | 225 x = zeros (size (b)); |
| 5837 | 226 else |
| 227 x = x0; | |
| 228 endif | |
| 229 | |
| 8507 | 230 if (nargin < 5 || isempty (m1)) |
| 231 exist_m1 = 0; | |
| 232 else | |
| 233 exist_m1 = 1; | |
| 234 endif | |
| 6747 | 235 |
| 8507 | 236 if (nargin < 6 || isempty (m2)) |
| 237 exist_m2 = 0; | |
| 238 else | |
| 239 exist_m2 = 1; | |
| 240 endif | |
| 5837 | 241 |
| 5838 | 242 if (nargin < 4 || isempty (maxit)) |
| 243 maxit = min (size (b, 1), 20); | |
| 5837 | 244 endif |
| 245 | |
| 5838 | 246 maxit += 2; |
| 5837 | 247 |
| 5838 | 248 if (nargin < 3 || isempty (tol)) |
| 5837 | 249 tol = 1e-6; |
| 250 endif | |
| 251 | |
| 252 preconditioned_residual_out = false; | |
| 253 if (nargout > 5) | |
| 5838 | 254 T = zeros (maxit, maxit); |
| 5837 | 255 preconditioned_residual_out = true; |
| 256 endif | |
| 257 | |
| 8506 | 258 ## Assume A is positive definite. |
| 259 matrix_positive_definite = true; | |
| 5837 | 260 |
| 5838 | 261 p = zeros (size (b)); |
| 5837 | 262 oldtau = 1; |
| 8507 | 263 if (isnumeric (a)) |
| 8506 | 264 ## A is a matrix. |
| 8507 | 265 r = b - a*x; |
| 8506 | 266 else |
| 267 ## A should be a function. | |
| 8507 | 268 r = b - feval (a, x, varargin{:}); |
| 5837 | 269 endif |
| 270 | |
| 5838 | 271 resvec(1,1) = norm (r); |
| 5837 | 272 alpha = 1; |
| 273 iter = 2; | |
| 274 | |
| 6747 | 275 while (resvec (iter-1,1) > tol * resvec (1,1) && iter < maxit) |
| 8507 | 276 if (exist_m1) |
| 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 |
| 5838 | 294 tau = z' * r; |
| 6747 | 295 resvec (iter-1,2) = sqrt (tau); |
| 5838 | 296 beta = tau / oldtau; |
| 5837 | 297 oldtau = tau; |
| 6747 | 298 p = z + beta * p; |
| 8507 | 299 if (isnumeric (a)) |
| 8506 | 300 ## A is a matrix. |
| 8507 | 301 w = a * p; |
| 8506 | 302 else |
| 303 ## A should be a function. | |
| 8507 | 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; |
| 5837 | 318 ## EVS = eig(T(2:iter-1,2:iter-1)); |
| 319 ## fprintf(stderr,"PCG condest: %g (iteration: %d)\n", max(EVS)/min(EVS),iter); | |
| 320 endif | |
| 6747 | 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]; |
| 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 | |
| 361 resvec (iter-1,2) = sqrt (r' * z); | |
| 5837 | 362 else |
| 5838 | 363 resvec = resvec(:,1); |
| 5837 | 364 endif |
| 365 | |
| 366 flag = 0; | |
| 6747 | 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 | |
| 390 %!demo | |
| 391 %! | |
| 10549 | 392 %! # Simplest usage of pcg (see also 'help pcg') |
| 5837 | 393 %! |
| 10549 | 394 %! N = 10; |
| 395 %! A = diag ([1:N]); b = rand (N, 1); y = A \ b; #y is the true solution | |
| 396 %! x = pcg (A, b); | |
| 397 %! printf('The solution relative error is %g\n', norm (x - y) / norm (y)); | |
| 5837 | 398 %! |
| 10549 | 399 %! # You shouldn't be afraid if pcg issues some warning messages in this |
| 400 %! # example: watch out in the second example, why it takes N iterations | |
| 401 %! # of pcg to converge to (a very accurate, by the way) solution | |
| 5837 | 402 %!demo |
| 403 %! | |
| 10549 | 404 %! # Full output from pcg, except for the eigenvalue estimates |
| 405 %! # We use this output to plot the convergence history | |
| 5837 | 406 %! |
| 10549 | 407 %! N = 10; |
| 408 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution | |
| 409 %! [x, flag, relres, iter, resvec] = pcg (A, b); | |
| 410 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X)); | |
| 411 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)'); | |
| 412 %! semilogy([0:iter], resvec / resvec(1),'o-g'); | |
| 6747 | 413 %! legend('relative residual'); |
| 5837 | 414 %!demo |
| 415 %! | |
| 10549 | 416 %! # Full output from pcg, including the eigenvalue estimates |
| 417 %! # Hilbert matrix is extremely ill conditioned, so pcg WILL have problems | |
| 5837 | 418 %! |
| 10549 | 419 %! N = 10; |
| 420 %! A = hilb (N); b = rand (N, 1); X = A \ b; #X is the true solution | |
| 421 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], 200); | |
| 422 %! printf('The solution relative error is %g\n', norm (x - X) / norm (X)); | |
| 423 %! printf('Condition number estimate is %g\n', eigest(2) / eigest (1)); | |
| 424 %! printf('Actual condition number is %g\n', cond (A)); | |
| 425 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); | |
| 426 %! semilogy([0:iter], resvec,['o-g';'+-r']); | |
| 6747 | 427 %! legend('absolute residual','absolute preconditioned residual'); |
| 5837 | 428 %!demo |
| 429 %! | |
| 10549 | 430 %! # Full output from pcg, including the eigenvalue estimates |
| 431 %! # We use the 1-D Laplacian matrix for A, and cond(A) = O(N^2) | |
| 432 %! # and that's the reasone we need some preconditioner; here we take | |
| 433 %! # a very simple and not powerful Jacobi preconditioner, | |
| 434 %! # which is the diagonal of A | |
| 5837 | 435 %! |
| 10549 | 436 %! N = 100; |
| 437 %! A = zeros (N, N); | |
| 438 %! for i=1 : N - 1 # form 1-D Laplacian matrix | |
| 439 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; | |
| 440 %! endfor | |
| 441 %! b = rand (N, 1); X = A \ b; #X is the true solution | |
| 442 %! maxit = 80; | |
| 443 %! printf('System condition number is %g\n', cond (A)); | |
| 444 %! # No preconditioner: the convergence is very slow! | |
| 5837 | 445 %! |
| 10549 | 446 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit); |
| 447 %! printf('System condition number estimate is %g\n', eigest(2) / eigest(1)); | |
| 448 %! title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)'); | |
| 449 %! semilogy([0:iter], resvec(:,1), 'o-g'); | |
| 6747 | 450 %! legend('NO preconditioning: absolute residual'); |
| 5837 | 451 %! |
| 10549 | 452 %! pause(1); |
| 453 %! # Test Jacobi preconditioner: it will not help much!!! | |
| 5837 | 454 %! |
| 10549 | 455 %! M = diag (diag (A)); # Jacobi preconditioner |
| 456 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); | |
| 457 %! printf('JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1)); | |
| 458 %! hold on; | |
| 459 %! semilogy([0:iter], resvec(:,1), 'o-r'); | |
| 6747 | 460 %! legend('NO preconditioning: absolute residual', ... |
| 461 %! 'JACOBI preconditioner: absolute residual'); | |
| 5837 | 462 %! |
| 10549 | 463 %! pause(1); |
| 464 %! # Test nonoverlapping block Jacobi preconditioner: it will help much! | |
| 5837 | 465 %! |
| 10549 | 466 %! M = zeros (N, N); k = 4; |
| 467 %! for i = 1 : k : N # form 1-D Laplacian matrix | |
| 468 %! M (i:i+k-1, i:i+k-1) = A (i:i+k-1, i:i+k-1); | |
| 469 %! endfor | |
| 470 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], maxit, M); | |
| 471 %! printf('BLOCK JACOBI preconditioned system condition number estimate is %g\n', eigest(2) / eigest(1)); | |
| 472 %! semilogy ([0:iter], resvec(:,1),'o-b'); | |
| 6747 | 473 %! legend('NO preconditioning: absolute residual', ... |
| 474 %! 'JACOBI preconditioner: absolute residual', ... | |
| 475 %! 'BLOCK JACOBI preconditioner: absolute residual'); | |
| 10549 | 476 %! hold off; |
| 5837 | 477 %!test |
| 478 %! | |
| 10549 | 479 %! #solve small diagonal system |
| 5837 | 480 %! |
| 10549 | 481 %! N = 10; |
| 482 %! A = diag ([1:N]); b = rand (N, 1); X = A \ b; #X is the true solution | |
| 483 %! [x, flag] = pcg (A, b, [], N+1); | |
| 484 %! assert(norm (x - X) / norm (X), 0, 1e-10); | |
| 485 %! assert(flag, 0); | |
| 5837 | 486 %! |
| 487 %!test | |
| 488 %! | |
| 10549 | 489 %! #solve small indefinite diagonal system |
| 490 %! #despite A is indefinite, the iteration continues and converges | |
| 491 %! #indefiniteness of A is detected | |
| 5837 | 492 %! |
| 10549 | 493 %! N = 10; |
| 494 %! A = diag([1:N] .* (-ones(1, N) .^ 2)); b = rand (N, 1); X = A \ b; #X is the true solution | |
| 495 %! [x, flag] = pcg (A, b, [], N+1); | |
| 496 %! assert(norm (x - X) / norm (X), 0, 1e-10); | |
| 497 %! assert(flag, 3); | |
| 5837 | 498 %! |
| 499 %!test | |
| 500 %! | |
| 10549 | 501 %! #solve tridiagonal system, do not converge in default 20 iterations |
| 5837 | 502 %! |
| 10549 | 503 %! N = 100; |
| 504 %! A = zeros (N, N); | |
| 505 %! for i = 1 : N - 1 # form 1-D Laplacian matrix | |
| 506 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; | |
| 507 %! endfor | |
| 508 %! b = ones (N, 1); X = A \ b; #X is the true solution | |
| 509 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12); | |
| 510 %! assert(flag); | |
| 511 %! assert(relres > 1.0); | |
| 512 %! assert(iter, 20); #should perform max allowable default number of iterations | |
| 5837 | 513 %! |
| 514 %!test | |
| 515 %! | |
| 10549 | 516 %! #solve tridiagonal system with 'prefect' preconditioner |
| 517 %! #converges in one iteration, so the eigest does not work | |
| 518 %! #and issues a warning | |
| 5837 | 519 %! |
| 10549 | 520 %! N = 100; |
| 521 %! A = zeros (N, N); | |
| 522 %! for i = 1 : N - 1 # form 1-D Laplacian matrix | |
| 523 %! A (i:i+1, i:i+1) = [2 -1; -1 2]; | |
| 524 %! endfor | |
| 525 %! b = ones (N, 1); X = A \ b; #X is the true solution | |
| 526 %! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b); | |
| 527 %! assert(norm (x - X) / norm (X), 0, 1e-6); | |
| 528 %! assert(flag, 0); | |
| 529 %! assert(iter, 1); #should converge in one iteration | |
| 530 %! assert(isnan (eigest), isnan ([NaN, NaN])); | |
| 5837 | 531 %! |