my $A = [ [ 1, 2 ],
              [ 3, 4 ] ];   # matrix A

    my $B = [ [ 5, 6 ],
              [ 7, 8 ] ];   # matrix B

    my $mats = [ $A, $B ];  # array of matrices

##</code><code>##

    # [ [ f(1,5), f(2,6) ],
    #   [ f(3,7), f(4,8) ] ]

##</code><code>##

    #                            c=0  c=1
    #                           ________
    #   Matrix A (m=0)         /   /   /
    #                    r=0  / 1 / 2 /.
    #                        /___/___/ .
    #                       /   /   /  .
    #                 r=1  / 3 / 4 /   .
    #                     /___/___/    .
    #                                  .
    #                           ________
    #   Matrix B (m=1)         /   /   /
    #                    r=0  / 5 / 6 /
    #                        /___/___/
    #                       /   /   /
    #                 r=1  / 7 / 8 /
    #                     /___/___/
    #
    #                   c=0  c=1

##</code><code>##

    my $result;
    for my $r (0 .. @{$mats->[0]} - 1) {
        for my $c (0 .. @{$mats->[0][0]} - 1) {
            $result->[$r][$c] = f( map $_->[$r][$c], @$mats );
        }
    }

##</code><code>##

    use Data::Dumper;

    sub f {
        "f(" . join(", ", map dumpf($_), @_) . ")";
    }

    sub dumpf {
        local $Data::Dumper::Indent = 0;
        local $Data::Dumper::Terse = 1;
        local $_ = Dumper(@_);
        tr/\'//d;
        $_;
    }

##</code><code>##

    print dumpf($result), $/;
    # [ [f(1,5), f(2,6)],
    #   [f(3,7), f(4,8)] ]

##</code><code>##

    sub mapat {
        my $depth = shift;
        my $f     = shift;
        $depth
            ? map [ map mapat($depth-1, $f, $_), @$_ ], @_
            : $f->(@_);
    }

##</code><code>##

    for (0..3) {
        print "mapat $_ f: ",
              dumpf(mapat($_, \&f, $mats)), $/;
    }
    # mapat 0 f: f([[[1,2],[3,4]],[[5,6],[7,8]]])
    # mapat 1 f: [f([[1,2],[3,4]]),f([[5,6],[7,8]])]
    # mapat 2 f: [[f([1,2]),f([3,4])],[f([5,6]),f([7,8])]]
    # mapat 3 f: [[[f(1),f(2)],[f(3),f(4)]],[[f(5),f(6)],[f(7),f(8)]]]

##</code><code>##

    # Dimensions
    #
    #   (r,c)        <=== transpose ===>    (c,r)
    #
    # Data
    #
    # [ [ 1, 2 ],                         [ [ 1, 3, 5 ], 
    #   [ 3, 4 ],    <=== transpose ===>    [ 2, 4, 6 ] ]
    #   [ 5, 6 ] ] 
    # 

##</code><code>##

    sub zip {
        my $len = min( map scalar @$_, @_ );
        map [ do { my $i=$_; map $_->[$i], @_ } ], 0..$len-1;
    }

    sub transpose {
        map [ zip(@$_) ], @_;
    }

##</code><code>##

    #  0 1 2      0 1 2
    # (m,r,c) -> (r,m,c)  # via transpose(...)
    #
    # ([[[1,2],[3,4]],[[5,6],[7,8]]]
    #         -> [[[1,2],[5,6]],[[3,4],[7,8]]]

##</code><code>##

    #  0 1 2      0 1 2
    # (r,m,c) -> (r,c,m)  # via mapat(1, \&transpose, ...)
    #
    # [[[1,2],[5,6]],[[3,4],[7,8]]]
    #         -> [[[1,5],[2,6]],[[3,7],[4,8]]]

##</code><code>##

    mapat(2, \&f, [[[1,5],[2,6]],[[3,7],[4,8]]]);
    # [[f([1,5]),f([2,6])],[f([3,7]),f([4,8])]]

##</code><code>##

    sub matrixwise_map(&@) {
        my $f = shift;
        # force array context since we're always dealing w/ matricies
        ( mapat( 2, $f, mapat(1, \&transpose, transpose(@_)) ) )[0];
    }

##</code><code>##

    $result = matrixwise_map(\&f, $mats);

    print dumpf($result), $/;
    # [[f([1,5]),f([2,6])],[f([3,7]),f([4,8])]]

##</code><code>##

    use List::Util qw( sum );
    sub mean { sum(@_) / @_ }

    print dumpf( matrixwise_map {mean(@$_)} $mats ), $/;
    # [[3,4],[5,6]]