I have settled this question [here][1]; in brief, a [FORTRAN implementation][2] of algorithms for solving this problem have been published in the Collected Algorithms of the ACM.

  [1]: http://math.stackexchange.com/questions/55687/55742#55742
  [2]: http://www.netlib.org/toms/655