By using the generating function of Stirling numbers (https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind) we can fairly easily find the following identity: \begin{equation} \sum\limits_{\begin{array}{r} n_1+n_2=n \\n_1\ge1,n_2\ge1\end{array}} \frac{2 H_{n_1-1}}{n_1} \cdot \frac{1}{n_2} = 3 \frac{[H_{n-1}]^2 - H^{(2)}_{n-1}}{n} \end{equation} This follows from computing the power series expansion of $[\log(1-x)]^3$ in two different ways firstly from the generating function and secondly using convolutions.
Now, the question is how do we compute more generic convolutions like the following: \begin{equation} {\mathfrak S}_q^{(p)}(n):=\sum\limits_{\begin{array}{r} n_1+n_2=n \\n_1\ge1,n_2\ge1\end{array}} \frac{H^{(p)}_{n_1-1}}{n_1^q} \frac{1}{n_2} = ? \end{equation} where $p,q$ are positive integers. Note that if we didn't have the harmonic numbers in the numerator we can always compute the convolution using partial fraction decomposition and the definition of generalized harmonic numbers. However in our case this technique does not seem to work. How do we proceed in here?
We demonstrate that the convolution above can be computed by elementary methods. We use the following identities: \begin{eqnarray} H_{n-1}^{(p)} &=& \sum\limits_{k=1}^\infty \left( \frac{1}{k^p} - \frac{1}{(n-1+k)^p} \right) \\ \frac{1}{n_1^q (n-n_1)} &=& \left(\sum\limits_{l=1}^q \frac{1}{n_1^{q+1-l}} \frac{1}{n^l}\right) + \frac{1}{(n-n_1)} \frac{1}{n^q} \end{eqnarray} and \begin{eqnarray} &&\frac{1}{n_1^q (n-n_1) (n_1-1+k)^p}=\\ && \sum\limits_{l_1=1}^q \frac{(-1)^{q-l_1}}{(k-1)^{q+p-l_1}} \binom{q+p-1-l_1}{p-1} \frac{1}{(n-n_1)} \frac{1}{n^{l_1}}+ \sum\limits_{l_1=1}^p \frac{(-1)^{q}}{(k-1)^{q+p-l_1}} \binom{q+p-1-l_1}{q-1} \frac{1}{(n-n_1)} \frac{1}{(n+k-1)^{l_1}}+ \sum\limits_{1 \le l_2 \le l_1 \le q} \frac{(-1)^{q-l_1}}{(k-1)^{q+p-l_1}} \binom{q+p-1-l_1}{p-1} \frac{1}{n_1^{l_2}} \frac{1}{n^{l_1+1-l_2}} + \sum\limits_{1 \le l_2 \le l_1 \le p} \frac{(-1)^{q}}{(k-1)^{q+p-l_1}} \binom{q+p-1-l_1}{q-1} \frac{1}{(n_1+k-1)^{l_2}} \frac{1}{(n+k-1)^{l_1+1-l_2}} \end{eqnarray} The last two equations are just simple fraction decompositions of the quantity on the left hand side in the variable $n_1$. Inserting the identities above into the definition of our convolution allows us to do the sum over $n_1$ straight away. We have: \begin{eqnarray} &&{\mathfrak S}_q^{(p)}(n) = \\ &&\sum\limits_{l=1}^q H_{n-1}^{(q+1-l)} \frac{\zeta(p)}{n^l} - \sum\limits_{l=1}^{q+p} H_{n-1}^{(q+p+1-l)} \frac{1}{n^l}+ H_{n-1} \left( \frac{\zeta(p)}{n^q} - \frac{1}{n^{q+p}}\right) -\sum\limits_{2 \le l_2 \le l_1 \le q} (-1)^{q-l_1} \zeta(q+p-l_1) \binom{q+p-1-l_1}{p-1} \frac{H_{n-1}^{(l_2)}}{n^{l_1+1-l_2}} -{\mathfrak R}_q^{(p)}(n) - {\bar {\mathfrak R}}_q^{(p)}(n) \end{eqnarray} Here the remaining terms read: \begin{eqnarray} {\mathfrak R}_q^{(p)}(n) &:=& \sum\limits_{l_1=1}^p \sum\limits_{k=1}^\infty \frac{(-1)^q}{(k+n)^{l_1} k^{q+p-l_1}} \binom{q+p-1-l_1}{q-1} H_{n-1} \\ &=& \sum\limits_{l_1=1}^p \left( \binom{q+p-1-l_1}{q-1} H_{n-1} (-1)^{l_1+p}\right.\\ &&\left. \sum\limits_{l_2=1}^{l_1-1} \left(\zeta (l_2+1)-H_n^{(l_2+1)}\right) n^{l_2-p-q+1} \binom{-l_2+p+q-2}{-l_1+p+q-1} + \sum\limits_{l_2=1}^{q+p-l_1-1} \zeta (l_2+1) (-1)^{-l_2-1} n^{l_2-p-q+1} \binom{-l_2+p+q-2}{l_1-1} \right) \\ {\bar {\mathfrak R}}_q^{(p)}(n) &:=& \sum\limits_{1 \le l_2 \le l_1 \le p} \sum\limits_{k=1}^\infty \frac{(-1)^q}{(k+n)^{l_1+1-l_2} k^{q+p-l_1}} \binom{q+p-1-l_1}{q-1} \left(H_{n-1+k}^{(l_2)} - H_k^{(l_2)}\right) \end{eqnarray} The second quantity is more complicated but we believe it can be handled in a similar way.