gaussj1( a, b)

function x = gaussj1( a, b)
%gaussj1        rsolution  de ax=b par pivot partiel
%               x = gaussj1( a, b)
%               a matrice du systme
%               b second(s) membres
%               Cette fonction est donne dans un but pdagogique
%               x = a\b est plus efficace dans tous les cas
disp('La fonction GAUSSJ1 est donne dans un but pdagogique');
  disp('la division matricielle gauche est 10 fois plus efficace');
  [n, m] = size(a);
  if n ~= m
    error('a doit tre carre');
  end;
  [nb mb] = size(b);
  if n ~= nb
    error('b doit avoir autant de lignes que a');
  end;
  
  a = [a b];
  for l = 1:n-1
    [apiv lpiv] = max(abs(a(l:n,l)));
    lpiv = lpiv+l-1;
    la = a(lpiv,:);
    a(lpiv,:) = a(l,:);
    a(l,:) = la./la(l);
    t1 = a(l,:);
    t1 = t1(ones(size(l+1:n)),:);
    t2 = a(l+1:n,l);
    t2 = t2(:,ones(size(1:m+mb)));
    a(l+1:n, :) = a(l+1:n, :) - t1.*t2;
  end;
  a(n,:) = a(n,:)./a(n,n);
  
  x = zeros(size(b));
  for l = n:-1:1
    x(l,:) = a(l,m+1:m+mb);
    for j = l+1:n
      x(l,:) = x(l,:)- x(j,:)*a(l,j);
    end;
  end;