function [x,err,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 N-by-N matrix % M N-by-N preconditioner matrix % Q N-by-n auxiliary matrix [q_1,...,q_n] with columns % q_1,...,q_n % x initial guess % b right hand side vector % max_it (integer) maximum number of iterations % tol error tolerance % % output x solution % err error norm % iter (integer) number of iterations performed % flag (integer): 0 = solution found to tolerance % 1 = no convergence given max_it % -1 = breakdown % % storage: % alphaa, betaa, rhoo, bnrm2, err, tol, d, e, f --- real/complex numbers % i, j, s, n, N, flag, iter, max_it --- integers % c --- real/complex array of length n % x, r, g, g_tld, u, u_tld, w, z, b --- vectors of length N % A, M --- N-by-N matrices % Q, G, W --- N-by-n matrices % D --- N-by-(n-1) matrix. % D is not required when n = 1. % N = size(A,2); n = size(Q,2); G = zeros(N,n); % reserve spaces for G, W, D and c W = zeros(N,n); if n > 1 D = zeros(N,n-1); end c = zeros(1,n); % end reservation 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(:,n) = r; g_tld = M\r; W(:,n) = A*g_tld; c(n) = Q(:,1)'*W(:,n); if c(n) == 0 flag = -1; return end e = Q(:,1)'*r; for j = 0:max_it alphaa = e/c(n); u = r - alphaa*W(:,n); u_tld = M\u; z = A*u_tld; rhoo = z'*z; if rhoo == 0 flag = -1; return end rhoo = - z'*u / rhoo; x = x - rhoo*u_tld + alphaa*g_tld; r = rhoo*z + u; err = norm(r)/bnrm2; iter = iter + 1; if err < tol flag = 0; return end if iter >= max_it return end d = rhoo*c(n); if d == 0 flag = -1; return end for i = 1:n-1 z = u; g = zeros(N,1); w = zeros(N,1); f = Q(:,i+1)'*u; if j >= 1 betaa = - f / c(i); z = z + betaa*D(:,i); g = g + betaa*G(:,i); w = w + betaa*W(:,i); for s = i+1:n-1 betaa = -Q(:,s+1)'*z / c(s); z = z + betaa*D(:,s); g = g + betaa*G(:,s); w = w + betaa*W(:,s); end end w = r + rhoo*w; betaa = -Q(:,1)'*w/ d; w = w + (rhoo*betaa)*W(:,n); g = g + w + betaa*G(:,n); for s = 1:i-1 betaa = -Q(:,s+1)'*w / c(s); g = g + betaa*G(:,s); w = w + betaa*D(:,s); end D(:, i) = w - u; c(i) = Q(:,i+1)'*D(:,i); if c(i) == 0 flag = -1; return end G(:,i) = g; g_tld = M\g; W(:,i) = A*g_tld; alphaa = f/c(i); if i < n - 1 u = u - alphaa*D(:,i); end alphaa = rhoo*alphaa; x = x + alphaa*g_tld; r = r - alphaa*W(:,i); err = norm(r)/bnrm2; iter = iter + 1; if err < tol flag = 0; return end if iter >= max_it return end end e = Q(:,1)'*r; betaa = - e / d; w = r + (rhoo*betaa)*W(:,n); g = w + betaa*G(:, n); for s = 1 : n-1 betaa = - Q(:,s+1)'*w / c(s); g = g + betaa*G(:,s); w = w + betaa*D(:, s); end G(:, n) = g; g_tld = M \ g; W(:, n) = A*g_tld; c(n) = Q(:,1)'*W(:,n); if c(n) == 0 flag = -1; return end end