As an exercise in learning about "divergent series" and their classical summation-procedures I played with a method which is based on a matrix-transformation using the "Eulerian numbers". It is often "stronger" than Euler-summation (of any order), can sum not only the geometric series with $q<1$ but also the alternating series of factorials (it is similar to Borel-summation) - but cannot sum the nonalternating versions of that series when they are classically divergent, and thus cannot sum not the harmonic series.
The pure series transformation using the matrix $E$ from a series $a_1+a_2+a_3+...$ into a vector $T$ containing the transforms looks like
$$ [a_1,a_2,a_3,...] \cdot E = [t_1, t_2,t_3,...] = T $$
and for convergent and summable cases the partial sums of the $t_k$ approximate the sum of the series-terms on the lhs and we have then a usable series-summation method (it gives also convergence acceleration for the summable cases). For reference, I call this provisorically the "Eulerian transform" or "Eulerian summation".
For the given case of the harmonic series the method gives for the partial sums of the $t_k$ infinity and thus the harmonic series is not summable by the "Eulerian summation".
But besides this that partial sums give the following remarkable approximation (writing the $N$'th harmonic number as $H_N$)
$$ \sum_{k=1}^N t_k - H_N = \log 2 +\varepsilon_N
\qquad \qquad \lim_{N \to \infty} \varepsilon_N=0$$
which reminds of the formula for the Euler-Mascheroni-constant, but has a different constant term.
Remark: the method has not yet a sufficient theoretical framework but because the entries in $E$ have a tractable composition one can do relatively simple partial analyses for the self-study.
See
initial exploration of E and the basic idea of summation method and for the current problem the
investigation of the transform for some non-summable cases