Double harmonic sum $\sum_{n\geq 1}\frac{H^{(p)}_nH_n}{n^q}$ Are there any general formula for the following series 
$$\tag{1}\sum_{n\geq 1}\frac{H^{(p)}_nH_n}{n^q}$$
Where we define 
$$H^{(p)}_n= \sum_{k=1}^n \frac{1}{k^p}\,\,\,\,\,H^{(1)}_n\equiv H_n =\sum_{k=1}^n\frac{1}{k} $$
For the special case $p=q=2$ in (1) I found the following paper
Stating that 
$$\sum_{n\geq 1}\frac{H^{(2)}_nH_n}{n^2}=\zeta(5)+\zeta(2)\zeta(3)$$
See equation (3a) .
Is there any other paper in the literature discussing (1) or any special cases ?

In here we provide a generating function of the quantities in question. Let us define:
\begin{equation}
{\bf H}^{(p,r)}_q(t) := \sum\limits_{m=1}^\infty H_m^{(p)} H_m^{(r)} \frac{t^m}{m^q}
\end{equation}
In here we take $q\ge1$. We have:
\begin{eqnarray}
&&{\bf H}^{(p,1)}_q(t) = Li_p(1) \cdot \frac{1}{2} [\log(1-t)]^2 \cdot 1_{q=1}+\\
&&\frac{(-1)^{q}}{2} \sum\limits_{l=(q-2)}^{p+q-3}  \left(\binom{l}{q-2} 1_{l < p+q-3} + ({\mathcal A}^{(p)}_{q-2}) 1_{l=p+q-3}\right) \cdot \underbrace{\int\limits_0^1 \frac{[Li_1(t \xi)]^2}{\xi} Li_{l+1}(\xi)  \frac{[\log(1/\xi)]^{p+q-3-l}}{(p+q-l-3)!}d\xi}_{I_1}+\\
&& \frac{(-1)^{q-1}}{2} \sum\limits_{j=0}^{q-3} \left({\mathcal A}^{(p)}_{q-2-j}\right) \cdot \zeta(p+q-2-j) \underbrace{\int\limits_0^1 \frac{[Li_1(t \xi)]^2}{\xi} \frac{[\log(\xi)]^j}{j!} d\xi}_{I_2}+\\
&& \sum\limits_{l=1}^p  \underbrace{\int\limits_0^1 \frac{Li_q(t \xi)}{\xi} Li_l(\xi) \frac{[\log(1/\xi)]^{p-l}}{(p-l)!} d\xi}_{I_3}
\end{eqnarray}
Here $t\in (-1,1)$ and $p=1,2,\cdots$ and
\begin{equation}
{\mathcal A}^{(p)}_{q} := p+\sum\limits_{j=2}^{q} \binom{p+j-2}{j}
= p \cdot 1_{p=1} + \frac{p+q-1}{p-1} \binom{p+q-2}{q}\cdot 1_{p > 1}
\end{equation}
Note 1: The quantities in the right hand side all contain products of poly-logarithms and a power of logarithm. Those quantities, in principle, have been already dealt with in An integral involving product of poly-logarithms and a power of a logarithm. for example.
Note 2: Now that we have the generating functions we will find recurrence relations for the sums in question  and hopefully provide some closed form expressions .
Now we have:
\begin{eqnarray}
&&I_1 =\\
&&\sum\limits_{l_1=2}^{l+1} \binom{p+q-2-l_1}{l+1-l_1} (-1)^{l+1-l_1} \zeta(l_1) \left({\bf H}^{(1)}_{p+q-l_1}(t) - Li_{p+q+1-l_1}(t) \right)+\\
&&\sum\limits_{l_1=2}^{p+q-2-l} \binom{p+q-2-l_1}{l} (-1)^{l-1} \zeta(l_1) \left({\bf H}^{(1)}_{p+q-l_1}(t) - Li_{p+q+1-l_1}(t) \right)+\\
&&\sum\limits_{l_1=1}^{p+q-2-l} \binom{p+q-2-l_1}{l} (-1)^{l-0} \left( {\bf H}^{(l_1,1)}_{p+q-l_1}(t) - {\bf H}^{(l_1)}_{p+q+1-l_1}(t) \right)
\end{eqnarray}
and
\begin{eqnarray}
&&I_2=2 (-1)^j \left[ {\bf H}^{(1)}_{j+2}(t) - Li_{j+3}(t)\right]
\end{eqnarray}
and
\begin{eqnarray}
&&I_3=\\
&&\sum\limits_{l_1=2}^l \binom{p-l_1}{p-l}(-1)^{l-l_1} \zeta(l_1) Li_{p+q+1-l_1}(t) +\\
&& \sum\limits_{l_1=2}^{p-l+1} \binom{p-l_1}{l-1}(-1)^{l} \zeta(l_1) Li_{p+q+1-l_1}(t)-
\sum\limits_{l_1=1}^{p-l+1} \binom{p-l_1}{l-1} (-1)^l {\bf H}^{(l_1)}_{p+q+1-l_1}(t)
\end{eqnarray}
