function [x, error, iter, flag] = mlbicgstab(A, x, b, Q, M, max_it, tol) % mlbicgstab.m solves the linear system Ax=b using the % multiple Lanczos BiCGSTAB Method with preconditioning. % % input A matrix % M preconditioner matrix % Q initial matrix [q_1,...,q_k] with columns q_1,...,q_k % x initial guess vector % b right hand side vector % max_it (integer) maximum number of iterations % tol error tolerance % % output x solution vector % error error norm % iter (integer) number of iterations performed % flag (integer): 0 = solution found to tolerance % 1 = no convergence given max_it % -1 = breakdown n = size(A,2); k = size(Q,2); iter = 0; flag = 1; jmax_it = fix(max_it./k); G = zeros(n,k); W = zeros(n,k); D = zeros(n,k); bnrm2 = norm(b); if bnrm2 == 0.0 bnrm2 = 1.0; end r = b - A*x; G(:,k) = r; error = norm( r ) / bnrm2; if error < tol flag = 0; return end for j = 0:jmax_it if (iter >= max_it) | (flag == -1) | (flag == 0) break end g_tld = M\G(:,k); W(:,k) = A*g_tld; c(k) = Q(:,1)'*W(:,k); if c(k) == 0 flag = -1; break end alpha = Q(:,1)'*r/c(k); u = r - alpha*W(:,k); u_tld = M\u; temp = A*u_tld; rho =temp'*temp; if rho == 0 flag = -1; break end rho = -u'*temp/rho; x = x - rho*u_tld + alpha*g_tld; iter = iter + 1; r = rho*temp + u; error = norm(r)/bnrm2; if error < tol flag = 0; break end for i = 1:k if iter >= max_it break end zd = u; zg = r; zw = zeros(n,1); if j ~= 0 for s = i:k-1 beta = -Q(:,s+1)'*zd/c(s); zd = zd + beta*D(:,s); zg = zg + beta*G(:,s); zw = zw + beta*W(:,s); end end beta = rho*c(k); if beta == 0 flag = -1; break end beta = -Q(:,1)'*(r + rho*zw)/beta; zg = zg + beta*G(:,k); zw = rho*(zw + beta*W(:,k)); zd = r + zw; for s = 1:i-1 beta = -Q(:,s+1)'*zd/c(s); zd = zd + beta*D(:,s); zg = zg + beta*G(:,s); end D(:,i) = zd - u; G(:,i) = zg + zw; if i < k c(i) = Q(:,i+1)'*D(:,i); if c(i) == 0 flag = -1; break end alpha = Q(:,i+1)'*u/c(i); u = u - alpha*D(:,i); g_tld = M\G(:,i); x = x + rho*alpha*g_tld; iter = iter + 1; W(:,i) = A*g_tld; r = r - rho*alpha*W(:,i); error = norm(r)/bnrm2; if error < tol flag = 0; break end end end end