function P = upbinv(alpha,beta)
%UPBINV P = upbinv(alpha,beta)
%Compute the inverse of the upper bidiagonal
%matrix with diagonal alpha and superdiagonal beta

n = length(alpha);
P = zeros(n,n);
alphainv = 1./alpha;
tau = zeros(n-1,1);
for i = 1:n-1
   tau(i) = beta(i)/alpha(i);
end

for i = 1:n
   sign = 1;
   for k = i:n
      prod = 1.0;
      for l = i:k-1
         prod = prod*tau(l);
      end
      P(i,k) = sign*alphainv(k)*prod;
      sign = -sign;
   end
end
