//
// resolution de A.u=b par la methode de Gauss 
// avec stratégie de pivot.
//
function u=fct_gauss(A,b)
//
[imax,jmax]=size(A);
//
// début boucle principale
for j=2:imax
// choix du pivôt    
    for i=j:imax
      pivot=A(j-1,j-1);
      kp=j-1;     
/// stratégie de pivôt///////////      
      for k=j-1:imax
        if abs(A(k,j-1))>abs(pivot) then
          pivot=A(k,j-1);
          kp=k;
        end
/////////////////////////////////        
      end
      a3=A(j-1,j-1:imax);
      A(j-1,j-1:imax)=A(kp,j-1:imax);
      A(kp,j-1:imax)=a3;
      b3=b(j-1);
      b(j-1)=b(kp);
      b(kp)=b3;
	m=A(i,j-1)/A(j-1,j-1);
      A(i,:)=A(i,:)-A(j-1,:)*m;
      b(i)=b(i)-m*b(j-1);
    end
// fin boucle principale///////////
end
// remontée//
u=zeros(imax,1);
u(imax)=b(imax)/A(imax,imax);
for j=imax-1:-1:1,
  u(j)=(b(j)-A(j,j+1:imax)*u(j+1:imax))/A(j,j);
end
/// fin remontée///
endfunction
//
tic();
n=16;
A=rand(n,n);  // cas aléatoire
//
for i=1:n   // cas qui ne marche pas (matrice de Hilbert)
    for j=1:n
       A(i,j)=1/(i+j-1) 
    end
 end
b=A*ones(n,1);
u=fct_gauss(A,b);
disp('erreur norme 2');disp(norm(u-ones(n,1)));
disp('solution:');disp(u)
disp('temps de calcul:');disp(toc())

