Everything I wrote in the previous answer you cite still holds: doable, but possibly painful as you may need to write converters. Contributions would still be welcome.
For what it is worth, I use the (Rcpp)Armadillo containers for three-dimensional data as they do have the slicing operators. Note that you can't easily convert them to something R likes ,ie I think we still automated converters for cube
to lists of matrices.
Edit: For what it is worth, here is a short loop from a recent GitHub project of mine:
for (unsigned int j=k-1-1; j>0; j--) {
arma::mat Ppred = AA.slice(j) * P.slice(j) * AA.slice(j).t() + QQ.slice(j);
arma::mat lhs = (P.slice(j) * AA.slice(j).t());
arma::mat rhs = Ppred;
D.slice(j) = arma::solve(rhs.t(), lhs.t()).t();
M.col(j) = M.col(j) + D.slice(j) * (M.col(j+1) - AA.slice(j) * M.col(j));
P.slice(j) = P.slice(j) + D.slice(j) *
(P.slice(j+1) - Ppred) * D.slice(j).t();
}
This uses Armadillo slicing on both the left and right-hand sides. And this works rather well from R thanks to RcppArmadillo (modulo the aforementioned issue that because a R has no real native 3-d structure, so we can't pass a 3-d matrix back easily).