% MATLAB script to experiment with descent methods for solving Ax=b
%
% Written by: Douglas B. Meade
% Created: 17 April 2006

% Define A, b, and initial guess (x)

A = [ 5 2 1 1; 2 6 2 1; 1 2 7 2; 1 1 2 8 ];
b = [ 29 31 26 19 ]';
x = zeros(n,1);

[m,n] = size(A);

M = eye(n,n);         % no preconditioning
%M = diag(diag(A));    % Jacobi preconditioner
%M = tril(A);          % Gauss-Seidel preconditioner

E = eye(n,n);

L = 100;
eps = 1e-4;

% Loop to compute next ten iterates

r = b-A*x;
k = 1;
while norm(r,2)>eps,
   %p = E(:,mod(k-1,n)+1);    % cycle through the standard coordinate unit vectors
   %p = r;                    % steepest descent
   p = inv(M)*r;             % preconditioned steepest descent
   x = descent1(A,b,x,p);
   [k norm(r,2)], [p x]
%   disp('Press any key to continue')
%   pause                             % hit any key to continue
   r = b-A*x;
   if k==L,
      break;
   end
   k = k+1;
end