function [x,kk,err]=grad(A,b,x0,kmax,tol)
%
% Gradient method for linear systems
%
% Input
%  A     Matrix (square matrix)
%  b     r.h.s (column vector)
%  x0    initial guess (column vector)
%  kmax  maximum number of iterations
%  tol  Tolerance for the stopping test
%
% Output
%  x    solution (column vector)
%  kk  Iterations
%  err   array with errors at each iteration
%
%    [x,kk,err]=grad(A,b,x0,kmax,tol)
%
x=x0;
[n,m]=size(A);
r=b-A*x;
bb=norm(b,2);
k=0;err1=norm(r)/bb;
err=[err1];
while k<kmax & err1> tol
q=A*r;
nrk=r'*r;
alpha=nrk/(r'*q);
%
x=x+alpha*r;
%
r=r-alpha*q;
res=norm(r,2);
err1=res/bb;
err=[err;err1];
k=k+1;
end
kk=k;
fprintf('\n The Gradient method has satisfied the stopping test in')
fprintf('%3g iterations,\n with tolerance = %6.3g\n', kk, tol)

return

