Solving Dynamic Programming Problem by
Pipeline Implementation on GPU

In this study, we implement a typical DP program on GPU. To simplify the exposition, we focus on programs that solve such a simplified DP problem defined as follows:

For example, Fibonacci number problem is equal to the S-DP problem where $k=2,a_{1}=2,a_{2}=1$, $\otimes=+$, and ST$[0]$=ST$[1]$=$1$.

To begin with, we show a straightforward sequential algorithm that solves the S-DP problem. Fig. 1 shows the algorithm.

The outer loop computes values from ST$[a_{1}]$ to ST$[n-1]$ in order, and the inner loop computes ST[i] for each i by equation (1). Since the outer loop takes $n-a_{1}+1=O(n)$ steps and the inner loop requires $O(k-1)$ steps, this algorithm takes $O(nk)$ steps in total.

Next, we consider parallelizing the algorithm for S-DP problem. The straightforward approach is to parallelize the inner loop by using GPU cores. We can easily write a multi-thread program that executes the inner loop-body, $\mbox{{\tt ST}}{\tt[i]}=\mbox{{\tt ST}}{\tt[i]}\;\otimes\;\mbox{{\tt ST}}{\tt[i}-a_{\tt j}{\tt]}$, for each j in parallel using $k-1$ threads at one time. Such an implementation, however, does not improve the time cost, because every thread has access to the same ST[i] and thus memory access conflicts occur. As a result, those memory conflicts should be automatically solved at run-time by the serializing mechanism of GPU, and consequently the whole time cost stays in $O(nk)$ steps, which is the same time cost as that of the sequential implementation.

To avoid the memory access conflicts, we can use a well-known standard parallel prefix computation algorithm [6, 7], in which the computations of $\;\otimes\;$ over $k$ values are executed in a tournament fashion. Since the parallel prefix computation runs in $O(\log k)$ steps for $k$ values, obviously the entire time cost can be improved to $O(n\log k)$ steps by using $k$ threads.

Although we can successfully reduce the time cost from $O(nk)$ to $O(n\log k)$ by using the parallel prefix computation, it is not work-time optimal because there are many idle threads during the computations in a tournament fashion. In the next section we propose other parallel implementation strategy and show that we can improve the time cost further.

In our implementation, we use a group of $k$ threads to establish $k$-stage pipeline, and this thread group treats $k$ consecutive elements at one time in parallel. Fig. 2 describes our pipeline algorithm for the S-DP problem. The index variable i of the outer loop stands for the head position of the $k$-thread group. The inner loop controls each thread’s behaviour in such a way that the j-th thread executes computation for ST[i-j+1] using the value stored in ST[i-j+1-$a_{\tt j}$].

An execution example is shown in Fig. 3, where $k=3$, $a_{1}=5$, $a_{2}=3$, and $a_{3}=1$ hold and the initial values are already stored in ST[0], ST[1],$\ldots$ , ST[4].

In Step 1, the head position i of the thread group is 5. In this step the only one thread is activated and executes ST[5] $\leftarrow$ ST[0]. In Step 2, the head position is incremented to 6, and two threads are activated. The first thread treats ST[6] and the second thread works on ST[5]. In Step 3, the head position becomes 7, and now all $k=3$ threads actively execute operations for ST[7], ST[6], and ST[5] respectively. It should be noted that finally in Step 3 the content of ST[5] is completely determined while those of ST[7] and ST[6] are partially computed and not yet determined. From Step 3, all the $k=3$ threads are active until Step $n-a_{1}$ when the head position i of thread group reaches $n-1$, and after that step the number of active threads decreases one by each step. As you can see there is no memory access conflict in this example.

As for the time-complexity of our pipeline implementation, from a theoretical viewpoint, it takes only $O(n)$ steps, because the outer loop takes $n+k-a_{1}-1=O(n)$ cycles and the inner loop requires $O(1)$ time if there is no adjacent offset pair $(a_{m},a_{m+1})$ such that $a_{m}=a_{m+1}+1$.

However, from a practical viewpoint, because of the memory access conflicts, the inner loop may take more time steps. Actually, in the worst case when consecutive offset numbers are given, those ST$[i_{\tt j}-a_{\tt j}]$, in the right hand side of the assignment statement, coincidentally become the same element of array ST  and hence the worst memory access conflicts occur. In such a case, all threads in the inner loop are serialized and it takes time proportional to $k$. See Fig. 4 for such a worst case example. In this example, all four threads try to have access to ST[${\tt i}-4$] at the same time in the inner loop.

Let $seq=(a_{p},a_{p+1},\ldots,a_{q})$ be one of the longest sub-sequences of given offset numbers $(a_{1},a_{2},\ldots,a_{k})$ satisfying $a_{r}=a_{r+1}+1$ for all $p\leq r<q$. Then, it is easy to check that in the inner loop every thread $r$ ($p\leq r<q$) has access to the same element of array ST, and as a result those conflicts are serialized at run time by GPU’s serializing mechanism and hence the memory access time becomes $(q-p+1)$ times slower than that of conflict-free case. For such a case, in [5], we proposed a 2-by-2 pipeline implementation technique where each thread invoked in the inner loop executes two computations for each element of array ST. The details can be found in [5].

Before going to the next section, we show the performance of our pipeline implementation on GPU. We use a computer with 3.40 [GHz] Intel Xeon CPU E3-1245 v3 and an NVIDIA GeForce GTX TITAN Black. The OS is Ubuntu 17.10 and we use g++ 7.2.0 for compiling cpp programs and CUDA 9.2 [8]. The experimental results are shown in TABLE I.

TABLE I shows the average execution time for each implementation on GPU. In the table, SEQUENTIAL, NAIVE-PARALLEL, and PIPELINE respectively stand for the sequential implementation, the naive multi-thread implementation, and our pipeline implementation (for a general case). Here, we use $\min$ operation for $\otimes$. The average is computed among 100 executions for each setting.

As for the comparison between the sequential implementation and the parallel ones, parallel implementations are much faster even though it is NAIVE-PARALLEL. Although there is no difference in time between NAIVE-PARALLEL and PIPELINE until $n\leq 2^{17}$, PIPELINE is faster than NAIVE-PARALLEL when $n\geq 2^{18}$.

It is well-known that the MCM problem can be efficiently solved by DP with a two-dimensional solution table of triangular shape, and that each element is computed along the diagonal direction. See Fig. 5 for an example. In the figure, each number represents the order of elements to be computed.

The content of each element is computed by a binary operator $\downarrow$ returning a smaller operand and a binary function $f(*,*)$ as in Fig. 6. The detailed definition of the MCM problem can be found in [1]. In the figure, the element marked 13 is computed by the elements marked by 1, 6, 10, 11, 8, and 4. If we write ST[$x$] for the element marked $x$ here, the computation is expressed as

Since the elements of two-dimensional solution table are computed in a total order (linear order), we can line up them into a linear array according to that order. Once the solution table is transformed into a linear array, we can apply the pipeline technique to the MCM problem as well. An execution example is shown in Fig. 7.

Here, we explain how to solve the MCM problem by pipeline implementation in details.

Firstly, we assume that the original two-dimensional solution table of diagonal shape is already mapped to a one dimensional solution table ST appropriately. That is, we line up all the elements of the two-dimensional table into a linear array ST according to the total (linear) order in which each element is computed along the diagonal direction as in Fig. 5 and 7. Since the number of input matrices is $n$, the number of elements in the solution table is $n(n+1)/2$. First $n$ elements of array ST are preset with initial values, those correspond to the elements located in the diagonal line of the original two-dimensional table.

Let each computation for ST[$i$] be represented as

It should be noted that here each $k_{i}$ may differ from others and that it is $n-1$ at largest. See Fig. 5 and 7 for an example. Here, we have $n=5$, and ST[13] and ST[12] can be represented as follows:

and

To solve such defined MCM problem by the pipeline implementation on GPU, we have to modify the pipeline algorithm designed for the S-DP problem, because the offset numbers for each ST[$i$] may differ from others. Thus, for the MCM problem, we propose a modified pipeline algorithm in Fig. 8. In the MCM-pipeline algorithm, each element of ST is computed by equation (2), and ST$[1]$, ST[2],$\ldots$ ,ST$[n]$ are preset with initial values. In substep 1, 2, and 3, the computation of $f(l_{(i,j)},r_{(i,j)})$ is executed, and in substep 4, that obtained value is used for the computation by $\downarrow$ and the result is stored to ST$[i]$.

In this subsection, we prove that no memory access conflict occurs during the execution of the MCM-pipeline algorithm described in Fig. 8.

To begin with, we prove the following lemma.

Proof: Assume that threads $p$ and $q$ try to read the same element of ST in substep 1 and that $p<q$ holds. Let $(row_{p},col_{p})$ (resp. $(row_{q},col_{q})$) be the pair of row and column indexes of the elements in the original two-dimensional solution table of triangle shape of the MCM problem for which thread $p$ (resp. $q$) is now computing. Since threads $p$ and $q$ try to read the same element of ST and in substep 1 they read the value for the left argument of the function $f$, the relation $row_{p}=row_{q}$ must hold. Then, since threads $p$ and $q$ respectively read the $p$-th and $q$-th elements from the left in the same row of the original two-dimensional solution table, the relation $p=q$ must hold if the two threads read the same element of ST, which contradicts the assumption $p<q$.  

Next, we prove the following lemma, which can be proved in a similar way to the proof of Lemma 1.

Proof: Assume that threads $p$ and $q$ try to read the same element of ST in substep 2 and that $p<q$ holds. Let $(row_{p},col_{p})$ (resp. $(row_{q},col_{q})$) be the pair of row and column indexes of the elements in the original two-dimensional solution table of triangle shape of the MCM problem for which thread $p$ (resp. $q$) is now computing. Since threads $p$ and $q$ try to read the same element of ST and in substep 2 they read the value for the right argument of the function $f$, the relation $col_{p}=col_{q}$ must hold. On the other hand, thread $p$ reads the $p$-th elements below of the row $row_{p}$ of the original two-dimensional solution table, and thread $q$ does the $q$-th elements below of the row $row_{q}$ of the table. Since threads $p$ and $q$ read row $(row_{p}+p)$ and row $(row_{q}+q)$ respectively, the relation $row_{p}+p=row_{q}+q$ must hold if the two threads read the same element of ST. Since $p<q$ and $col_{p}=col_{q}$ hold, the relation $row_{p}<row_{q}$ must hold from the way of mapping from the original two-dimensional solution table to the linear array ST. This leads to the relation $row_{p}+p\not=row_{q}+q$, which contradicts the assumption that threads $p$ and $q$ read the same element of array ST.  

In substep 3, each thread simply executes computation using local variables. In substep 4, it is obvious that each thread has access to a distinct element of ST. Therefore, by Lemma 1 and Lemma 2, we obtain the following theorem.

As for the time-complexity of the MCM-pipeline implementation, from a theoretical viewpoint, it takes only $O(n^{2})$ steps with $(n-1)$ threads, because the outer loop takes $O(n^{2})$ cycles and the inner loop requires $O(1)$ time (Theorem 1).