You can do it in as little as two lines of code with blockproc() from the Image Processign Toolbox:
U = randi(9, 124, 7);
meanFilterFunction = @(theBlockStructure) mean(theBlockStructure.data(:));
blockSize = [4 1];
blockMeanArray = blockproc(U, blockSize, meanFilterFunction);
blockproc() moves along in "jumps" of the block size that you specify - it does not slide over by one element at a time like conv2() does. If you'd used conv2() you would have to subsample the result.