title: Matrix Multiplication along specified array dimension with R/Rcpp
body: <p>Given an <code>n</code> dimensional array <code>X</code>, a <code>d</code> by <code>d-1</code> dimensional matrix <code>V</code> and two specified dimensions <code>(p1, p2) &lt;= (n, n)</code>; I would like a function that preforms matrix multiplication of <code>V</code> along the dimensions <code>(p1, p2)</code> of <code>X</code>. </p>
  
  <p>That is given <code>X</code>:</p>
  
  <pre><code>library(abind)
  set.seed(4)
  X &lt;- matrix(runif(4), 2, 2)
  X &lt;- abind(x, x+5, along = 3)
  &gt; a
  , , 1
  
       [,1] [,2]
  [1,]    1    3
  [2,]    2    4
  
  , , 2
  
       [,1] [,2]
  [1,]    6    8
  [2,]    7    9
  </code></pre>
  
  <p>and given a matrix <code>V</code></p>
  
  <pre><code>V &lt;- matrix(c(1, 2))
        [,1]
  [1,]    1
  [2,]    2
  </code></pre>
  
  <p>For example, if <code>p1=2</code> and <code>p2=1</code>  I would like to remove the following for loop</p>
  
  <pre><code>p1 &lt;- 1
  p2 &lt;- 2
  a.out &lt;- array(0, c(2, 1, 2))
  for (i in 1:dim(a)[2]){
    a.out[,,i] &lt;- a[,,i]%*%V # note indexed along other dimension 
  }
  
  &gt; a.out
  , , 1
  
       [,1]
  [1,]    7
  [2,]   10
  
  , , 2
  
       [,1]
  [1,]   22
  [2,]   25
  </code></pre>
  
  <p>The hard part here is that I want to allow for arbitrary dimensional arrays (i.e., <code>n</code> could be greater than 3). </p>
  
  <p><strong>1st Edit:</strong>
  This problem is not the same as <a href="https://stackoverflow.com/questions/23378980/indexing-slice-from-3d-rcpp-numericvector">Indexing slice from 3D Rcpp NumericVector</a> as I am discussing arbitrary number of dimensions >=2 and the question is not only about indexing. </p>
  
  <p><strong>2nd Edit:</strong>
  Just to be a little more clear here is another example of what I am trying to do.  Here the dimension of <code>X</code> is 4, <code>p1 = 2</code>, <code>p3=3</code>, and the dimension of <code>X</code> along the <code>p1</code> dimension is 12. The following code computes the desired result as <code>X.out</code> for random <code>X</code> and <code>V</code>. </p>
  
  <pre><code>X &lt;- array(rnorm(672), c(4, 7, 12, 2))
  V &lt;- matrix(rnorm(132), 12, 11) # p1 = 2, p2 = 3, V is of dimension D x D-1
  d &lt;- dim(X)
  X.out &lt;- array(0, dim=c(d[1:2], d[3]-1, d[4]))
  for(i in 1:d[1]){
    for (j in 1:d[4]){
      X.out[i,,,j] &lt;- X[i,,,j]%*%V # p1 = 2, p2 = 3
    }
  }
  </code></pre>