//
// TP1: algorithme de gradient à pas constant ou optimal
// cas des fonctionnelles quadratiques
//
function [u, nb_iter, residu] = GPO (A,b,u0,epsilon,max_iter)
u=u0
nb_iter=0;
residu=norm(A*u-b)
while (nb_iter<max_iter)&(residu>epsilon)
  nb_iter=nb_iter+1
  d=A*u-b;rho=d'*d/((A*d)'*d);
  u=u-rho*(A*u-b)
  residu=norm(A*u-b)
end
endfunction
//
function [u, nb_iter, residu] = GPF (A,b,u0,epsilon,max_iter,rho)
u=u0
nb_iter=0;
residu=norm(A*u-b)
while (nb_iter<max_iter)&(residu>epsilon)
  nb_iter=nb_iter+1
  u=u-rho*(A*u-b)
  residu=norm(A*u-b)
end
endfunction

A=zeros(5,5)
for i=1:4
  A(i,i+1)=-1;
  A(i+1,i)=-1;  
end
A(1,1)=3,A(2,2)=12;A(3,3)=24;A(4,4)=48;A(5,5)=96
//
//for i=1:5
// disp(det(A(1:i,1:i)))   
//end
//
//
b=[1;2;3;4;5]
max_iter=20000;
epsilon=1E-5;
u0=zeros(5,1)
v=spec(A);
rho=2/(v(1)+v(5))
[u, nb_iter, residu] = GPF (A,b,u0,epsilon,max_iter,rho)
disp('minimum GPF');disp(u)
disp('precision GPF');disp(residu)
disp('nb iterations GPF');disp(nb_iter)
[u, nb_iter, residu] = GPO (A,b,u0,epsilon,max_iter)
disp('minimum GPO');disp(u)
disp('precision GPO');disp(residu)
disp('nb iterations GPO');disp(nb_iter)
