% This is matlab code of preconditioned mlbicgstabt (see Appendix section). % Last updated on 6-13-2011. function [x,err,iter,flag] = mlbicgstabt(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 shadow matrix [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: % kappa = 0, standard minimization % kappa > 0, 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 iterations % = -1, breakdown % storage: F: N-by-(n-1) matrix. % G, Q, W: N-by-n matrices. % A, M: N-by-N matrices. % x, r, g_h, z, b: N-by-1 matrices. % c: 1-by-n matrix. N = size(A,2); n = size(Q,2); G = zeros(N,n); % initialize work spaces W = zeros(N,n); if n > 1 F = zeros(N,n-1); end 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 if n > 1 F = M' \ (A' * Q(:,1: n-1)); end G(:,1) = r; g_h = M \ r; W(:,1) = A*g_h; 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_h; 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) = W(:,i+1) - G(:,i+1)./omega; for s = 0:i-1 beta = -F(:,s+1)'*G(:,i+1) / c(s+1); G(:,i+1) = G(:,i+1) + beta*G(:,s+1); end else beta = -F(:,1)'*r / c(1); G(:,i+1) = r + beta*G(:,1); for s = 1:i-1 beta = -F(:,s+1)'*G(:,i+1) / c(s+1); G(:,i+1) = G(:,i+1) + beta*G(:,s+1); end end g_h = M \ G(:,i+1); W(:,i+1) = A*g_h; 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_h; r = r - alpha*W(:,n); err = norm(r)/bnrm2; if err < tol flag = 0; iter = iter + 1; return end g_h = M \ r; z = A*g_h; 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_om = abs(rho); if (abs_om < kappa) & (abs_om ~= 0) omega = omega*kappa/abs_om; end end if omega == 0 flag = -1; return end x = x + omega*g_h; 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) = W(:,1) - G(:,1)./omega; g_h = M \ G(:,1); W(:,1) = A * g_h; c(1) = Q(:,1)'*W(:,1); if c(1) == 0 flag = -1; return end end