function [P,dP,Z] = gseid(A,B,P,delta,max1)
% [X,dX] = gseid(A,B,P,delta,max1)
% [X,dX,Z] = gseid(A,B,P,delta,max1)
% A is an nxn diagonally dominant matrix, input.
% B is a nx1 vector, input.
% P is an nX1 starting vector, input.
% The tolerance is delta, input.
% The maximum number of iterations is max1, input.
% X is an nx1 vector, output.
% dX is an nx1 vector, output.
% Z is a list of iterations, output.
Z = P';
n = length(B);
Pold = P;
for k=1:max1,
for r = 1:n,
Sum1 = B(r) - A(r,[1:r-1,r+1:n])*P([1:r-1,r+1:n]);
P(r) = Sum1/A(r,r);
end
dP = abs(Pold-P);
err = norm(dP);
relerr = err/(norm(P)+eps);
Pold = P;
Z = [Z;P'];
if (err<delta)|(relerr<delta), break, end
end