Overwriting entries from a matrix in cpp11

Question

I have some Rcpp code to fill entries in a matrix, here's a MWE of what it does:

#include <Rcpp.h>

using namespace Rcpp;

void invert_tri(NumericMatrix &M, int K) {
    for(int i=0 ; i<K ; ++i) {
        for(int j=i+1 ; j<K ; ++j) {
            M(j, i) = M(i, j);
        }
    }
}

// [[Rcpp::export]]
NumericMatrix multiply_entries(NumericMatrix M, int factor) {
    int K = M.ncol();

    for(int i=0 ; i<K ; ++i) {
        for(int j=0 ; j<K ; ++j) {
            M(i, j) = M(i, j) * factor;
        }
    }

    invert_tri(M, K);

    return M;
}

Then in R I would use it like this:

M = matrix(c(1,NA,2,4), nrow = 2, ncol = 2)

multiply_entries(M, 2L)

This changes

     [,1] [,2]
[1,]    1    2
[2,]   NA    4

into

     [,1] [,2]
[1,]    2    4
[2,]    4    8

Then, with cpp11 I would sketch the following code

#include <cpp11.hpp>
#include <cpp11/doubles.hpp>

using namespace cpp11;

void invert_tri(writable::doubles_matrix<> &M, int K) {
    for(int i=0 ; i<K ; ++i) {
        for(int j=i+1 ; j<K ; ++j) {
            M(j, i) = M(i, j);
        }
    }
}

[[cpp11::register]]
doubles_matrix<> multiply_entries(writable::doubles_matrix<> M, int factor) {
    int K = M.ncol();

    for(int i=0 ; i<K ; ++i) {
        for(int j=0 ; j<K ; ++j) {
            M(i, j) = M(i, j) * factor;
        }
    }

    invert_tri(M, K);

    return M;
}

But here doing M(j, i) += M(i, j) would fail and with = fails because of implicit deletion. Do you have any ideas for an alternative approach to translate the working code.

score 1 · Answer 1 · answered Apr 15 '23 at 03:40

The original code fails because the compiler can't resolve which type of assignment to perform in the statement M(j, i) = M(i ,j)

Explanation:

Both the left and right hand side are of type r_vector<double>::proxy.
The right-hand-side value can be cast to a double, for which there is a valid assignment to the left hand side.
The copy and move assignment operators for this type (r_vector<double>proxy>) are both deleted due to const-qualifiers on members.
While you would think that because there is only one non-deleted assignment that there would be no problem - however - deleted assignment operators are still considered during overload resolution.

The simplest solution I can think of is to force the compiler into a corner by casting the right hand side to a double, e.g.

void invert_tri(writable::doubles_matrix<> &M, int K) {
    for(int i=0 ; i<K ; ++i) {
        for(int j=i+1 ; j<K ; ++j) {
            M(j, i) = (double)M(i, j);
        }
    }
}

score 0 · Accepted Answer · answered Dec 26 '22 at 20:18

maybe not the best approach for large matrices, but this just works

void invert_tri(writable::doubles_matrix<> &M, int K) {
    for(int i=0 ; i<K ; ++i) {
        for(int j=i+1 ; j<K ; ++j) {
            M(j, i) += (-1.0 * M(j, i)) + M(i, j);
        }
    }
}

Overwriting entries from a matrix in cpp11

2 Answers2