gj( a, b)

function x = gj( a, b)
%gaussjordan    pivot partiel
%               x = gaussjordan( a, b)
%               version dbutant
%               a matrice du systme
%               b second(s) membre(s)
%               Cette fonction est donne dans un but pdagogique
%               x = a\b est plus efficace dans tous les cas
  disp('La fonction gj 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                         % mise sous forme triangulaire
    [apiv lpiv] = max(abs(a(l:n,l)));
    lpiv = lpiv+l-1;
    a([l lpiv],:) = a([lpiv,l],:);    % change de deux lignes
    a(l,:) = a(l,:)./a(lpiv,l);       % divise  (*)
    for j = l+1:n                     % soustrait multiplie (**)
      a(j,:) = a(j,:) - a(j,l) .*  a(l,:);
    end;
  end;
  x = zeros(size(b));                 % rsolution du systme triangulaire
  for i = n:-1:1
    x(i,:) = a(i,n+1:n+mb);
    for j = i+1:n
      x(i,:) = x(i,:)- x(j,:)*a(i,j); % (***)
    end;
  end;