In the intro above the question is asked how one might go about computing $$\sum_{n=1}^\infty \frac{1}{3n(3n-1)}.$$
I can contribute an asymptotic expansion and an infinite series to this discussion, shown and proved below. In addition to posting this answer I am also asking two questions I. Can you verify this expansion using a different proof technique? II. Might there even be a reference to the below formula somewhere?
This is the proof, which is procedural, with the result at the end. Suppose we seek to evaluate
$$T(q) = \sum_{n\ge 1} \frac{1}{qn(qn-1)} =
\sum_{n=1}^p \frac{1}{qn(qn-1)} + \sum_{n\ge 1} \frac{1}{(qn+pq)(qn+pq-1)} \\=
\sum_{n=1}^p \frac{1}{qn(qn-1)} + \frac{1}{q^2}\sum_{n\ge 1} \frac{1}{(n+p)(n+(pq-1)/q)}$$
with $p, q\ge 2$ two integers.
The sum term, call it $S(x)$, is harmonic and may be evaluated by inverting its Mellin transform.
Recall the Mellin transform identity for harmonic sums with base function $g(x)$, which is
$$\mathfrak{M}\left(\sum_{k\ge 1}\lambda_k g(\mu_k x); s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s}\right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad
\mu_k = k \quad \text{and} \quad
g(x) = \frac{1}{(x+p)(x+(pq-1)/q)}.$$
The Mellin transform of $g(x)$ is
$$\int_0^\infty \frac{1}{(x+p)(x+(pq-1)/q)} x^{s-1} dx.$$
We chose $p\ge 2$ so that the base function does not have a pole at $x=-1$, which would cause the expansion about infinity that we eventually obtain not to converge. (As we evaluate the harmonic sum at $x=1$ this would be situated right on the boundary between the two disks about zero and about infinity, causing an evaluation at $x=1$ to fail.)
We evaluate the Mellin transform next. Use a circular contour to get
$$ g^*(s) (1 - e^{2\pi i (s-1)}) \\=
2\pi i \left(
\operatorname{Res}\left(\frac{x^{s-1}}{(x+p)(x+\frac{pq-1}{q})}; s=-p\right) +
\operatorname{Res}\left(\frac{x^{s-1}}{(x+p)(x+\frac{pq-1}{q})}; s=-\frac{pq-1}{q}\right)
\right) \\=
2\pi i \left(-q(-p)^{s-1} + q(-(pq-1)/q)^{s-1}\right)
= 2\pi i e^{\pi i (s-1)} q \left(((pq-1)/q)^{s-1}-p^{s-1}\right).$$
This implies that
$$g^*(s) = 2\pi i q \frac{-e^{\pi i s}}{1 - e^{2\pi i s}}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right)\\=
-2\pi i q \frac{1}{e^{-\pi i s} - e^{\pi i s}}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right) \\=
q\pi\frac{2i}{e^{\pi i s} - e^{-\pi i s}}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right) =
q\frac{\pi}{\sin(\pi s)} \left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right).$$
It follows that the Mellin transform of the sum term for $S(x)$ is given by
$$ Q(s) =
q \frac{\pi}{\sin(\pi s)}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right) \zeta(s)
\quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(s).$$
Finally we invert the Mellin transform with the inversion integral
$$\frac{1}{2\pi i}\int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$
to obtain an expansion about infinity of $S(x)$ starting with the pole at $s=2$ and getting
$$ S(x) \sim
- q \sum_{m\ge 2} (-1)^m
\left(\left(\frac{pq-1}{q}\right)^{m-1} - p^{m-1}\right) \frac{\zeta(m)}{x^m}.$$
Setting $x=1$ and substituting into the original equation we finally have
$$ T(q) = \sum_{n=1}^p \frac{1}{qn(qn-1)}
- \frac{1}{q} \sum_{m\ge 2} (-1)^m
\left(\left(\frac{pq-1}{q}\right)^{m-1} - p^{m-1}\right) \zeta(m).$$
Addendum Wed Apr 2 20:14:50 CEST 2014. The following MSE link shows how to work with divergent Mellin transforms where there is a pole on the positive real axis (even ones).