Mercurial > octave-antonio

annotate scripts/sparse/bicg.m @ 19700:00e31f316a3a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix Matlab incompatibility of "ismatrix" (bug #42422). * data.cc (isvector): new tests * data.cc (isrow): documentation improved, new tests * data.cc (iscolumn): documentation improved, new tests * data.cc (ismatrix): is matrix now only checks the dimension due to Matlab compatibility, documentation improved, new tests * accumarray.m: use more appropriate function, than ismatrix * gradient.m: use more appropriate function, than ismatrix * num2str.m: use more appropriate functions, than ismatrix * ntsc2rgb.m: use more appropriate function, than ismatrix * condest.m: use more appropriate function, than ismatrix * expm.m: use more appropriate function, than ismatrix * onenormest.m: use more appropriate function, than ismatrix * isocolors.m: use more appropriate function, than ismatrix * isonormals.m: use more appropriate function, than ismatrix * isosurface.m: use more appropriate function, than ismatrix * __errcomm__.m: use more appropriate function, than ismatrix * __interp_cube__.m: use more appropriate function, than ismatrix * __marching_cube__.m: use more appropriate function, than ismatrix * __stem__.m: use more appropriate function, than ismatrix * stairs.m: use more appropriate function, than ismatrix * validsetargs.m: use more appropriate functions, than ismatrix * unique.m: use more appropriate functions, than ismatrix * bicg.m: additional tests for numerical value * bicgstab.m: additional tests for numerical value * cgs.m: additional tests for numerical value * gmres.m: additional tests for numerical value * qmr.m: additional tests for numerical value
author Kai T. Ohlhus <k.ohlhus@gmail.com>
date Thu, 12 Feb 2015 18:34:56 +0100
parents 4197fc428c7d
children 5f2c0ca0ef51
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
16768
e2de3c8882be copyright notice fixes
John W. Eaton <jwe@octave.org>
parents: 14853
diff changeset
1 ## Copyright (C) 2006 Sylvain Pelissier
19697
4197fc428c7d maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents: 17795
diff changeset
2 ## Copyright (C) 2012-2015 Carlo de Falco
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
3 ##
17795
0a8c35ae5ce1 maint: Fix various problems with GPL copyright statements.
Rik <rik@octave.org>
parents: 17744
diff changeset
4 ## This file is part of Octave.
0a8c35ae5ce1 maint: Fix various problems with GPL copyright statements.
Rik <rik@octave.org>
parents: 17744
diff changeset
5 ##
0a8c35ae5ce1 maint: Fix various problems with GPL copyright statements.
Rik <rik@octave.org>
parents: 17744
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
7 ## it under the terms of the GNU General Public License as published by
17795
0a8c35ae5ce1 maint: Fix various problems with GPL copyright statements.
Rik <rik@octave.org>
parents: 17744
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
9 ## (at your option) any later version.
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
10 ##
17795
0a8c35ae5ce1 maint: Fix various problems with GPL copyright statements.
Rik <rik@octave.org>
parents: 17744
diff changeset
11 ## Octave is distributed in the hope that it will be useful,
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
12 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
14 ## GNU General Public License for more details.
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
15 ##
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
16 ## You should have received a copy of the GNU General Public License
17795
0a8c35ae5ce1 maint: Fix various problems with GPL copyright statements.
Rik <rik@octave.org>
parents: 17744
diff changeset
17 ## along with Octave; If not, see <http://www.gnu.org/licenses/>.
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
18
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
19 ## -*- texinfo -*-
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
20 ## @deftypefn {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
21 ## @deftypefnx {Function File} {@var{x} =} bicg (@var{A}, @var{b}, @var{rtol}, @var{maxit}, @var{P})
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
22 ## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} bicg (@var{A}, @var{b}, @dots{})
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
23 ## Solve @code{A x = b} using the Bi-conjugate gradient iterative method.
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
24 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
25 ## @itemize @minus
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
26 ## @item @var{rtol} is the relative tolerance, if not given
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
27 ## or set to [] the default value 1e-6 is used.
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
28 ##
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
29 ## @item @var{maxit} the maximum number of outer iterations,
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
30 ## if not given or set to [] the default value
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
31 ## @code{min (20, numel (b))} is used.
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
32 ##
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
33 ## @item @var{x0} the initial guess, if not given or set to []
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
34 ## the default value @code{zeros (size (b))} is used.
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
35 ## @end itemize
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
36 ##
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
37 ## @var{A} can be passed as a matrix or as a function handle or
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
38 ## inline function @code{f} such that @code{f(x, "notransp") = A*x}
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
39 ## and @code{f(x, "transp") = A'*x}.
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
40 ##
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
41 ## The preconditioner @var{P} is given as @code{P = M1 * M2}.
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
42 ## Both @var{M1} and @var{M2} can be passed as a matrix or as
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
43 ## a function handle or inline function @code{g} such that
14359
7277fe922e99 doc: Use Octave preference for double quote in docstrings in scripts/
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
44 ## @code{g(x, "notransp") = M1 \ x} or @code{g(x, "notransp") = M2 \ x} and
7277fe922e99 doc: Use Octave preference for double quote in docstrings in scripts/
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
45 ## @code{g(x, "transp") = M1' \ x} or @code{g(x, "transp") = M2' \ x}.
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
46 ##
13931
9de488c6c59c doc: Spellcheck documentation before 3.6.0 release
Rik <octave@nomad.inbox5.com>
parents: 13929
diff changeset
47 ## If called with more than one output parameter
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
48 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
49 ## @itemize @minus
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
50 ## @item @var{flag} indicates the exit status:
14853
72b8b39e12be doc: Periodic grammarcheck of documentation.
Rik <octave@nomad.inbox5.com>
parents: 14609
diff changeset
51 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
52 ## @itemize @minus
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
53 ## @item 0: iteration converged to the within the chosen tolerance
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
54 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
55 ## @item 1: the maximum number of iterations was reached before convergence
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
56 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
57 ## @item 3: the algorithm reached stagnation
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
58 ## @end itemize
16816
12005245b645 doc: Periodic grammarcheck of documentation.
Rik <rik@octave.org>
parents: 16768
diff changeset
59 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
60 ## (the value 2 is unused but skipped for compatibility).
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
61 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
62 ## @item @var{relres} is the final value of the relative residual.
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
63 ##
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
64 ## @item @var{iter} is the number of iterations performed.
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
65 ##
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
66 ## @item @var{resvec} is a vector containing the relative residual at each iteration.
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
67 ## @end itemize
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
68 ##
13929
9cae456085c2 Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents: 13796
diff changeset
69 ## @seealso{bicgstab, cgs, gmres, pcg}
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
70 ##
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
71 ## @end deftypefn
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
72
16768
e2de3c8882be copyright notice fixes
John W. Eaton <jwe@octave.org>
parents: 14853
diff changeset
73 ## Author: Sylvain Pelissier <sylvain.pelissier@gmail.com>
e2de3c8882be copyright notice fixes
John W. Eaton <jwe@octave.org>
parents: 14853
diff changeset
74 ## Author: Carlo de Falco
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
75
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
76 function [x, flag, res1, k, resvec] = bicg (A, b, tol, maxit, M1, M2, x0)
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
77
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
78 if (nargin >= 2 && isvector (full (b)))
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
79
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
80 if (ischar (A))
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
81 fun = str2func (A);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
82 Ax = @(x) feval (fun, x, "notransp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
83 Atx = @(x) feval (fun, x, "transp");
19700
00e31f316a3a Fix Matlab incompatibility of "ismatrix" (bug #42422).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents: 19697
diff changeset
84 elseif (isnumeric (A) && ismatrix (A))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
85 Ax = @(x) A * x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
86 Atx = @(x) A' * x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
87 elseif (isa (A, "function_handle"))
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
88 Ax = @(x) feval (A, x, "notransp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
89 Atx = @(x) feval (A, x, "transp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
90 else
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
91 error (["bicg: first argument is expected to " ...
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
92 "be a function or a square matrix"]);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
93 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
94
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
95 if (nargin < 3 || isempty (tol))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
96 tol = 1e-6;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
97 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
98
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
99 if (nargin < 4 || isempty (maxit))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
100 maxit = min (rows (b), 20);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
101 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
102
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
103 if (nargin < 5 || isempty (M1))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
104 M1m1x = @(x, ignore) x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
105 M1tm1x = M1m1x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
106 elseif (ischar (M1))
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
107 fun = str2func (M1);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
108 M1m1x = @(x) feval (fun, x, "notransp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
109 M1tm1x = @(x) feval (fun, x, "transp");
19700
00e31f316a3a Fix Matlab incompatibility of "ismatrix" (bug #42422).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents: 19697
diff changeset
110 elseif (isnumeric (M1) && ismatrix (M1))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
111 M1m1x = @(x) M1 \ x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
112 M1tm1x = @(x) M1' \ x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
113 elseif (isa (M1, "function_handle"))
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
114 M1m1x = @(x) feval (M1, x, "notransp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
115 M1tm1x = @(x) feval (M1, x, "transp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
116 else
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
117 error (["bicg: preconditioner is expected to " ...
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
118 "be a function or matrix"]);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
119 endif
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
120
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
121 if (nargin < 6 || isempty (M2))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
122 M2m1x = @(x, ignore) x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
123 M2tm1x = M2m1x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
124 elseif (ischar (M2))
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
125 fun = str2func (M2);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
126 M2m1x = @(x) feval (fun, x, "notransp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
127 M2tm1x = @(x) feval (fun, x, "transp");
19700
00e31f316a3a Fix Matlab incompatibility of "ismatrix" (bug #42422).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents: 19697
diff changeset
128 elseif (isnumeric (M2) && ismatrix (M2))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
129 M2m1x = @(x) M2 \ x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
130 M2tm1x = @(x) M2' \ x;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
131 elseif (isa (M2, "function_handle"))
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
132 M2m1x = @(x) feval (M2, x, "notransp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
133 M2tm1x = @(x) feval (M2, x, "transp");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
134 else
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
135 error (["bicg: preconditioner is expected to " ...
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
136 "be a function or matrix"]);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
137 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
138
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
139 Pm1x = @(x) M2m1x (M1m1x (x));
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
140 Ptm1x = @(x) M1tm1x (M2tm1x (x));
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
141
13091
e5aaba072d2b maint: style fixes in iterative linear solvers
Carlo de Falco <kingcrimson@tiscali.it>
parents: 13073
diff changeset
142 if (nargin < 7 || isempty (x0))
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
143 x0 = zeros (size (b));
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
144 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
145
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
146 y = x = x0;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
147 c = b;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
148
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
149 r0 = b - Ax (x);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
150 s0 = c - Atx (y);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
151
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
152 d = Pm1x (r0);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
153 f = Ptm1x (s0);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
154
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
155 bnorm = norm (b);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
156 res0 = Inf;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
157
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
158 if (any (r0 != 0))
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
159
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
160 for k = 1:maxit
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
161
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
162 a = (s0' * Pm1x (r0)) ./ (f' * Ax (d));
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
163
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
164 x += a * d;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
165 y += conj (a) * f;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
166
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
167 r1 = r0 - a * Ax (d);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
168 s1 = s0 - conj (a) * Atx (f);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
169
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
170 beta = (s1' * Pm1x (r1)) ./ (s0' * Pm1x (r0));
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
171
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
172 d = Pm1x (r1) + beta * d;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
173 f = Ptm1x (s1) + conj (beta) * f;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
174
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
175 r0 = r1;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
176 s0 = s1;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
177
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
178 res1 = norm (b - Ax (x)) / bnorm;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
179 if (res1 < tol)
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
180 flag = 0;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
181 if (nargout < 2)
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
182 printf ("bicg converged at iteration %i ", k);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
183 printf ("to a solution with relative residual %e\n", res1);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
184 endif
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
185 break;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
186 endif
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
187
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
188 if (res0 <= res1)
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
189 flag = 3;
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
190 printf ("bicg stopped at iteration %i ", k);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
191 printf ("without converging to the desired tolerance %e\n", tol);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
192 printf ("because the method stagnated.\n");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
193 printf ("The iterate returned (number %i) ", k-1);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
194 printf ("has relative residual %e\n", res0);
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
195 break
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
196 endif
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
197 res0 = res1;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
198 if (nargout > 4)
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
199 resvec(k) = res0;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
200 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
201 endfor
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
202
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
203 if (k == maxit)
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
204 flag = 1;
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
205 printf ("bicg stopped at iteration %i ", maxit);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
206 printf ("without converging to the desired tolerance %e\n", tol);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
207 printf ("because the maximum number of iterations was reached. ");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
208 printf ("The iterate returned (number %i) has ", maxit);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
209 printf ("relative residual %e\n", res1);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
210 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
211
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
212 else
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
213 flag = 0;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
214 if (nargout < 2)
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
215 printf ("bicg converged after 0 interations\n");
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
216 endif
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
217 endif
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
218
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
219 else
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
220 print_usage ();
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
221 endif
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
222
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
223 endfunction;
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
224
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
225
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
226 %!test
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
227 %! n = 100;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
228 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
229 %! b = sum (A, 2);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
230 %! tol = 1e-8;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
231 %! maxit = 15;
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
232 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
233 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
234 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, maxit, M1, M2);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
235 %! assert (x, ones (size (b)), 1e-7);
13796
53674ceb9133 maint: fix function definition in test
John W. Eaton <jwe@octave.org>
parents: 13141
diff changeset
236
53674ceb9133 maint: fix function definition in test
John W. Eaton <jwe@octave.org>
parents: 13141
diff changeset
237 %!function y = afun (x, t, a)
16933
e39f00a32dc7 maint: Use parentheses around condition for switch(),while(),if() statements.
Rik <rik@octave.org>
parents: 16816
diff changeset
238 %! switch (t)
17174
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 16933
diff changeset
239 %! case "notransp"
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 16933
diff changeset
240 %! y = a * x;
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 16933
diff changeset
241 %! case "transp"
c3c1ebfaa7dc maint: Use common indentation for switch statement.
Rik <rik@octave.org>
parents: 16933
diff changeset
242 %! y = a' * x;
13796
53674ceb9133 maint: fix function definition in test
John W. Eaton <jwe@octave.org>
parents: 13141
diff changeset
243 %! endswitch
53674ceb9133 maint: fix function definition in test
John W. Eaton <jwe@octave.org>
parents: 13141
diff changeset
244 %!endfunction
53674ceb9133 maint: fix function definition in test
John W. Eaton <jwe@octave.org>
parents: 13141
diff changeset
245 %!
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
246 %!test
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
247 %! n = 100;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
248 %! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n);
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
249 %! b = sum (A, 2);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
250 %! tol = 1e-8;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
251 %! maxit = 15;
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
252 %! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n);
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
253 %! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
254 %!
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
255 %! [x, flag, relres, iter, resvec] = bicg (@(x, t) afun (x, t, A),
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
256 %! b, tol, maxit, M1, M2);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
257 %! assert (x, ones (size (b)), 1e-7);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
258
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
259 %!test
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
260 %! n = 100;
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
261 %! tol = 1e-8;
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
262 %! a = sprand (n, n, .1);
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
263 %! A = a' * a + 100 * eye (n);
13141
e81ddf9cacd5 maint: untabify and remove trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 13091
diff changeset
264 %! b = sum (A, 2);
13021
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
265 %! [x, flag, relres, iter, resvec] = bicg (A, b, tol, [], diag (diag (A)));
d55d396a9a55 Add an implementation of the biconjugate gradient iterative method
Carlo de Falco <kingcrimson@tiscali.it>
parents:
diff changeset
266 %! assert (x, ones (size (b)), 1e-7);
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14359
diff changeset
267