//
// programme 5
//
write(%io(2),'programme 5: calcul de l''inverse d''une matrice')
write(%io(2),'par la methode de Gauss')
write(%io(2),'//')
//
write(%io(2),'taille de la matrice (aleatoire) a inverser?')
imax=read(%io(1),1,1);
//
// tirage de la matrice a inverser
//
a=rand(imax,imax);
//
timer();
ainv=zeros(imax,imax);
for l=1:imax
//
// second membre
//
   b=zeros(imax,1);
   b(l)=1;
   a2=a;
//
   for j=2:imax
    for i=j:imax
      pivot=a2(j-1,j-1);
      kp=j-1;
      for k=j-1:imax 				// pivot partiel 
        if abs(a2(k,j-1))>abs(pivot) then
          pivot=a2(k,j-1);
          kp=k;
        end
      end
      a3=a2(j-1,j-1:imax);
      a2(j-1,j-1:imax)=a2(kp,j-1:imax);
      a2(kp,j-1:imax)=a3;
      b3=b(j-1);
      b(j-1)=b(kp);
      b(kp)=b3;
	m=a2(i,j-1)/a2(j-1,j-1);
      a2(i,:)=a2(i,:)-a2(j-1,:)*m;
      b(i)=b(i)-m*b(j-1);
    end
   end
//
// resolution du systeme lineaire
//
   x=zeros(imax,1);
   x(imax)=b(imax)/a2(imax,imax);
   for j=imax-1:-1:1,
     x(j)=(b(j)-a2(j,j+1:imax)*x(j+1:imax))/a2(j,j);
   end
//
   ainv(:,l)=x;
//
end
//
// ecriture des resultats
//
write(%io(2),'inverse1 (pivot de Gauss):'),disp(ainv)
t1=timer();
write(%io(2),'temps de calcul:')
disp(t1)
write(%io(2),'inverse2 (instruction inv):'),disp(inv(a))
//comparaison avec l'instruction inv
t2=timer();
write(%io(2),'temps de calcul:')
disp(t2)
write(%io(2),'erreur relative:')
disp(norm(ainv-inv(a))/norm(ainv))



