有没有办法使用矩阵的非连续视图

即像

arma::mat A = arma::zeros(3,3);
arma::uvec idx = {0,2};
A(idx,idx) += 2;


但是使用矩阵A的子视图?



arma::subview<double> swA = A.submat(0,0,2,2);
swA(idx,idx) += 2.5;


这最后一点不编译为

 error: no match for call to ‘(arma::subview<double>) (arma::uvec&, arma::uvec&)’
swA(idx,idx) += 2.5;


只是提供一些上下文信息,可以帮助我将arma::subviewc用作函数的参数。由于A.submat(0,0,2,2)rtype,所以我不能将其作为非常量参数传递给函数,并且在函数内部,我需要能够以非连续的方式访问元素。

一个最小的(非工作)示例也明白我的意思可能是以下内容

#include <iostream>
#include <armadillo>

void f(arma::subview<double> x)
{
  arma::uvec i = {0,2};
  arma::uvec j = {1,2};
  x(i,j)  += 2.5;
}

int main ()
{

  arma::mat A = arma::zeros(5,5);
  std::cout << A << std::endl;

  f(A.submat(0,0,2,2));
  std::cout << A << std::endl;

  return 0;
}


gcc再次返回

 error: no match for call to ‘(arma::subview<double>) (arma::uvec&, arma::uvec&)’


解决此问题的愚蠢做法是复制矩阵,将其作为参考传递,然后复制回A:

#include <iostream>
#include <armadillo>

void f(arma::mat& x)
{
  arma::uvec i = {0,2};
  arma::uvec j = {1,2};
  x(i,j)  += 2.5;
}

int main ()
{

  arma::mat A = arma::zeros(5,5);
  std::cout << A << std::endl;

  arma::mat B = A.submat(0,0,2,2);
  f(B);
  A.submat(0,0,2,2) = B;
  std::cout << A << std::endl;

  return 0;
}


可以编译并完美运行,但是我需要避免不惜一切代价复制矩阵(A比5x5大得多!)

再次明确,我无法执行以下操作

[...]
void f(arma::mat& x)
{
  arma::uvec i = {0,2};
  arma::uvec j = {1,2};
  x(i,j)  += 2.5;
}
[...]
f(A.submat(0,0,2,2));


因为子视图将是rtype,而我将从gcc获取

invalid initialization of non-const reference of type ‘arma::mat& {aka arma::Mat<double>&}’ from an rvalue of type ‘arma::mat {aka arma::Mat<double>}


我有麻烦吗(只是一个实现问题,不在开发人员的TODO列表中),还是比我聪明的人有一个优雅的解决方案?

谢谢!

旁注:


如果在那里做这样的事情微不足道,我很乐意将线性代数库更改为其他库(例如Eigen或其他),但是我宁愿不这样做,因为我已经使用armadillo多年了,而且我一直非常满意...
我知道在我展示的简单代码中使用循环和更简单的子视图获得相同结果的可能性,但是实际代码更加复杂,并且此子视图将用于矩阵运算,因此我必须循环并将子矩阵复制到一个临时对象中,我想避免

最佳答案

当您在矩阵上调用方法subview时,他/她/它使用了不同的类(arma::matrix<T >operator())....并且您发现(还有我来这里回答问题的我),这是一个Rvalue,因此它可以未被引用...但是类arma::matrix<T>的元素可以由类子视图的对象初始化(所以为什么他们不创建将子视图用作左值的访问函数?为什么还要为其创建类子视图? )...我认为他们可以使用一些模板和Move()来解决此问题,但从底层继承到诸如arma::mat之类的较高层次。
(您要求的问题不仅发生在稀疏矩阵上,甚至发生在稠密矩阵上,仅当您输入f(const arma:mat& )但不能对其进行修改时,该问题仍然有效)

解决方案:什么都没有,我尝试使用std :: move()将子视图的右值欺骗为一个空矩阵,将执行这3个循环的时间基准化为100000次。

/ ************************************ 1 ***************** ***** /

void f(arma::mat &B)
{
}

// [[Rcpp::export]]
int func()
{
  arma::mat A = arma::zeros(50,50);
  //std::cout <<"A:\n"<< A << std::endl;

  arma::uvec i = arma::regspace<arma::uvec>( 0, 1, 10);
  arma::uvec j = arma::regspace<arma::uvec>( 0, 1, 10);
  unsigned p=1;
  do{
    arma::mat B=A(i,j);
    f(B);
    A(i,j)= B;
    p++;
  }while(p<100000);
  return 0;
}


/ ************************** 2 ********************** * /

void f(arma::mat &B)
{
}

// [[Rcpp::export]]
int func2()
{
  arma::mat A = arma::zeros(50,50);
  arma::uvec i = arma::regspace<arma::uvec>( 0, 1, 10);
  arma::uvec j = arma::regspace<arma::uvec>( 0, 1, 10);
  unsigned p=1;
  arma::mat B;
  do{
    B=A(i,j);
    f(B);
    A(i,j) = B;
    p++;
  }while(p<100000);
  return 0;
}


/ ************************ 3 ********************** /

void f(arma::mat &B)
{
}

// [[Rcpp::export]]
int func3()
{
  arma::mat A = arma::zeros(50,50);
  arma::uvec i = arma::regspace<arma::uvec>( 0, 1, 10);
  arma::uvec j = arma::regspace<arma::uvec>( 0, 1, 10);
  unsigned p=1;
  arma::mat B;
  do{
    B=std::move(A(i,j));
    f(B);
    A(i,j)= std::move(B);
    p++;
  }while(p<100000);
  return 0;
}


从R编译和启动显示与构造函数copy或std :: move()没有区别。即使arma::mat B(std::move(A(i,j))没有差异也尝试过,时间与func()相同(由于变量声明,该时间略高于func2()func3())。

  Unit: milliseconds
    expr      min       lq     mean   median       uq      max neval cld
  func() 20.56159 20.68200 21.22679 20.98329 21.55746 27.62831  1000   b
 func2() 19.21450 19.37398 19.88782 19.66397 20.14731 27.86020  1000  a
 func3() 19.19892 19.37547 19.89098 19.67835 20.17973 27.19525  1000  a


这似乎是最快的(以前,A(i,j)= std :: move(B)将B的大小设为0,并且当B = A(i,j)时需要重新分配):

arma::mat B;
  do{
    B=std::move(A(i,j));
    f(B);
    A(i,j)= B;
    p++;
  }while(p<100000);
  return 0;


相同的基准

Unit: milliseconds
    expr      min       lq    mean   median       uq      max neval
 func4() 19.22223 19.30687 19.3907 19.33198 19.36293 23.20771  1000


std::move(A(i,j))似乎不会影响矩阵A(不会松散其内容),因为A(i,j)是子视图的r值,而不是矩阵本身。

关于c++ - Armadillo :arma::subview的非连续 View ,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/41123589/

10-12 23:53