//! fichier PUL.sci:

function [U,L]=uldecomp(A)
//---------------------------------------------------
//- DECOMPOSITION U * L
//- cas A= tridiag(.,.,.);  = UL 
//- avec L=tridiag(l,d,0), U=tridiag(0,1,u)
//- ATTENTION A de type sparse, U et L aussi
//---------------------------------------------------
n=size(A,1);
d=zeros(n,1);
l=zeros(n-1,1);
u=zeros(n-1,1);
d(n)=full(A(n,n));
for i=n-1:-1:1
  l(i)=full(A(i+1,i));
  u(i)=full(A(i,i+1))/d(i+1);
  d(i)=full(A(i,i))-l(i)*u(i); 
end
L=diag(sparse(d))+diag(sparse(l),-1);
U=speye(L)+diag(sparse(u),1);
endfunction

function x=montee(U,y);
//-- Resolution de Ux=y: Remontee
  n=length(y);
  x=zeros(n,1); 
  x(n)=y(n)/U(n,n);
  for k=n-1:-1:1; x(k)=(y(k)-U(k,k+1)*x(k+1))/U(k,k); end;
endfunction

function x=p_descente(L,b,g)
//-- Resolution de Lx-b >=0, x>=g, (Lx-b,x-g)=0,
//-- dans le cas L_{ii}>0 et L triang. superieure:
  n=length(b);
  x=zeros(b);
  x(1)=b(1)/L(1,1); x(1)=max(g(1),x(1));
  for k=2:n
    x(k)=(b(k)-L(k,k-1)*x(k-1))/L(k,k);
    x(k)=max(g(k),x(k));
  end
endfunction