Double harmonic sum $\sum_{n\geq 1}\frac{H^{(p)}_nH_n}{n^q}$

A couple of places to start for further looking:

The paper Further summation formulae related to generalized harmonic numbers by Zheng mentions the $p = q = 2$ case in example 2.3, and has several related results, but not explicitly of the form you're looking for.

Appendix B of the paper On some log-cosine integrals related to $\zeta(3), \zeta(4)$ and $\zeta(6)$ by Mark Coffey has some similar things, and maybe looking through its references and papers that cite it would yield more.

Here we provide an answer for $p=1$ and arbitrary $q \ge 2$. Let us denote the quantity in question as follows: \begin{equation} {\mathfrak S}^{(p,1)}_q := \sum\limits_{m=1}^\infty \frac{H^{(p)}_m H_m}{m^q} \end{equation} In order to calculate the above we start from the generating function. We have: \begin{eqnarray} \sum\limits_{m=1}^\infty H^{(p)}_m H_m \frac{t^m}{m} &= & - \int\limits_0^t \frac{\log(1-\xi/t)}{\xi} \cdot \frac{Li_p(\xi)}{1-\xi} d\xi \\ &=& \left\{ \begin{array}{rr} -\frac{1}{3} [\log(1-t)]^3 -\log(1-t) Li_2(t) + Li_3(t) & \mbox{for $p=1$} \\ \cdots \end{array} \right. \end{eqnarray} Now we use the common trick,i.e. we divide by $t$ and multiply by an appropriate power of $\log(1/t)$ and integrate appropriately. We have: \begin{eqnarray} {\mathfrak S}^{(1,1)}_q &=& \underbrace{-\frac{1}{3} \int\limits_0^1 \frac{[\log(1/\xi)]^{q-2}}{(q-2)!} \cdot \frac{[\log(1-\xi)]^3}{\xi} d\xi}_{I_1} +\\ &&\underbrace{\int\limits_0^1 \frac{[\log(1/\xi)]^{q-2}}{(q-2)!} \cdot \frac{Li_1(\xi) Li_2(\xi)}{\xi} d\xi}_{I_2}+\\ &&\underbrace{\int\limits_0^1 \frac{[\log(1/\xi)]^{q-2}}{(q-2)!} \cdot \frac{ Li_3(\xi)}{\xi} d\xi}_{I_3} \end{eqnarray} Clearly all the inthree integrals in here have already been dealt with in MSE. The last integral is trivial . We have: \begin{equation} I_3= Li_{q+2}(+1) \end{equation} The second integral has been evaluated in An integral involving product of poly-logarithms and a power of a logarithm.. We have: \begin{equation} I_2= \frac{1}{2}\left[ \zeta(2) \zeta(q) + \sum\limits_{j=1}^{q-3} j \zeta(q-j-1) \zeta(j+3)-\sum\limits_{j=1}^{q-2} j {\bf H}^{(q-j-1)}_{j+3}(+1) \right] \end{equation} where the quantities highlighted in bold font have been calculated in Calculating alternating Euler sums of odd powers .They all reduce to single zeta values for $q+2\le 7$ and otherwise involve a few additional two-dimensional zeta values. Finally the first integral has been evaluated in Compute an integral containing a product of powers of logarithms. .It always reduces to single zeta values. We have: \begin{eqnarray} &&I_1 = \\ &&\frac{(-1)^{q-1}}{(q-1)!} \left[ -\frac{1}{3} \Psi^{(q+1)}(1) + \frac{1}{2} \sum\limits_{j=1}^{q-2} \binom{q-1}{j} \left\{\Psi^{(j+1)}(1)\Psi^{(q-1-j)}(1) + \Psi^{(j+0)}(1)\Psi^{(q-j)}(1)\right\}+\right.\\ && \left. -\frac{1}{3} \sum\limits_{1 \le j < j_1 \le q-2} \frac{(q-1)!}{j! (j_1-j)!(q-1-j_1)!} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-1-j_1)}(1) \right] \end{eqnarray} where $\Psi^{(j)}(1)= (-1)^{j+1} j! \zeta(j+1)$ for $j=1,2,3,\cdots$.