Mercurial > hg > octave-nkf
annotate scripts/sparse/bicgstab.m @ 9665:1dba57e9d08d
use blas_trans_type for xgemm
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sat, 26 Sep 2009 10:41:07 +0200 |
| parents | 1673a0dc019f |
| children | be55736a0783 |
| 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 -*- | |
| 20 ## @deftypefn {Function File} {} bicgstab (@var{A}, @var{b}) | |
| 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. | |
| 25 ## The @var{tol} specifies the tolerance of the method, the default value is 1e-6. | |
| 26 ## The @var{maxit} specifies the maximum number of iterations, the default value is min(20,N). | |
| 27 ## The @var{M1} specifies a preconditioner, can also be a function handler which returns M\X. | |
| 28 ## The @var{M2} combined with @var{M1} defines preconditioner as preconditioner=M1*M2. | |
| 29 ## The @var{x0} is the initial guess, the default value is zeros(N,1). | |
| 30 ## | |
| 31 ## The value @var{x} is a computed result of this procedure. | |
| 32 ## The value @var{flag} can be 0 when we reach tolerance in @var{maxit} iterations, 1 when | |
| 33 ## we don't reach tolerance in @var{maxit} iterations and 3 when the procedure stagnates. | |
| 34 ## The value @var{relres} is a relative residual - norm(b-A*x)/norm(b). | |
| 35 ## The value @var{iter} is an iteration number in which x was computed. | |
| 36 ## The value @var{resvec} is a vector of @var{relres} for each iteration. | |
| 37 ## | |
| 38 ## @end deftypefn | |
| 39 | |
| 40 function [x, flag, relres, iter, resvec] = bicgstab (A, b, tol, maxit, M1, M2, x0) | |
| 41 | |
| 42 if (nargin < 2 || nargin > 7 || nargout > 5) | |
| 43 print_usage (); | |
| 44 elseif (!isnumeric (A) || rows (A) != columns (A)) | |
| 45 error ("bicgstab: the first argument must be a n-by-n matrix"); | |
| 46 elseif (!isvector (b)) | |
| 47 error ("bicgstab: b must be a vector"); | |
| 48 elseif (!any (b)) | |
| 49 error ("bicgstab: b shuldn't be a vector of zeros"); | |
| 50 elseif (rows (A) != rows (b)) | |
| 51 error ("bicgstab: the first and second argument must have the same number of rows"); | |
| 52 elseif (nargin > 2 && !isscalar (tol)) | |
| 53 error ("bicgstab: tol must be a scalar"); | |
| 54 elseif (nargin > 3 && !isscalar (maxit)) | |
| 55 error ("bicgstab: maxit must be a scalar"); | |
| 56 elseif (nargin > 4 && ismatrix (M1) && (rows (M1) != rows (A) || columns (M1) != columns (A))) | |
| 57 error ("bicgstab: M1 must have the same number of rows and columns as A"); | |
| 58 elseif (nargin > 5 && (!ismatrix (M2) || rows (M2) != rows (A) || columns (M2) != columns (A))) | |
| 59 error ("bicgstab: M2 must have the same number of rows and columns as A"); | |
| 60 elseif (nargin > 6 && !isvector (x0)) | |
| 61 error ("bicgstab: x0 must be a vector"); | |
| 62 elseif (nargin > 6 && rows (x0) != rows (b)) | |
| 8507 | 63 error ("bicgstab: x0 must have the same number of rows as b"); |
| 8416 | 64 endif |
| 65 | |
| 8507 | 66 ## Default tolerance. |
| 8416 | 67 if (nargin < 3) |
| 68 tol = 1e-6; | |
| 69 endif | |
| 70 | |
| 8507 | 71 ## Default maximum number of iteration. |
| 8416 | 72 if (nargin < 4) |
| 73 maxit = min (rows (b), 20); | |
| 74 endif | |
| 75 | |
| 8507 | 76 ## Left preconditioner. |
| 8416 | 77 if (nargin == 5) |
|
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
78 if (isnumeric (M1)) |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
79 precon = @(x) M1 \ x; |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
80 endif |
| 8416 | 81 elseif (nargin > 5) |
|
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
82 if (issparse (M1) && issparse (M2)) |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
83 precon = @(x) M2 \ (M1 \ x); |
| 8416 | 84 else |
|
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
85 M = M1*M2; |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
86 precon = @(x) M \ x; |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
87 endif |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
88 else |
|
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
89 precon = @(x) x; |
| 8416 | 90 endif |
| 91 | |
| 92 ## specifies initial estimate x0 | |
| 93 if (nargin < 7) | |
| 94 x = zeros (rows (b), 1); | |
| 95 else | |
| 96 x = x0; | |
| 97 endif | |
| 98 | |
| 99 norm_b = norm (b); | |
| 100 | |
| 101 res = b - A*x; | |
| 102 rr = res; | |
| 103 | |
| 8507 | 104 ## Vector of the residual norms for each iteration. |
| 9278 | 105 resvec = [norm(res)/norm_b]; |
| 8416 | 106 |
| 8507 | 107 ## Default behaviour we don't reach tolerance tol within maxit iterations. |
| 8416 | 108 flag = 1; |
| 109 | |
| 110 for iter = 1:maxit | |
| 111 rho_1 = res' * rr; | |
| 112 | |
| 113 if (iter == 1) | |
| 114 p = res; | |
| 115 else | |
| 116 beta = (rho_1 / rho_2) * (alpha / omega); | |
| 117 p = res + beta * (p - omega * v); | |
| 118 endif | |
| 119 | |
|
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
120 phat = precon (p); |
| 8416 | 121 |
| 122 v = A * phat; | |
| 123 alpha = rho_1 / (rr' * v); | |
| 124 s = res - alpha * v; | |
| 125 | |
|
9279
1673a0dc019f
fix & improve preconditioning strategy in cgs & bicgstab
Jaroslav Hajek <highegg@gmail.com>
parents:
9278
diff
changeset
|
126 shat = precon (s); |
| 8416 | 127 |
| 128 t = A * shat; | |
| 129 omega = (t' * s) / (t' * t); | |
| 130 x = x + alpha * phat + omega * shat; | |
| 131 res = s - omega * t; | |
| 132 rho_2 = rho_1; | |
| 133 | |
| 134 relres = norm (res) / norm_b; | |
| 135 resvec = [resvec; relres]; | |
| 136 | |
| 137 if (relres <= tol) | |
| 8507 | 138 ## We reach tolerance tol within maxit iterations. |
| 8416 | 139 flag = 0; |
| 140 break; | |
| 8507 | 141 elseif (resvec (end) == resvec (end - 1)) |
| 142 ## The method stagnates. | |
| 8416 | 143 flag = 3; |
| 144 break; | |
| 145 endif | |
| 146 endfor | |
| 147 | |
| 9278 | 148 if (nargout < 2) |
| 149 if (flag == 0) | |
| 150 printf (["bicgstab converged at iteration %i ", | |
| 151 "to a solution with relative residual %e\n"],iter,relres); | |
| 152 elseif (flag == 3) | |
| 153 printf (["bicgstab stopped at iteration %i ", | |
| 154 "without converging to the desired tolerance %e\n", | |
| 155 "because the method stagnated.\n", | |
| 156 "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); | |
| 157 else | |
| 158 printf (["bicgstab stopped at iteration %i ", | |
| 159 "without converging to the desired toleranc %e\n", | |
| 160 "because the maximum number of iterations was reached.\n", | |
| 161 "The iterate returned (number %i) has relative residual %e\n"],iter,tol,iter,relres); | |
| 162 endif | |
| 163 endif | |
| 164 | |
| 8416 | 165 endfunction |
| 166 | |
| 167 %!demo | |
| 168 %! % Solve system of A*x=b | |
| 169 %! A = [5 -1 3;-1 2 -2;3 -2 3] | |
| 170 %! b = [7;-1;4] | |
| 171 %! [x, flag, relres, iter, resvec] = bicgstab(A, b) | |
| 172 |