function [u,k] = nma_SD( A,f,tol,max_iterations)
%function solves $Au=f$ using the method of steepest descent.

%------------------------------------------------------
%file name: nma_SD.m
% This function solves Au=f using the method of
% steepest descent.
%
% INPUT:
%  A: System matrix,  must be PSD matrix
%  f: RHS, column vector
%  tol: tolerance on the residual ||r||. When the residual (f-Au) becomes
%       smaller than this value, then the iterative process will stop and
%       u is returned.
%  max_iterations: if this number is reached before tolerance is reached
%       then process is stopped and error is thrown.
%
% Output:
%  u: solution to Au=f
%  k: number of iterations used to reach the solution
%
% see nma_generate_A_f_poisson_2D_dirichlet_zero_square.m for an
% example on how to generate A,f for the 2D poisson pde on unit square
%
% EXAMPLE:
% This example shows using this function to solve the 2D Poisson PDE
% problem on unit square
%
% clear all
% myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 );
% n = 31;   % for grid 33x33
% [A,f]=nma_generate_A_f_poisson_2D_dirichlet_zero_square(n,myforce);
% A=full(A);
% tol=10^-2; max_iterations=10000;
% [u,k]=nma_SD( -A,f,tol,max_iterations)
% u=reshape(-u,n,n);
% mesh(u)
% compare with exact solution
% uu=-A\f;
% size(uu)
% uu=reshape(uu,n,n);
% figure
% mesh(uu)
%
% by Nasser M. Abbasi
% Matlab 2010 a
% Free to use at your own risk


if nargin ~= 4
   error('nma_steepest_descent_solver:: number of arguments must be 2');
end

if not(nma_is_spd_matrix(A))
   error('nma_steepest_descent_solver:: A is not SPD');
end

[nRowF,nColF] = size(f(:));

if nColF ~= 1
   error('nma_steepest_descent_solver:: f must be a vector');
end

[n,~] = size(A);
if n ~= nRowF
   error('nma_steepest_descent_solver:: f size not compatible with A size');
end

u = zeros(n,1);
r = f;
k = 0;

while (norm(r,2)>tol) && k<=max_iterations
   fprintf('k=%d, norm=%f\n',k,norm(r,2));
   w     = A*r;
   alpha = (r'*r)/(r'*w);
   u     = u + alpha*r;
   r     = r - alpha*w;
   k     = k + 1;
end

if k>max_iterations
   error('nma_steepest_descent_solver:: no convergence before reaching max_iterations');
end

end