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

  n = length(b)
  x = b

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

  return (x)
}