Mercurial > hg > octave-lyh
annotate scripts/sparse/cgs.m @ 11100:cdf940db26a0
provide mxIsFunctionHandle MEX interface function
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Sat, 16 Oct 2010 13:27:18 -0400 |
parents | be55736a0783 |
children | 1740012184f9 |
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 -*- |
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
20 ## @deftypefn {Function File} {} cgs (@var{A}, @var{b}) |
8338
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. |
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
26 ## The @var{maxit} specifies the maximum number of iteration, default value is |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
27 ## MIN(20,N). |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
28 ## The @var{M1} specifies a preconditioner, can also be a function handler which |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
29 ## returns M\X. |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
30 ## The @var{M2} combined with @var{M1} defines preconditioner as |
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9279
diff
changeset
|
31 ## preconditioner=M1*M2. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
32 ## 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
|
33 ## |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
34 ## @end deftypefn |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
35 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
36 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
|
37 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
38 if (nargin < 2 || nargin > 7 || nargout > 5) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
39 print_usage (); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
40 elseif (!isnumeric (A) || rows (A) != columns (A)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
41 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
|
42 elseif (!isvector (b)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
43 error ("cgs: b must be a vector"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
44 elseif (rows (A) != rows (b)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
45 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
|
46 elseif (nargin > 2 && !isscalar (tol)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
47 error ("cgs: tol must be a scalar"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
48 elseif (nargin > 3 && !isscalar (maxit)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
49 error ("cgs: maxit must be a scalar"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
50 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
|
51 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
|
52 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
|
53 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
|
54 elseif (nargin > 6 && !isvector (x0)) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
55 error ("cgs: x0 must be a vector"); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
56 elseif (nargin > 6 && rows (x0) != rows (b)) |
8507 | 57 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
|
58 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
59 |
8507 | 60 ## Default tolerance. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
61 if (nargin < 3) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
62 tol = 1e-6; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
63 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
64 |
8507 | 65 ## Default maximum number of iteration. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
66 if (nargin < 4) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
67 maxit = min (rows (b),20); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
68 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
69 |
8507 | 70 ## Left preconditioner. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
71 if (nargin == 5) |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
72 if (isnumeric (M1)) |
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
73 precon = @(x) M1 \ x; |
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
74 endif |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
75 elseif (nargin > 5) |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
76 if (issparse (M1) && issparse (M2)) |
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
77 precon = @(x) M2 \ (M1 \ x); |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
78 else |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
79 M = M1*M2; |
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
80 precon = @(x) M \ x; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
81 endif |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
82 else |
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
83 precon = @(x) x; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
84 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
85 |
8507 | 86 ## Specifies initial estimate x0. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
87 if (nargin < 7) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
88 x = zeros (rows (b), 1); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
89 else |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
90 x = x0; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
91 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
92 |
9278 | 93 res = b - A * x; |
94 norm_b = norm (b); | |
8507 | 95 ## Vector of the residual norms for each iteration. |
9278 | 96 resvec = [ norm(res)/norm_b ]; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
97 ro = 0; |
8507 | 98 ## 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
|
99 flag = 1; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
100 for iter = 1 : maxit |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
101 |
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
102 z = precon (res); |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
103 |
8507 | 104 ## Cache. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
105 ro_old = ro; |
9278 | 106 ro = res' * z; |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
107 if (iter == 1) |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
108 p = z; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
109 else |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
110 beta = ro / ro_old; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
111 p = z + beta * p; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
112 endif |
8507 | 113 ## Cache. |
114 q = A * p; | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
115 alpha = ro / (p' * q); |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
116 x = x + alpha * p; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
117 |
9278 | 118 res = res - alpha * q; |
119 relres = norm(res) / norm_b; | |
120 resvec = [resvec;relres]; | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
121 |
9278 | 122 if (relres <= tol) |
8507 | 123 ## We reach tolerance tol within maxit iterations. |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
124 flag = 0; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
125 break; |
8507 | 126 elseif (resvec (end) == resvec (end - 1)) |
127 ## The method stagnates. | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
128 flag = 3; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
129 break; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
130 endif |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
131 endfor; |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
132 |
9278 | 133 if (nargout < 1) |
134 if ( flag == 0 ) | |
135 printf (["cgs converged at iteration %i ", | |
136 "to a solution with relative residual %e\n"],iter,relres); | |
137 elseif (flag == 3) | |
138 printf (["cgs stopped at iteration %i ", | |
139 "without converging to the desired tolerance %e\n", | |
140 "because the method stagnated.\n", | |
141 "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); | |
142 else | |
143 printf (["cgs stopped at iteration %i ", | |
144 "without converging to the desired tolerance %e\n", | |
145 "because the maximum number of iterations was reached.\n", | |
146 "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); | |
147 endif | |
148 endif | |
149 | |
8338
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
150 endfunction |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
151 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
152 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
153 |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
154 %!demo |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
155 %! % Solve system of A*x=b |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
156 %! 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
|
157 %! b=[7;-1;4] |
a35bf360b919
Add the cgs and treelayout functions
Radek Salac <salac.r@gmail.com>
parents:
diff
changeset
|
158 %! [a,b,c,d,e]=cgs(A,b) |