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