Mercurial > octave

comparison scripts/sparse/pcg.m @ 30893:e1788b1a315f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use "fcn" as preferred abbreviation for "function" in m-files. * accumarray.m, accumdim.m, quadl.m, quadv.m, randi.m, structfun.m, __is_function__.m, uigetfile.m, uimenu.m, uiputfile.m, doc_cache_create.m, colorspace_conversion_input_check.m, imageIO.m, argnames.m, vectorize.m, vectorize.m, normest1.m, inputname.m, nthargout.m, display_info_file.m, decic.m, ode15i.m, ode15s.m, ode23.m, ode23s.m, ode45.m, odeset.m, check_default_input.m, integrate_adaptive.m, ode_event_handler.m, runge_kutta_23.m, runge_kutta_23s.m, runge_kutta_45_dorpri.m, runge_kutta_interpolate.m, starting_stepsize.m, __all_opts__.m, fminbnd.m, fminsearch.m, fminunc.m, fsolve.m, fzero.m, sqp.m, fplot.m, plotyy.m, __bar__.m, __ezplot__.m, flat_entry.html, profexport.m, movfun.m, bicg.m, bicgstab.m, cgs.m, eigs.m, gmres.m, pcg.m, __alltohandles__.m, __sprand__.m, qmr.m, tfqmr.m, dump_demos.m: Replace "func", "fun", "fn" in documentation and variable names with "fcn".
author Rik <rik@octave.org>
date Mon, 04 Apr 2022 18:14:56 -0700
parents 796f54d4ddbf
children 597f3ee61a48
comparison
equal deleted inserted replaced
30892:1a3cc2811090 30893:e1788b1a315f
34 ## The input arguments are: 34 ## The input arguments are:
35 ## 35 ##
36 ## @itemize 36 ## @itemize
37 ## @item @var{A} is the matrix of the linear system and it must be square. 37 ## @item @var{A} is the matrix of the linear system and it must be square.
38 ## @var{A} can be passed as a matrix, function handle, or inline function 38 ## @var{A} can be passed as a matrix, function handle, or inline function
39 ## @code{Afun} such that @code{Afun(x) = A * x}. Additional parameters to 39 ## @code{Afcn} such that @code{Afcn(x) = A * x}. Additional parameters to
40 ## @code{Afun} may be passed after @var{x0}. 40 ## @code{Afcn} may be passed after @var{x0}.
41 ## 41 ##
42 ## @var{A} has to be Hermitian and Positive Definite (@nospell{HPD})@. If 42 ## @var{A} has to be Hermitian and Positive Definite (@nospell{HPD})@. If
43 ## @code{pcg} detects @var{A} not to be positive definite, a warning is printed 43 ## @code{pcg} detects @var{A} not to be positive definite, a warning is printed
44 ## and the @var{flag} output is set. 44 ## and the @var{flag} output is set.
45 ## 45 ##
161 ## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n)); 161 ## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
162 ## b = A * ones (n, 1); 162 ## b = A * ones (n, 1);
163 ## M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)' 163 ## M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
164 ## M2 = M1'; 164 ## M2 = M1';
165 ## M = M1 * M2; 165 ## M = M1 * M2;
166 ## Afun = @@(x) A * x; 166 ## Afcn = @@(x) A * x;
167 ## Mfun = @@(x) M \ x; 167 ## Mfcn = @@(x) M \ x;
168 ## M1fun = @@(x) M1 \ x; 168 ## M1fcn = @@(x) M1 \ x;
169 ## M2fun = @@(x) M2 \ x; 169 ## M2fcn = @@(x) M2 \ x;
170 ## @end group 170 ## @end group
171 ## @end example 171 ## @end example
172 ## 172 ##
173 ## @sc{Example 1:} Simplest use of @code{pcg} 173 ## @sc{Example 1:} Simplest use of @code{pcg}
174 ## 174 ##
178 ## 178 ##
179 ## @sc{Example 2:} @code{pcg} with a function which computes 179 ## @sc{Example 2:} @code{pcg} with a function which computes
180 ## @code{@var{A} * @var{x}} 180 ## @code{@var{A} * @var{x}}
181 ## 181 ##
182 ## @example 182 ## @example
183 ## x = pcg (Afun, b) 183 ## x = pcg (Afcn, b)
184 ## @end example 184 ## @end example
185 ## 185 ##
186 ## @sc{Example 3:} @code{pcg} with a preconditioner matrix @var{M} 186 ## @sc{Example 3:} @code{pcg} with a preconditioner matrix @var{M}
187 ## 187 ##
188 ## @example 188 ## @example
190 ## @end example 190 ## @end example
191 ## 191 ##
192 ## @sc{Example 4:} @code{pcg} with a function as preconditioner 192 ## @sc{Example 4:} @code{pcg} with a function as preconditioner
193 ## 193 ##
194 ## @example 194 ## @example
195 ## x = pcg (Afun, b, 1e-6, 100, Mfun) 195 ## x = pcg (Afcn, b, 1e-6, 100, Mfcn)
196 ## @end example 196 ## @end example
197 ## 197 ##
198 ## @sc{Example 5:} @code{pcg} with preconditioner matrices @var{M1} 198 ## @sc{Example 5:} @code{pcg} with preconditioner matrices @var{M1}
199 ## and @var{M2} 199 ## and @var{M2}
200 ## 200 ##
203 ## @end example 203 ## @end example
204 ## 204 ##
205 ## @sc{Example 6:} @code{pcg} with functions as preconditioners 205 ## @sc{Example 6:} @code{pcg} with functions as preconditioners
206 ## 206 ##
207 ## @example 207 ## @example
208 ## x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun) 208 ## x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn)
209 ## @end example 209 ## @end example
210 ## 210 ##
211 ## @sc{Example 7:} @code{pcg} with as input a function requiring an argument 211 ## @sc{Example 7:} @code{pcg} with as input a function requiring an argument
212 ## 212 ##
213 ## @example 213 ## @example
216 ## y = x; 216 ## y = x;
217 ## for i = 1:p 217 ## for i = 1:p
218 ## y = A * y; 218 ## y = A * y;
219 ## endfor 219 ## endfor
220 ## endfunction 220 ## endfunction
221 ## Apfun = @@(x, p) Ap (A, x, p); 221 ## Apfcn = @@(x, p) Ap (A, x, p);
222 ## x = pcg (Apfun, b, [], [], [], [], [], 2); 222 ## x = pcg (Apfcn, b, [], [], [], [], [], 2);
223 ## @end group 223 ## @end group
224 ## @end example 224 ## @end example
225 ## 225 ##
226 ## @sc{Example 8:} explicit example to show that @code{pcg} uses a 226 ## @sc{Example 8:} explicit example to show that @code{pcg} uses a
227 ## split preconditioner 227 ## split preconditioner
273 endif 273 endif
274 274
275 ## Check if the input data A,b,m1,m2 are consistent (i.e. if they are 275 ## Check if the input data A,b,m1,m2 are consistent (i.e. if they are
276 ## matrix or function handle) 276 ## matrix or function handle)
277 277
278 [Afun, M1fun, M2fun] = __alltohandles__ (A, b, M1, M2, "pcg"); 278 [Afcn, M1fcn, M2fcn] = __alltohandles__ (A, b, M1, M2, "pcg");
279 279
280 maxit += 2; 280 maxit += 2;
281 n_arg_out = nargout; 281 n_arg_out = nargout;
282 282
283 ## Set Initial data 283 ## Set Initial data
299 x = x_pr = x_min = x0; 299 x = x_pr = x_min = x0;
300 300
301 ## x_pr (x previous) needs to check the stagnation 301 ## x_pr (x previous) needs to check the stagnation
302 ## x_min needs to save the iterated with minimum residual 302 ## x_min needs to save the iterated with minimum residual
303 303
304 r = b - feval (Afun, x, varargin{:}); 304 r = b - feval (Afcn, x, varargin{:});
305 iter = 2; 305 iter = 2;
306 iter_min = 0; 306 iter_min = 0;
307 flag = 1; 307 flag = 1;
308 resvec = zeros (maxit + 1, 2); 308 resvec = zeros (maxit + 1, 2);
309 resvec(1, 1) = norm (r); 309 resvec(1, 1) = norm (r);
318 318
319 while (resvec(iter-1,1) > tol * b_norm && iter < maxit) 319 while (resvec(iter-1,1) > tol * b_norm && iter < maxit)
320 if (iter == 2) # Check whether M1 or M2 are singular 320 if (iter == 2) # Check whether M1 or M2 are singular
321 try 321 try
322 warning ("error","Octave:singular-matrix","local"); 322 warning ("error","Octave:singular-matrix","local");
323 z = feval (M1fun, r, varargin{:}); 323 z = feval (M1fcn, r, varargin{:});
324 z = feval (M2fun, z, varargin{:}); 324 z = feval (M2fcn, z, varargin{:});
325 catch 325 catch
326 flag = 2; 326 flag = 2;
327 break; 327 break;
328 end_try_catch 328 end_try_catch
329 else 329 else
330 z = feval (M1fun, r, varargin{:}); 330 z = feval (M1fcn, r, varargin{:});
331 z = feval (M2fun, z, varargin{:}); 331 z = feval (M2fcn, z, varargin{:});
332 endif 332 endif
333 333
334 tau = z' * r; 334 tau = z' * r;
335 resvec(iter - 1, 2) = sqrt (tau); 335 resvec(iter - 1, 2) = sqrt (tau);
336 beta = tau / old_tau; 336 beta = tau / old_tau;
337 old_tau = tau; 337 old_tau = tau;
338 p = z + beta * p; 338 p = z + beta * p;
339 w = feval (Afun, p, varargin{:}); 339 w = feval (Afcn, p, varargin{:});
340 340
341 ## Needed only for eigest. 341 ## Needed only for eigest.
342 342
343 old_alpha = alpha; 343 old_alpha = alpha;
344 den = p' * w; 344 den = p' * w;
376 endwhile 376 endwhile
377 377
378 if (n_arg_out > 5) 378 if (n_arg_out > 5)
379 ## Apply the preconditioner once more and finish with the precond 379 ## Apply the preconditioner once more and finish with the precond
380 ## residual. 380 ## residual.
381 z = feval (M1fun, r, varargin{:}); 381 z = feval (M1fcn, r, varargin{:});
382 z = feval (M2fun, z, varargin{:}); 382 z = feval (M2fcn, z, varargin{:});
383 endif 383 endif
384 384
385 ## (Eventually) computes the eigenvalue of inv(m2)*inv(m1)*A 385 ## (Eventually) computes the eigenvalue of inv(m2)*inv(m1)*A
386 if (n_arg_out > 5) 386 if (n_arg_out > 5)
387 if (flag != 4) 387 if (flag != 4)
464 %! b = A * ones (n, 1); 464 %! b = A * ones (n, 1);
465 %! M1 = ichol (A); # for this tridiagonal case it corresponds to chol (A)' 465 %! M1 = ichol (A); # for this tridiagonal case it corresponds to chol (A)'
466 %! M2 = M1'; 466 %! M2 = M1';
467 %! M = M1 * M2; 467 %! M = M1 * M2;
468 %! x = pcg (A, b); 468 %! x = pcg (A, b);
469 %! Afun = @(x) A * x; 469 %! Afcn = @(x) A * x;
470 %! x = pcg (Afun, b); 470 %! x = pcg (Afcn, b);
471 %! x = pcg (A, b, 1e-6, 100, M); 471 %! x = pcg (A, b, 1e-6, 100, M);
472 %! x = pcg (A, b, 1e-6, 100, M1, M2); 472 %! x = pcg (A, b, 1e-6, 100, M1, M2);
473 %! Mfun = @(x) M \ x; 473 %! Mfcn = @(x) M \ x;
474 %! x = pcg (Afun, b, 1e-6, 100, Mfun); 474 %! x = pcg (Afcn, b, 1e-6, 100, Mfcn);
475 %! M1fun = @(x) M1 \ x; 475 %! M1fcn = @(x) M1 \ x;
476 %! M2fun = @(x) M2 \ x; 476 %! M2fcn = @(x) M2 \ x;
477 %! x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun); 477 %! x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn);
478 %! function y = Ap (A, x, p) # compute A^p * x 478 %! function y = Ap (A, x, p) # compute A^p * x
479 %! y = x; 479 %! y = x;
480 %! for i = 1:p 480 %! for i = 1:p
481 %! y = A * y; 481 %! y = A * y;
482 %! endfor 482 %! endfor
483 %! endfunction 483 %! endfunction
484 %! Afun = @(x, p) Ap (A, x, p); 484 %! Afcn = @(x, p) Ap (A, x, p);
485 %! ## solution of A^2 * x = b 485 %! ## solution of A^2 * x = b
486 %! x = pcg (Afun, b, [], [], [], [], [], 2); 486 %! x = pcg (Afcn, b, [], [], [], [], [], 2);
487 487
488 %!demo 488 %!demo
489 %! n = 10; 489 %! n = 10;
490 %! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n)); 490 %! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
491 %! b = A * ones (n, 1); 491 %! b = A * ones (n, 1);
504 %! ## Check that all type of inputs work 504 %! ## Check that all type of inputs work
505 %! A = toeplitz (sparse ([2, 1 ,0, 0, 0])); 505 %! A = toeplitz (sparse ([2, 1 ,0, 0, 0]));
506 %! b = A * ones (5, 1); 506 %! b = A * ones (5, 1);
507 %! M1 = diag (sqrt (diag (A))); 507 %! M1 = diag (sqrt (diag (A)));
508 %! M2 = M1; # M1 * M2 is the Jacobi preconditioner 508 %! M2 = M1; # M1 * M2 is the Jacobi preconditioner
509 %! Afun = @(z) A*z; 509 %! Afcn = @(z) A*z;
510 %! M1_fun = @(z) M1 \ z; 510 %! M1_fcn = @(z) M1 \ z;
511 %! M2_fun = @(z) M2 \ z; 511 %! M2_fcn = @(z) M2 \ z;
512 %! [x, flag, ~, iter] = pcg (A,b); 512 %! [x, flag, ~, iter] = pcg (A,b);
513 %! assert (flag, 0); 513 %! assert (flag, 0);
514 %! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2); 514 %! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2);
515 %! assert (flag, 0); 515 %! assert (flag, 0);
516 %! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2); 516 %! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2);
517 %! assert (flag, 0); 517 %! assert (flag, 0);
518 %! [x, flag] = pcg (A, b, [], [], M1_fun, M2_fun); 518 %! [x, flag] = pcg (A, b, [], [], M1_fcn, M2_fcn);
519 %! assert (flag, 0); 519 %! assert (flag, 0);
520 %! [x, flag] = pcg (A, b,[],[], M1_fun, M2); 520 %! [x, flag] = pcg (A, b,[],[], M1_fcn, M2);
521 %! assert (flag, 0); 521 %! assert (flag, 0);
522 %! [x, flag] = pcg (A, b,[],[], M1, M2_fun); 522 %! [x, flag] = pcg (A, b,[],[], M1, M2_fcn);
523 %! assert (flag, 0); 523 %! assert (flag, 0);
524 %! [x, flag] = pcg (Afun, b); 524 %! [x, flag] = pcg (Afcn, b);
525 %! assert (flag, 0); 525 %! assert (flag, 0);
526 %! [x, flag] = pcg (Afun, b,[],[], M1 * M2); 526 %! [x, flag] = pcg (Afcn, b,[],[], M1 * M2);
527 %! assert (flag, 0); 527 %! assert (flag, 0);
528 %! [x, flag] = pcg (Afun, b,[],[], M1, M2); 528 %! [x, flag] = pcg (Afcn, b,[],[], M1, M2);
529 %! assert (flag, 0); 529 %! assert (flag, 0);
530 %! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2); 530 %! [x, flag] = pcg (Afcn, b,[],[], M1_fcn, M2);
531 %! assert (flag, 0); 531 %! assert (flag, 0);
532 %! [x, flag] = pcg (Afun, b,[],[], M1, M2_fun); 532 %! [x, flag] = pcg (Afcn, b,[],[], M1, M2_fcn);
533 %! assert (flag, 0); 533 %! assert (flag, 0);
534 %! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2_fun); 534 %! [x, flag] = pcg (Afcn, b,[],[], M1_fcn, M2_fcn);
535 %! assert (flag, 0); 535 %! assert (flag, 0);
536 536
537 %!test 537 %!test
538 %! ## solve a small diagonal system 538 %! ## solve a small diagonal system
539 %! N = 10; 539 %! N = 10;
656 %! b = single (1); 656 %! b = single (1);
657 %! [x, flag] = pcg (A, b); 657 %! [x, flag] = pcg (A, b);
658 %! assert (class (x), "single"); 658 %! assert (class (x), "single");
659 659
660 %!test 660 %!test
661 %!function y = Afun (x) 661 %!function y = Afcn (x)
662 %! A = toeplitz ([2, 1, 0, 0]); 662 %! A = toeplitz ([2, 1, 0, 0]);
663 %! y = A * x; 663 %! y = A * x;
664 %!endfunction 664 %!endfunction
665 %! [x, flag] = pcg ("Afun", [3; 4; 4; 3]); 665 %! [x, flag] = pcg ("Afcn", [3; 4; 4; 3]);
666 %! assert (x, ones (4, 1), 1e-6); 666 %! assert (x, ones (4, 1), 1e-6);
667 667
668 %!test 668 %!test
669 %! ## unpreconditioned residual 669 %! ## unpreconditioned residual
670 %! A = toeplitz (sparse ([4, 1, 0, 0, 0])); 670 %! A = toeplitz (sparse ([4, 1, 0, 0, 0]));