如何将Rcpp矩阵覆写逻辑转换为cpp11实现？

原有Rcpp代码可实现矩阵元素缩放后，将上三角元素复制到对应下三角位置：

#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;
}

在R中调用验证功能：

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

转换为cpp11代码时，直接沿用Rcpp的赋值逻辑会报错（赋值操作被隐式删除），尝试的cpp11代码如下：

#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;
}

解决方案 #

cpp11的矩阵类设计与Rcpp不同：operator()仅用于读取元素（返回const引用），修改元素必须通过set(row, col, value)成员方法完成。修正后的代码如下：

#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) {
            // 使用set()方法完成下三角元素赋值
            M.set(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) {
            // 使用set()方法完成元素缩放
            M.set(i, j, M(i, j) * factor);
        }
    }

    invert_tri(M, K);

    return M;
}

说明 #

• 对于writable::doubles_matrix<>类型，所有元素修改操作必须调用set()方法，传入行索引、列索引和目标值
• 读取元素仍可使用M(i,j)的方式，无需修改
• 修正后代码功能与原Rcpp代码完全一致，在R中调用方式不变

内容的提问来源于stack exchange，提问作者pachadotdev