First, we repeat the arguments from this [stackexchange answer][1].  $P^{-1}$ is an $M$-matrix, and can thus be written as $s(I-A)$ for some positive $s$ and some $A$ with non-negative entries.  As $P^{-1}$ is positive definite, the spectrum of $A$ lies to the left of $\{ z: \hbox{Re}(z) = 1 \}$, and hence by Perron-Frobenius the spectral radius of $A$ is less than $1$.  Thus we have the absolutely convergent Neumann series
$$ P = s^{-1} (1 + A + A^2 + \dots )$$
and hence
$$ PVPVP = s^{-3} \sum_{i=0}^\infty \sum_{j=0}^\infty \sum_{k=0}^\infty A^i V A^j V A^k.$$
It thus suffices to show that
$$ \sum_{i+j+k=m} A^i V A^j V A^k \quad (1)$$
has non-negative coefficients for each $m \geq 0$ (where $i,j,k$ are understood to be non-negative integers).    By change of variables, this is
$$ \sum_{0 \leq q \leq r \leq m} A^q V A^{r-q} V A^{m-r}.$$
Writing $A = (a_{st})_{1 \leq s,t \leq n}$ and $V = \hbox{diag}(v_1,\dots,v_n)$, the $st$ coefficient of (1) can be expanded as
$$ 
\sum_{s=s_0,s_1,\dots,s_m=t} a_{s_0 s_1} \dots a_{s_{m-1} s_m}\sum_{0 \leq q \leq r \leq m} v_{s_q} v_{s_r}.$$
But the quadratic form
$$ \sum_{0 \leq q \leq r \leq m} x_q x_r = \frac{1}{2}(x_0+\dots+x_m)^2 + \frac{1}{2} x_0^2 + \dots + \frac{1}{2} x_m^2$$
is positive definite, and the $a_{st}$ are non-negative, and the claim follows.


  [1]: http://math.stackexchange.com/a/989748