Advanced constructor for non-contiguous matrices

Question

In my implementations I work a lot with submatrices and blocks of matrices. I am wondering if there is a way in Armadillo that would allow me to extract a block of a larger matrix and use the same memory for this submatrix as for the block within the original matrix. My problem is that I do not know how to go about this as the positions in the original matrix are not contiguous.

Here is one simple example that illustrates what I want to do when my original matrix is A = [A1 A2]:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat foo(arma::mat A, arma::uword block_start, arma::uword block_length) {
  arma::uword rows = A.n_rows;
  arma::mat B = arma::mat(A.begin_col(block_start), rows, block_length, false);
// fill B here for illustration; 
// in the real piece of code I do various manipulations, multiplications, etc using B
  B.fill(1.0);
  return A;
}

/*** R
A <- matrix(0, 4, 4)
foo(A, 0, 2)
> A <- matrix(0, 4, 4)
> foo(A, 0, 2)
     [,1] [,2] [,3] [,4]
[1,]    1    1    0    0
[2,]    1    1    0    0
[3,]    1    1    0    0
[4,]    1    1    0    0
*/

In this case, the positions of the submatrix are contiguous and I can use the advanced constructor to link the memory.

Suppose now that I want the submatrix to be A[1:2, 1:2]. Can I obtain a copy B in Armadillo that uses the same memory as the original elements in A? (Ideally, a solution to this question would also generalize to the case where columns are also non-contiguous, e.g. A[c(1, 3), c(1, 3)].)

Edit: To clarify, I really need the matrix B in the function above to exist on its own. I don't fill it in my real code, but use it (and multiple other submatrices) in various matrix multiplications etc. So what I'm after is a way to create B as a non-contiguous submatrix of A while making sure that they use the same memory.

I think you want a square peg in a round hole. Think of 5x5 matrix. Now you want tthe inner 3x3 ... but you cannot have that contiguously without a copy because _the outer matrix's size determines how the 25 cells in the 5x5 are used_. — Dirk Eddelbuettel, Jun 01 '19 at 02:40
@DirkEddelbuettel You’re probably right. But it’s all about pointers, isn’t it? And everything can be pointed to. Ideally, one would be able to use `submat(rows, cols).begin()` but that does not work. Of course, it might be a lot unsafer because you have a more intricate dependence. I feel a little dirty extracting and copying submatrices time and time again, so avoiding that would be terrific. — hejseb, Jun 01 '19 at 03:52
Hm, I fear you are just repeating what you would like to have while ignoring how memory is used. Start with a notepad and "draw" a 5x5 matrix. Sketch out what the first, second, .., 20th address is then. Then try to find a multi-row / multi-col subview and convince yourself that it cannot a be contiguous sub-vector. I could be wrong, but I still think you are asking for a thing that cannot exist. — Dirk Eddelbuettel, Jun 01 '19 at 16:35
@DirkEddelbuettel I see your point. Just demonstrating the infeasibility of what I’m asking for is valuable to me. There is usually good reasons for features that are lacking and I guess this is one of them. — hejseb, Jun 01 '19 at 20:38

score 1 · Answer 1 · answered May 31 '19 at 10:44

You can use submatrix views to write to contiguous or non-contiguous memory:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat foo(arma::mat A) {
    A.submat(0, 0, 1, 1).fill(1.0);
    return A;
}

// [[Rcpp::export]]
arma::mat bar(arma::mat A) {
    arma::uvec rows;
    rows << 0 << 2;
    arma::uvec cols;
    cols << 0 << 2;
    A.submat(rows, cols).fill(2.0);
    return A;
}

/*** R
A <- matrix(0, 4, 4)
foo(A)
bar(A)
*/

Output:

> A <- matrix(0, 4, 4)

> foo(A)
     [,1] [,2] [,3] [,4]
[1,]    1    1    0    0
[2,]    1    1    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0

> bar(A)
     [,1] [,2] [,3] [,4]
[1,]    2    0    2    0
[2,]    0    0    0    0
[3,]    2    0    2    0
[4,]    0    0    0    0

Thanks, Ralf. These submatrix views are what I use all the time, but maybe I mislead you by making a too minimal example. I do various matrix computations using the submatrices that I extract, so `fill` was just to do something in my example. So I do really need the submatrix to exist as its own object, and what I'm after is to use the same memory for my submatrix as for the original elements in the big matrix. Sorry if I was unclear, I'll update my post to clarify. — hejseb, May 31 '19 at 11:42

Advanced constructor for non-contiguous matrices

1 Answers1