% This is matlab code of algorithm 9.2. Last updated on 6-17-2011. function [x,err,iter,flag] = mlbicgstab(A,x,b,Q,M,max_it,tol,kappa) % input: A: N-by-N matrix % M: N-by-N preconditioner matrix % Q: N-by-n auxiliary matrix with columns q_1,..., q_n % x: initial guess % b: right hand side vector % max_it: maximum number of iterations % tol: error tolerance % kappa: (real number) minimization step controller: % zero, standard minimization % positive, Sleijpen-van der Vorst minimization % output: x: solution computed % err: error norm % iter: number of iterations performed % flag: = 0, solution found to tolerance % = 1, no convergence given max_it % = -1, breakdown % storage: c: 1-by-n matrix % x, r, b, u_t, z: N-by-1 matrices % A, M: N-by-N matrices % Q, G, W: N-by-n matrices N = size(A,2); n = size(Q,2); G = zeros(N,n); % initialize work spaces W = zeros(N,n); c = zeros(1,n); % end initialization iter = 0; flag = 1; bnrm2 = norm(b); if bnrm2 == 0.0 bnrm2 = 1.0; end r = b - A*x; err = norm( r ) / bnrm2; if err < tol flag = 0; return end G(:,1) = M \ r; W(:,1) = A*G(:,1); c(1) = Q(:,1)'*W(:,1); if c(1) == 0 flag = -1; return end e = Q(:,1)'*r; for j = 0:max_it for i = 1:n-1 alpha = e / c(i); x = x + alpha*G(:, i); r = r - alpha*W(:, i); err = norm(r)/bnrm2; iter = iter + 1; if err < tol flag = 0; return end if iter >= max_it return end e = Q(:,i+1)'*r; if j >= 1 beta = -e / c(i+1); W(:,i+1) = r + beta*W(:,i+1); G(:,i+1) = beta*G(:,i+1); for s = i+1 : n-1 beta = -Q(:,s+1)'*W(:,i+1)/c(s+1); W(:,i+1) = W(:,i+1) + beta*W(:,s+1); G(:,i+1) = G(:,i+1) + beta*G(:,s+1); end G(:,i+1) = (M \ W(:,i+1)) - (1/omega)*G(:,i+1); else G(:,i+1) = M \ r; end W(:,i+1) = A*G(:,i+1); for s = 0:i-1 beta = -Q(:,s+1)'*W(:,i+1) / c(s+1); W(:,i+1) = W(:,i+1) + beta*W(:,s+1); G(:,i+1) = G(:,i+1) + beta*G(:,s+1); end c(i+1) = Q(:,i+1)'*W(:,i+1); if c(i+1) == 0 flag = -1; return end end alpha = e / c(n); x = x + alpha*G(:, n); r = r - alpha*W(:,n); err = norm(r)/bnrm2; if err < tol flag = 0; iter = iter + 1; return end u_t = M \ r; z = A*u_t; omega = z'*z; if omega == 0 flag = -1; return end rho = z'*r; omega = rho / omega; if kappa > 0 rho = rho / (norm(z)*norm(r)); abs_rho = abs(rho); if (abs_rho < kappa) & (abs_rho ~= 0) omega = omega*kappa/abs_rho; end end if omega == 0 flag = -1; return end x = x + omega*u_t; r = r - omega*z; err = norm(r)/bnrm2; iter = iter + 1; if err < tol flag = 0; return end if iter >= max_it return end e = Q(:,1)'*r; beta = - e/c(1); W(:,1) = r + beta*W(:,1); G(:,1) = beta*G(:,1); for s = 1:n-1 beta = -Q(:,s+1)'*W(:,1)/c(s+1); W(:,1) = W(:,1) + beta*W(:,s+1); G(:,1) = G(:,1) + beta*G(:,s+1); end G(:,1) = (M \ W(:,1)) - (1/omega)*G(:,1); W(:,1) = A*G(:,1); c(1) = Q(:,1)'*W(:,1); if c(1) == 0 flag = -1; return end end