Doing %o% that operates on columns instead of atomics

Maybe something like the code below?
There is surely a more elegant way of handling the indices, maybe with plyr?

           -s

xx <- array(1:6,c(3,2)); xx

[,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

yy <- array(1:12,c(3,4)); yy

[,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

array(

sapply(1:ncol(xx),
             function(i)
               sapply(1:ncol(yy),
                      function(j)
                        xx[,i]*yy[,j])),
      c(nrow(xx),ncol(xx),ncol(yy))
      )

, , 1

     [,1] [,2]
[1,]    1    4
[2,]    4   10
[3,]    9   18

, , 2

     [,1] [,2]
[1,]    7   10
[2,]   16   22
[3,]   27   36

, , 3

     [,1] [,2]
[1,]    4   16
[2,]   10   25
[3,]   18   36

, , 4

     [,1] [,2]
[1,]   28   40
[2,]   40   55
[3,]   54   72

xx

[,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6