How to calculate the regression variance for a GLS model? I need to calculate the regression variance ($\sigma^2$) in order to estimate both the confidence intervals and the prediction intervals in a gls regression analysis.  For the analysis, the covariance matrix ($V$) of the response variable ($y$) is known in advance, and so I use it directly as the weighting matrix (=$V^{-1}$) in the gls regression analysis.
The regression variance is a weighted sum of the residual error:
$\sigma^2 = \frac{ (Y – X\beta)^T C^{-1} (Y – X\beta)}{n – p}$
My question/problem is how to determine the weighting matrix $C^{-1}$?  $C$ cannot be set equal to $V$ since (according to the above equation) $C$ must be dimensionless while $V$ has the same units as $\sigma^2$.
Based on my reading of the literature and available texts, it seems that $C$ is the correlation matrix and is a scaled or normalized form of the covariance matrix $V$.  i.e., $V = Var(\epsilon^2) = \sigma^2 C$.  But my problem is that $\sigma^2$ is not yet known, and so I need another way find $C$ from $V$.
R functions such as gls() will compute the regression variance (if I knew how gls() does this, it would answer my question).  However I cannot use gls() in this case since I am specifying a user-defined covariance (weighting) matrix, and gls() only accepts a limited set of specific correlation structures.
In fact a possible solution can be found in this earlier post where an equation for the SEE (or sigma2) for a GLS regression was cited :
GLS calc of SEE: sqrt( sum( ( residuals from linear model) ^ 2 * glsWeight ) ) / sum( glsWeight ) * length( glsWeight ) / residualDegreeFreedom )
However I am unable to ascertain the validity of this equation and cannot find its source reference.
 A: You are a bit off I think because the formula for epsilon you seem to use relies on the neccesity to decompose your error variance into a constant term ($\sigma^2$) and a correlation term $\Omega$, which is the $C$ in your post. You actually have to use this $C$ as a weight for the regression. Let me explain.  
It is assumed you know $\Omega$, otherwise you have to do FGLS instead. The formula you have is a result of a transformed model. First, let's look at our GLS estimators:
Under the assumption that $V[\epsilon|X]=\sigma^2\Omega$ you can transform the model and estimate the following
$$\widehat{\beta}^{GLS} = [X'\Omega^{-1}X]^{-1}X'\Omega^{-1}y$$
This, generally, is a BLUE estimator under the assumptions.
Which has, under the above and usual assumptions this variance
$$V[\widehat{\beta}^{GLS}] = \sigma^2[X'\Omega^{-1}X]^{-1} $$
Since you know $\Omega$, you can now calculate everything and estimate the $\sigma^2$. To understand how have a look at how the GLS estimator is reached.
You reached this result by applying a matrix $P$ to the model formula as such
$$PY = PX\beta + P\epsilon$$
What is P? If you know omega, you can calculate the $P$ (Cholesky decomposition) with
$$ P'P = \Omega^{-1}$$
And then you have a consistent estimator in
$$ \widehat{\sigma}^2 = \frac{(P\hat{\epsilon}^{OLS})'P\hat{\epsilon}^{OLS}}{n-k} \\
=  
\frac{ \hat{\epsilon}'^{OLS} \Omega^{-1} \hat{\epsilon}^{OLS}}{n-k} $$
This is the formula you have listed, but it implies a possible decomposition  $V[\epsilon|X]=\sigma^2\Omega$. Note how it uses the residuals of OLS. This is because of said existence of a transformed and non-transformed model.  
You said you used $V[Y]=V$ as weight for GLS. Can you tell us the exact formula? Is it
$$\widehat{\beta}^{GLS} = [X'V^{-1}X]^{-1}X'V^{-1}y$$
?
In that case you used $V$ instead of $C$. You mentioned that you think $V[\epsilon]=V$ and $V=\sigma^2C$ (I assume this is what you mean).
You could go ahead and plug in $V$ instead of $C$.
But of course this is not actually correct because $V$ is not equal to $C$ as you discovered.
This seems to be your problem.
If you know nothing about the composition of your error terms, you just can't get to that decomposition and the $C$ you need. That is the 'point' of GLS. It takes advantage of a known correlation structure $\Omega$ or $C$ to regain BLUE estimators.
So the equation you can not validate is the result of said transformed model. For this you also need to have your $\Omega$ or $C$ as you called it and by that the $P$ to transform the model.
If you do not have this, use FGLS instead.
