I am interested in doing general linear least squares (NR3, 15.4) using QR decomposition. NR mentions this possiblity, but only covers the normal equation and SVD methods in detail.
I have no problem QR decomposing the design matrix and getting the fit parameters. Since the covariance matrix C = (R^T R)^(-1), I can get the covariance matrix by forming R^T R and inverting it explicitly. But Golub and van Loan mention in their problem 5.3.6 that there is a 2/3 n^3 method for overwriting R directly with C. I'm not smart enough to solve their problem by deriving the algorithm. Can someone point out the algorithm in the literature, a textbook, or as code available for perusal?