  We extend several celebrated methods in classical analysis for summing series
of complex numbers to series of complex matrices. These include the summation
methods of Abel, Borel, Ces\'aro, Euler, Lambert, N\"orlund, and
Mittag-Leffler, which are frequently used to sum scalar series that are
divergent in the conventional sense. One feature of our matrix extensions is
that they are fully noncommutative generalizations of their scalar counterparts
-- not only is the scalar series replaced by a matrix series, positive weights
are replaced by positive definite matrix weights, order on $\mathbb{R}$
replaced by Loewner order, exponential function replaced by matrix exponential
function, etc. We will establish the regularity of our matrix summation
methods, i.e., when applied to a matrix series convergent in the conventional
sense, we obtain the same value for the sum. Our second goal is to provide
numerical algorithms that work in conjunction with these summation methods. We
discuss how the block and mixed-block summation algorithms, the Kahan
compensated summation algorithm, may be applied to matrix sums with similar
roundoff error bounds. These summation methods and algorithms apply not only to
power or Taylor series of matrices but to any general matrix series including
matrix Fourier and Dirichlet series. We will demonstrate the utility of these
summation methods: establishing a Fej\'{e}r's theorem and alleviating the Gibbs
phenomenon for matrix Fourier series; extending the domains of matrix functions
and accurately evaluating them; enhancing the matrix Pad\'e approximation and
Schur--Parlett algorithms; and more.
