function y = backsolve(tmat,x,upper,transpose)
% BACKSOLVE(TMAT,X,UPPER,TRANSPOSE) backsolve for triangular matrix
% BACKSOLVE(TMAT,X,UPPER,TRANSPOSE) returns tmat^{-1} x if transpose = 1
% BACKSOLVE(TMAT,X,UPPER,TRANSPOSE) returns tmat^{-T} x if transpose = 0
% TMAT is assumed upper triangular if upper = 1
% else lower triangular.
if transpose==1
   if upper==1
      y = flipud(backsolve0(flipud(fliplr(tmat')),flipud(x)));
   else
      y = backsolve0(tmat',x);
   end
else
   if upper==1
      y = backsolve0(tmat,x);
   else
      y = flipud(backsolve0(flipud(fliplr(tmat)),flipud(x)));
   end
end


function y = backsolve0(tmat,x)
% backsolve for upper triangular matrix
[n m ] = size(tmat);
y = zeros(size(x));
y(n,:) = x(n,:)/tmat(n,n);
for i=(n-1):(-1):1
   y(i,:) = x(i,:) - tmat(i,(i+1):n)*y((i+1):n,:);
   y(i,:) = y(i,:)/tmat(i,i);
end

