I'm interested in finding the quickest way of grouping the elements of a large matrix in sub-groups of NxM elements and them summing them together.
To be completely clear, I'm actually not interested in the "regrouped" matrice, but only in the final results where the elements are summed.

I'll show you an example below:

Say I have the following matrix 9x8:

test = Array[Subscript[a, ##] &, {8, 9}]

I regroup it in sub-matrices NxM, in this example 3x2:

subtest = Partition[test, {2, 3}]

and then I sum them together (as suggested in the comment by @:

out = MapAt[Total[#, -1] &, subtest, {All, All}];

I could use other ways of summing the subgroups, as for example:

    out = Total /@ Flatten /@ # & /@ subtest;

Or using two nested tables, or for loops, etc.

My question is what is the fastest method for doing this? I need to do it on a 48k x 48k matrix, so I'd really need something reasonably quick.
Should I look into compiling nested for loops in C (not sure, I haven't ever tried)?

Something worth mentioning is that the entries of the matrix are all integers larger or equal to 0.

I'll add a (redundant) example with numeric values, thac can be however modified to larger matrices:

test = RandomInteger[1, {8, 9}];

{{0, 0, 1, 0, 1, 1, 0, 0, 0}, {1, 1, 0, 1, 0, 0, 1, 0, 0}, {0, 0, 1,
1, 1, 1, 0, 0, 0}, {0, 0, 1, 0, 0, 1, 1, 0, 1}, {0, 0, 1, 1, 0, 0,
1, 1, 1}, {1, 0, 0, 1, 0, 1, 1, 1, 1}, {1, 0, 1, 1, 1, 0, 1, 1,
0}, {0, 0, 1, 1, 1, 0, 1, 0, 1}}

m = 3

n = 2

out = MapAt[Total[#, -1] &, Partition[test, {n, m}], {All, All}]

{{3, 3, 1}, {2, 4, 2}, {2, 3, 6}, {3, 4, 4}}

Partition[test, {2,3}] is quite slow in this case because it has to rearrange the elements in the data vector that represents the entries of a packed array in the backend:

Flatten[test] == Flatten[Partition[test, {2, 3}]]

False

Using Span (;;) as follows employs 6 monotonically increasing read operations; in this specific case, these operations are faster than using Partition:

n = 24000;

test = RandomInteger[1, {n, n}];



a = Total[Partition[test, {2, 3}], {3, 4}]; // 

  AbsoluteTiming // First

b = Sum[test[[i ;; ;; 2, j ;; ;; 3]], {i, 1, 2}, {j, 1, 3}]; // 

  AbsoluteTiming // First

a == b

245.89

117.943

True

However, this performance advantage seems to decay when the matrix test becomes bigger (so swapping is required). E.g., for $n = 4800$, method b is aboutten times faster thana`, but for $n = 24000$, it's only a factor of 4.6 and here it has degraded to a factor of 2 or so...

SparseArray method

Have I said already that I love SparseArrays?

AbsoluteTiming[

  c = Dot[

    KroneckerProduct[

     IdentityMatrix[n/2, SparseArray], 

     ConstantArray[1, {1, 2}]

     ],

    Dot[

     test,

     KroneckerProduct[

      IdentityMatrix[n/3, SparseArray], 

      ConstantArray[1, {3, 1}]

      ]

     ]

    ]

  ][[1]]

a == c

76.3822

True

The story goes on...

A combination of the SparseArray method from above with a CompiledFunction):

cf = Compile[{{x, _Integer, 1}, {k, _Integer}},

   Table[

    Sum[Compile`GetElement[x, i + j], {j, 1, k}],

    {i, 0, Length[x] - 1, k}],

   CompilationTarget -> "C",

   RuntimeAttributes -> {Listable},

   Parallelization -> True,

   RuntimeOptions -> "Speed"

   ];

d = KroneckerProduct[

     IdentityMatrix[n/2, SparseArray], 

     ConstantArray[1, {1, 2}]

     ].cf[test, 3]; // AbsoluteTiming // First

a == d

33.5677

True