can BGFS be made to return the approximate hessian?

["Steven Scott" <steve_the_bayesian at yahoo dot com>]

In statistical applications, the hessian matrix at the mode of the
likelihood function is of interest because -H^{-1} is is the asymptotic
variance of the ML estimator.

For nonlinear statistical models the BGFS algorithm is often used to
find the mode of the likelihood function (and thus the MLE point
estimate).    The BGFS algorithm constructs an approximate hessian
matrix internally, but the GSL implementation of BGFS does not seem to
expose the approximate hessian to users.  

Is there an alternate interface to BFGS that provides an estimate of the
Hessian at the mode?  If not, that would be a valuable (and seemingly
easy) extension in a future GSL release.

Cheers, 
Steve Scott