function x =  forwsub(U, b)
    n = length(b);
    x = zeros(size(b));
    for k = 1:n
      x(k) = (b(k) - (U(k,1:k-1) * x(1:k-1)))/U(k,k);
    endfor
end;