This type of model is actually much more common in certain branches of physics and science than "normal" linear regression. So, in physics tools like ROOT, doing this type of fit is trivial, while linear regression is not natively implemented! Physicists tend to call this just a "fit" or a chi-square minimizing fit.
The normal linear regression model assumes that there is an overall variance $\sigma$ attached to every measurement. It then maximizes the likelihood $$ L \propto \prod_i e^{-\frac{1}{2} \left( \frac{y_i-(ax_i+b)}{\sigma} \right)^2} $$ or equivalently its logarithm $$ \log(L) = \mathrm{constant} - \frac{1}{2\sigma^2} \sum_i (y_i-(ax_i+b))^2 $$ Hence the name least-squares -- maximizing the likelihood is the wave as minimizing the sum of squares, and $\sigma$ is an unimportant constant, as long as it is constant. With measurements that have different known uncertainties, you'll want to maximize $$ L \propto \prod e^{-\frac{1}{2} \left( \frac{y-(ax+b)}{\sigma_i} \right)^2} $$ or equivalently its logarithm $$ \log(L) = \mathrm{constant} - \frac{1}{2} \sum \left( \frac{y_i-(ax_i+b)}{\sigma_i} \right)^2 $$ So, you actually want to weight the measurements by the inverse variance $1/\sigma_i^2$, not the variance. This makes sense -- a more accurate measurement has smaller uncertainty and should be given more weight. Note that if this weight is constant, it still factors out of the sum. So, it doesn't affect the estimated values, but it should affect the standard errors, taken from thee second derivative of $\log(L)$.
However, here we come to a another difference between physics and statistics at large. lm's interpretation of the uncertainties attached to your measurements and of the weights isn't quite the same as what you want. It'll take the weights, but it still thinks there is an overall $\sigma^2$ to account for regression error, which is not what you want -- you want your measurement errors to be the only kind of error there is. The question and answer here have more details:
lm weights and the standard error
There are a couple of possible solutions given in the answers there. In particular, an anonymous answer there suggests
vcov(mod)/summary(mod)$sigma^2
Basically, lm scales the covariance matrix based on its estimated $\sigma$, and you want to undo this. You can then get the information you want from the corrected covariance matrix. Try this, but try to double-check it if you can with manual linear algebra. And remember that the weights should the the inverse variances.