How covariance matrix is computed in generalized linear model?

In R I can compute a covariance matrix of glm by vcov function. I wonder how this computation works?

-
The vcov function is just using results computed in summary.glm, see stats:::vcov.glm and stats:::summary.glm. It appears that summary.glm (lines 26-29) relies on the inverse of the matrix whose Choleski decomposition was given to compute the scaled (by the dispersion parameter) or unscaled VC matrix. –  chl Mar 6 at 10:44