While going through the Kleiner and Lott notes "Notes on Perelman's papers", I encountered an argument that seems wrong to me, or (more likely) I do not understand something. It is about the $F$-funtional. It is defined as follows: $$F(g,f)=\int_{M}(R+\left | \bigtriangledown f \right |^2)e^{-f}dV,$$ where $g$ is Riemannian metric, $f$ is a smooth function over manifold $M$.
It is stated, that first variation of the functional can be expressed as follows: $$\delta F(v_{ij},h)=\int_{M}e^{-f}\left [ -v_{ij}(R_{ij}+\bigtriangledown_{i}\bigtriangledown _{j}f)+(\frac{v}{2}-h)(2\Delta f-\left | \bigtriangledown f \right |^2+R) \right ]dV,$$
where $v_{ij}=\delta g_{ij}$ is a symmetric covariant 2-tensor on $M$ (more specifically, it is an element of tangent space of infinite dimensional manifold of all smooth Riemannian metrics on $M$), $h=\delta f$ is a smooth function over $M$ (f is a smooth function as well), $v=g^{ij}v_{ij}$.
Part of the proof is as follows: $$\delta R=-\Delta v+\triangledown_{i}\triangledown_{j}v_{ij}-R_{ij}v_{ij}$$ $$\delta (e^{-f}dV)=(\frac{v}{2}-h)e^{-f}dV$$ $$\delta \left | \triangledown f \right |^2=-v^{ij}\triangledown_{i}f\triangledown_{j}f+2<\triangledown f,\triangledown h>$$ I understand these three lines very clearly - it is not very difficult. Putting all together I get $$\delta F=\int_{M}\left [-\Delta v +\triangledown_{i}\triangledown_{j}v_{ij}-R_{ij}v_{ij} +(R+\left | \triangledown f \right |^2)(\frac{v}{2}-h) +2<\triangledown f,\triangledown h>-v^{ij}\triangledown_{i}f\triangledown_{j}f \right ]e^{-f}dV$$.
And that is were my results do not match Kleiners and Lotts result. The difference is that in their derivation they have $v_{ij}\triangledown_{i}f\triangledown_{j}f$ instead of $v^{ij}\triangledown_{i}f\triangledown_{j}f$. But this part is just direct plug-in of those three variation formulas above. And those are exactly the same in my derivation and in authors.
So, what I am missing? What I'm not seeing?
I just realized, that I completely forgot one aspect. When I was thinking about repeated indices as a summation over them, I forgot, that the summation had to be done with respect to the metric $g$. For example $$div\alpha_{ij}=g^{jk}\triangledown_{j}\alpha_{ki}=\triangledown_{j}\alpha_{ji},$$
where $\alpha$ is a covariant 2 tensors, and the last equation reflected the convention not to bother to raise indices and to sum over repeated indices. Applying this in my case, we have $$v^{ij}\triangledown_{i}f\triangledown_{j}f=g^{ip}g^{jq}v_{pq}\triangledown_{i}f\triangledown_{j}f=v_{ij}\triangledown_{i}f\triangledown_{j}f$$
First expression is just a formula that I got, then lowering the indices we get the second expression and the convention that I described above leads to the third expression which is what Kleiner and Lott had.
Thus, my mistake was that I forgot to always have in mind the summation with respect to the metric.
Thanks for @mike and @Anthony Carapetis for forcing me to think more about the proper use of indices.