bksolve = function(A, b)
{
# Solves Ax = b for x, assuming that A is a nonsingular upper 
# triangular matrix.

  n = length(b)
  x = b

# add tests for A(i,i) close to 0!
  x[n] = b[n] / A[n,n]
  for (i in seq(n-1, 1, -1)) {
    x[i] = (b[i] - A[i,(i+1):n] %*% x[(i+1):n]) / A[i,i]
  }
  return (x)
}