How do I pass persistent references to submatrix views as function arguments?

Question

I am interested in passing subviews of matrices to other objects (or functions) so that the object/function can modify a defined segment of the column/matrix.

I am currently using the armadillo library and I can get and use a reference to a single matrix entry (function f below). I am looking for an effective way to pass around references to subviews. I think the code for function h below accomplishes this goal, but I'm not sure. My main question is: Assuming I want a reference to store in a class (not just use in a function as in h), is there a way to improve on function h below?

Three less important follow up questions are: 1) Can this be accomplished while writing the functions to accept arma::Col types of arguments rather than specifying a subview? 2) If not, is there another C++ matrix library that will allow this sort of thing without making it explicitly about subviews? (I looked at Eigen and its .blocks syntax but I think the issue is the same); and 3) am I right that the reference passed to function g will disappear when the expression completes?

#include <iostream>
#include <armadillo>

void f(double & x) { x = 5.33392;}
void g(arma::subview_col<double> && x, double y) { x(0) = y; }
void h(arma::subview_col<double>    x, double y) { x(0) = y; }

int main () {
  arma::Mat<double> A = arma::randu<arma::Mat<double> >(3,3);
  A(1,1) = 3;
  A.submat(0,0,1,1) = arma::zeros<arma::Mat<double> >(2,2);
  std::cout << A << std::endl;

  double & k = A(1,1);
  std::cout << k << std::endl;

  f(k);         // THis is ok, works.


  g(A.col(2), 153.0);    // It compiles and the value changes...
  // A.col(2) gets destroyed when the above expression is done
  // evaluating... so not a reference one can keep around.

  std::cout << A << std::endl;

  // Here the subview is passed by value, and it keeps its reference to 
  h(A.col(1), 153.0);

  std::cout << A << std::endl;

  return 0;
}

Compiles for me with:

 g++ -o submatrix_views -std=c++11 ./submatrix_views.cpp

Answer 1

I haven't used armadillo so I am not the most appropriate to respond, but I think this a more general C++ issue, plus, there's no answer yet :-)

So:

Assuming I want a reference to store in a class (not just use in a function as in h), is there a way to improve on function h below?

I think you're confusing a couple of basic C++ concepts unrelated to matrices and views. A reference essentially points to an object, so you should first worry about storing the object itself, then the reference.

The fact that there is a reference to an object may be a reason to have the lifetime of the object extended to the lifetime of the reference, but this is not always the case. See this answer to one of my previous questions for some explanation and rules.

So, the best way to be certain is to store the object first. I bet A.col(1) and A.submat(0,0,1,1) are temporary objects and you can store them in variables in any of the following ways:

arma::subview_col<double> col = A.col(1);
decltype(A.col(1)) col = A.col(1);
auto col = A.col(1);

If you know the type, you can simply have one such variable as a member of a class (if you don't, just use decltype). I don't know details of armadillo, but I bet A.col(1) doesn't contain more than a reference to A and an integer (I guess size_t) with value 1. It is a lightweight object designed to behave as a generalized kind of reference. Don't be afraid to store it. Just to practice, try

auto col = A.col(2);
col = 153.0;

any see it's the same as A.col(2) = 153.0; or g(A.col(2), 153.0);. At this point, nothing prevents you from storing col anywhere for later use.

Having the view stored somewhere, you then need to pass it to a function. If the view is to be used for reading only, write the following two overloads:

void read(      arma::subview_col<double> && x, double& y) { y = x(0); }
void read(const arma::subview_col<double> &  x, double& y) { y = x(0); }

If it is used for writing as well, then:

void write(arma::subview_col<double> && x, double y) { x(0) = y; }
void write(arma::subview_col<double> &  x, double y) { x(0) = y; }

My recommended universal solution is the following:

template <typename X> void read (X&& x, double& y) { y = std::forward<X>(x)(0); }
template <typename X> void write(X&& x, double  y) { std::forward<X>(x)(0) = y; }

or, even more generic:

template <typename X, typename Y>
void read (X&& x, Y&& y) { std::forward<Y>(y) = std::forward<X>(x)(0); }

template <typename X, typename Y>
void write(X&& x, Y&& y) { std::forward<X>(x)(0) = std::forward<Y>(y); }

Now X and Y can be anything, or more precisely, X is any kind of array or view, and Y is any type compatible with (convertible to/from) the element type of X. If they are not, the compiler will diagnose it. You could even have an array Y and an array (or view) of arrays X, if armadillo supports that.

Plus, you can generalize to n-ary functions with no worries about infinite const/non-const overload combinations.