The Roman domination in a graph $G$ is a variant of the classical domination, defined by means of a so-called Roman domination function $f\colon V(G)\to \{0,1,2\}$ such that if $f(v)=0$ then, the vertex $v$ is adjacent to at least one vertex $w$ with $f(w)=2$ . The weight $f(G)$ of a Roman dominating function of $G$ is the sum of the weights of all vertices of $G$ , that is, $f(G)=\sum _{u\in V(G)}f(u)$ . The Roman domination number $\gamma _{R}(G)$ is the minimum weight of a Roman dominating function of $G$ . In this paper we propose algorithms to compute this parameter involving the $(\min,+)$ powers of large matrices with high computational requirements and the GPU (Graphics Processing Unit) allows us to accelerate such operations. Specific routines have been developed to efficiently compute the $(\min,+)$ product on GPU architecture, taking advantage of its computational power. These algorithms allow us to compute the Roman domination number of cylindrical graphs $P_{m}\Box ~C_{n}$ i.e., the Cartesian product of a path and a cycle, in cases $m=7,8$ ,9 $n\geq 3$ and $m\geq $ 10 $, n\equiv 0\pmod 5$ . Moreover, we provide a lower bound for the remaining cases $m\geq $ 10 $, n\not \equiv 0\pmod 5$ .