function x=redBlackGSeidel(A,f,x,M,tau); % GSEIDEL: Solves Ax=f using Red-Black Gauss-Seidel iteration. % USAGE: % x=redBlackGSeidel(A,f,x,M,tau); % % INPUTS: % A : nxn matrix % f : nxm matrix % [x] : nxm matrix of starting values. Default: f % [M] : maximum number of iterations % [tau] : convergence tolerence % % OUTPUT: % x : nxm matrix of approximate solutions to Ax=f % % Written by: % -- % John L. Weatherwax 2006-12-19 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- if nargin<3 | isempty(x), x=f; end if nargin<4 | isempty(M), M=1000; end if nargin<5 | isempty(tau), tau=sqrt(eps); end % obtain the factors D, L and U in the splitting A = D - L - U: D=diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); I = eye(size(A)); Z = zeros(size(f)); % obtain the coefficient matrix when processing the % EVEN unknowns (the identity and zero avoid touching the ODD unknowns): PG_E = (D-L) \ U; PG_E(1:2:end,:) = I(1:2:end,:); fG_E = (D-L) \ f; fG_E(1:2:end,:) = Z(1:2:end,:); % obtain the coefficient matrix when processing the % ODD unknowns (the identity and zero avoid touching the EVEN unknowns): PG_O = (D-L) \ U; PG_O(2:2:end,:) = I(2:2:end,:); fG_O = (D-L) \ f; fG_O(2:2:end,:) = Z(2:2:end,:); divtol=100*max(abs(f-A*x)); for i=1:M % first update the even components x_new = PG_E * x + fG_E; % second update the odd components x_new = PG_O * x_new + fG_O; r=f-A*x_new; % the new residual dx=x_new - x; % convergence check if( all(abs(r)divtol ) disp('WARNING: Iterations diverging'); return end x = x_new; end if i==M, disp('WARNING: Maximum iterations exceeded'); end