The issue is in the objective function and is caused by two conventions. First, the weights need to be for the squared error. Second. typically an optional $1/2$ is added to avoid the resulting $2$'s in the gradient. So, you are dealing with the following problem:
\begin{align*}
\min_f \quad& \frac{1}{2} {(Xf - r)}^T W (Xf - r) \\
\text{s.t.} \quad& A^T f = 0.
\end{align*}
[Note: I would rather define the constraint as $Af = b$ but followed your original problem.]
The Lagrangian and its gradient w.r.t $f$ are
\begin{align*}
L(f,\lambda) &= \frac{1}{2} {(Xf - r)}^T W (Xf - r) + \lambda^T A^T f, \\
\nabla_f L(f,\lambda) &= X^T W (X f - r) + A \lambda \\
&= X^T W X f + A \lambda - X^T W r .
\end{align*}
where $\lambda$ is the dual variable for the constraint. Resulting stationarity and primal feasibility conditions you derive are found as
$$
\begin{bmatrix}
X^T W X & A \\
A^T & 0
\end{bmatrix}
\begin{bmatrix}
f \\
\lambda
\end{bmatrix}
=
\begin{bmatrix}
X^T W r \\
0
\end{bmatrix}.
$$
While the scalar $1/2$ is optional and only affects the dual variable, the other convention comes from how the weighted least squares is defined. The standard linear model assumes that the errors have constant variance. In weighted least squares the underlying model has errors with nonconstant variance. Indeed, the weight of an observation is proportional to the reciprocal of the error variance for that observation, $w_i = 1/\sigma_i^2$, thereby cancelling out the nonconstant variance. Therefore, $W$ is defined as
$$
W =
\begin{bmatrix}
w_1 & 0 & \ldots & 0 \\
0 & w_2 & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & w_n
\end{bmatrix}
$$
Finally, it is equivalent to transforming the particular linear model by multiplying it with $W^{1/2}$. So if you defined $\mathcal{W}=W^{1/2}$ you'd be fine, too.
