function x=symmetricGSeidel(A,f,x,M,tau); % GSEIDEL: Solves Ax=f using the symmetric Gauss-Seidel iteration. % USAGE: % x=symmetricGSeidel(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); % obtain the coefficient matrix (and right hand side) when sweeping in % ASSENDING order of the unknowns: PG_A = (D-L) \ U; f_A = (D-L) \ f; % obtain the coefficient matrix (and right hand side) when sweeping in % DECENDING order of the unknowns: PG_D = (D-U) \ L; f_D = (D-U) \ f; divtol=100*max(abs(f-A*x)); for i=1:M if( mod(i,2)~=0 ) % sweep in ascending order (D-L) x^{n+1} = U x^n + f x_new = PG_A * x + f_A; else % sweep in descending order (D-U) x^{n+1} = L x^n + f x_new = PG_D * x + f_D; end 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