Mercurial > hg > octave-lyh
annotate scripts/sparse/cgs.m @ 9278:2b35cb600d50
improve bicgstab and cgs
author | Radek Salac <salac.r@gmail.com> |
---|---|
date | Thu, 28 May 2009 07:03:11 +0200 |
parents | eb63fbe60fab |
children | 1673a0dc019f |
rev | line source |
---|---|
9278 | 1 ## Copyright (C) 2008 Radek Salac |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
2 ## |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
3 ## This file is part of Octave. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
4 ## |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
5 ## Octave is free software; you can redistribute it and/or modify it |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
6 ## under the terms of the GNU General Public License as published by |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
8 ## your option) any later version. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
9 ## |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
10 ## Octave is distributed in the hope that it will be useful, but |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
13 ## General Public License for more details. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
14 ## |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
15 ## You should have received a copy of the GNU General Public License |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
16 ## along with Octave; see the file COPYING. If not, see |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
17 ## <http://www.gnu.org/licenses/>. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
18 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
19 ## -*- texinfo -*- |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
20 ## @deftypefn {Function File} {} cgs (@var{A}, @var{b}) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
21 ## @deftypefnx {Function File} {} cgs (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
22 ## This procedure attempts to solve a system of linear equations A*x = b for x. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
23 ## The @var{A} must be square, symmetric and positive definite real matrix N*N. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
24 ## The @var{b} must be a one column vector with a length of N. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
25 ## The @var{tol} specifies the tolerance of the method, default value is 1e-6. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
26 ## The @var{maxit} specifies the maximum number of iteration, default value is MIN(20,N). |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
27 ## The @var{M1} specifies a preconditioner, can also be a function handler which returns M\X. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
28 ## The @var{M2} combined with @var{M1} defines preconditioner as preconditioner=M1*M2. |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
29 ## The @var{x0} is initial guess, default value is zeros(N,1). |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
30 ## |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
31 ## @end deftypefn |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
32 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
33 function [x, flag, relres, iter, resvec] = cgs (A, b, tol, maxit, M1, M2, x0) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
34 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
35 if (nargin < 2 || nargin > 7 || nargout > 5) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
36 print_usage (); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
37 elseif (!isnumeric (A) || rows (A) != columns (A)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
38 error ("cgs: first argument must be a n-by-n matrix"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
39 elseif (!isvector (b)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
40 error ("cgs: b must be a vector"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
41 elseif (rows (A) != rows (b)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
42 error ("cgs: first and second argument must have the same number of rows"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
43 elseif (nargin > 2 && !isscalar (tol)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
44 error ("cgs: tol must be a scalar"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
45 elseif (nargin > 3 && !isscalar (maxit)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
46 error ("cgs: maxit must be a scalar"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
47 elseif (nargin > 4 && ismatrix (M1) && (rows (M1) != rows (A) || columns (M1) != columns (A))) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
48 error ("cgs: M1 must have the same number of rows and columns as A"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
49 elseif (nargin > 5 && (!ismatrix (M2) || rows (M2) != rows (A) || columns (M2) != columns (A))) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
50 error ("cgs: M2 must have the same number of rows and columns as A"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
51 elseif (nargin > 6 && !isvector (x0)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
52 error ("cgs: x0 must be a vector"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
53 elseif (nargin > 6 && rows (x0) != rows (b)) |
8507 | 54 error ("cgs: x0 must have the same number of rows as b"); |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
55 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
56 |
8507 | 57 ## Default tolerance. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
58 if (nargin < 3) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
59 tol = 1e-6; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
60 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
61 |
8507 | 62 ## Default maximum number of iteration. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
63 if (nargin < 4) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
64 maxit = min (rows (b),20); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
65 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
66 |
8507 | 67 ## Left preconditioner. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
68 precon = []; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
69 if (nargin == 5) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
70 precon = M1; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
71 elseif (nargin > 5) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
72 if (isparse(M1) && issparse(M2)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
73 precon = @(x) M1 * (M2 * x); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
74 else |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
75 precon = M1 * M2; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
76 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
77 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
78 |
9278 | 79 if (nargin > 4 && isnumeric(precon) ) |
80 precon = inv(precon); | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
81 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
82 |
8507 | 83 ## Specifies initial estimate x0. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
84 if (nargin < 7) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
85 x = zeros (rows (b), 1); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
86 else |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
87 x = x0; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
88 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
89 |
9278 | 90 res = b - A * x; |
91 norm_b = norm (b); | |
8507 | 92 ## Vector of the residual norms for each iteration. |
9278 | 93 resvec = [ norm(res)/norm_b ]; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
94 ro = 0; |
8507 | 95 ## Default behavior we don't reach tolerance tol within maxit iterations. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
96 flag = 1; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
97 for iter = 1 : maxit |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
98 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
99 if (nargin > 4 && isnumeric (precon)) |
9278 | 100 z = precon * res; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
101 elseif (nargin > 4) |
8507 | 102 ## Our preconditioner is a function. |
9278 | 103 z = feval (precon, res); |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
104 else |
8507 | 105 ## We don't use preconditioning. |
9278 | 106 z = res; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
107 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
108 |
8507 | 109 ## Cache. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
110 ro_old = ro; |
9278 | 111 ro = res' * z; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
112 if (iter == 1) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
113 p = z; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
114 else |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
115 beta = ro / ro_old; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
116 p = z + beta * p; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
117 endif |
8507 | 118 ## Cache. |
119 q = A * p; | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
120 alpha = ro / (p' * q); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
121 x = x + alpha * p; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
122 |
9278 | 123 res = res - alpha * q; |
124 relres = norm(res) / norm_b; | |
125 resvec = [resvec;relres]; | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
126 |
9278 | 127 if (relres <= tol) |
8507 | 128 ## We reach tolerance tol within maxit iterations. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
129 flag = 0; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
130 break; |
8507 | 131 elseif (resvec (end) == resvec (end - 1)) |
132 ## The method stagnates. | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
133 flag = 3; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
134 break; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
135 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
136 endfor; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
137 |
9278 | 138 if (nargout < 1) |
139 if ( flag == 0 ) | |
140 printf (["cgs converged at iteration %i ", | |
141 "to a solution with relative residual %e\n"],iter,relres); | |
142 elseif (flag == 3) | |
143 printf (["cgs stopped at iteration %i ", | |
144 "without converging to the desired tolerance %e\n", | |
145 "because the method stagnated.\n", | |
146 "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); | |
147 else | |
148 printf (["cgs stopped at iteration %i ", | |
149 "without converging to the desired tolerance %e\n", | |
150 "because the maximum number of iterations was reached.\n", | |
151 "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); | |
152 endif | |
153 endif | |
154 | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
155 endfunction |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
156 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
157 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
158 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
159 %!demo |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
160 %! % Solve system of A*x=b |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
161 %! A=[5 -1 3;-1 2 -2;3 -2 3] |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
162 %! b=[7;-1;4] |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
163 %! [a,b,c,d,e]=cgs(A,b) |