function [x, error, iter, flag] = bicgstab(A, x, b, M, max_it, tol) % -- Iterative template routine -- % Univ. of Tennessee and Oak Ridge National Laboratory % October 1, 1993 % Details of this algorithm are described in "Templates for the % Solution of Linear Systems: Building Blocks for Iterative % Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra, % Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications, % 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps). % % [x, error, iter, flag] = bicgstab(A, x, b, M, max_it, tol) % % bicgstab.m solves the linear system Ax=b using the % BiConjugate Gradient Stabilized Method with preconditioning. % % input A REAL matrix % x REAL initial guess vector % b REAL right hand side vector % M REAL preconditioner matrix % max_it INTEGER maximum number of iterations % tol REAL error tolerance % % output x REAL solution vector % error REAL error norm % iter INTEGER number of iterations performed % flag INTEGER: 0 = solution found to tolerance % 1 = no convergence given max_it % -1 = breakdown: rho = 0 % -2 = breakdown: omega = 0 iter = 0; % initialization flag = 0; bnrm2 = norm( b ); if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end r = b - A*x; error = norm( r ) / bnrm2; if ( error < tol ) return, end omega = 1.0; r_tld = r; for iter = 1:max_it, % begin iteration rho = ( r_tld'*r ); % direction vector if ( rho == 0.0 ) break, end if ( iter > 1 ), beta = ( rho/rho_1 )*( alpha/omega ); p = r + beta*( p - omega*v ); else p = r; end p_hat = M \ p; v = A*p_hat; alpha = rho / ( r_tld'*v ); s = r - alpha*v; if ( norm(s) < tol ), % early convergence check x = x + alpha*p_hat; resid = norm( s ) / bnrm2; break; end s_hat = M \ s; % stabilizer t = A*s_hat; omega = ( t'*s) / ( t'*t ); x = x + alpha*p_hat + omega*s_hat; % update approximation r = s - omega*t; error = norm( r ) / bnrm2; % check convergence if ( error <= tol ), break, end if ( omega == 0.0 ), break, end rho_1 = rho; end if ( error <= tol | s <= tol ), % converged if ( s <= tol ), error = norm(s) / bnrm2; end flag = 0; elseif ( omega == 0.0 ), % breakdown flag = -2; elseif ( rho == 0.0 ), flag = -1; else % no convergence flag = 1; end % END bicgstab.m