[eigen] Question (bug?) on using TriangularView::solve with maps

[Perry de Valpine <pdevalpine@xxxxxxxxxxxx>]

I’m running into an issue with the following basic setup:

x  = A.triangularView<Eigen::Lower>().solve(b);

A reproducible example is below.

If I use a map for x or b it works.  But for some combinations of maps for x and b, depending on how they are set up, it fails.  In one case it runs but doesn’t put the answer in x.  In another (more dubious) case, it exits for Floating point exception 8.

If x is a map with default treatment of strides, it works for any of several ways of setting up the map for b.  But if x is a map with dynamic strides and b is also a map, then the answer x is not filled in.  I dug as far as seeing that this arises because of the if clause on line 258 of SolveTriangular.h:

template<typename Dest> inline void evalTo(Dest& dst) const

  {

    if(!(is_same<RhsNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_rhs))) {

      dst = m_rhs;

    }

    m_triangularMatrix.template solveInPlace<Side>(dst);

  }

In the cases where the solve works, the if condition evaluates to true, so dat = m_rhs is executed.  In the cases where the answer x is not filled in, the if condition evaluates to false.

There seems to be another separate issue.  In cases where I set up a map representing a one-column matrix, the outer stride should never be needed and it has worked fine to set it to 0 and use the map in other Eigen operations.  However when I do that for the map for x and use as above, I get a Floating point exception 8.  That appears to arise from the line    

 Index subcols = cols>0 ? l2/(4 * sizeof(Scalar) * otherStride) : 0;

in triangular_solve_matrix<…>::run in products/TriangularSolveMatrix.h.  It appears the otherStride is 0 and is being used in the denominator of a size calculation even when not necessary.  I realize one view is that the map should never have been set up with an outer stride of 0, but the habit has been supported by having that work ok in lots of other Eigen operations, so I’m not sure if that is generally intended to be ok or not.  

All of this arises from C++ Eigen code that we generate from the NIMBLE compiler for hierarchical model algorithms from R (r-nimble.org).  Since it is generated code it is not always ideally set up in the way one would hand-code the same idea.

Here is a reproducible example:

void fsTest() {

  Eigen::MatrixXd A(2,2);

  A(0,0) = 2;

  A(1,0) = 1;

  A(0,1) = 0;

  A(1,1) = 3;

  std::cout<<"showing A\n"<<A<<"\n";

  // Here is b as a MatrixXd

  Eigen::MatrixXd b(2, 1);

  b(0) = 1;

  b(1) = 5;

  std::cout<<"showing b\n"<<b<<"\n";

  // Here are two versions of b as maps

  // One issue is using zero as outer stride since it shouldn't be needed

  // Another issue is using the map at all

  vector<double> rawb(2);

  rawb[0] = 1;

  rawb[1] = 5;

  Eigen::Map<MatrixXd, Unaligned, Stride<Dynamic, Dynamic> > mapbVersion1(rawb.data(), 2, 1, Stride<Dynamic, Dynamic>(2, 1));

  Eigen::Map<MatrixXd, Unaligned, Stride<Dynamic, Dynamic> > mapbVersion2(rawb.data(), 2, 1, Stride<Dynamic, Dynamic>(0, 1));

  

  std::cout<<"showing mapbVersion1\n"<<mapbVersion1<<"\n";

  std::cout<<"showing mapbVersion2\n"<<mapbVersion2<<"\n";

  // Regular x

  Eigen::MatrixXd x;

  

 // Three versions of x as maps

  vector<double> rawxFromMaps(2); 

  Eigen::Map<MatrixXd> mapxVersion1(rawxFromMaps.data(), 2,1); // works

  Eigen::Map<MatrixXd, Unaligned, Stride<Dynamic, Dynamic> > mapxVersion2(rawxFromMaps.data(), 2, 1, Stride<Dynamic, Dynamic>(2, 1)); // runs, but gives 0, 0

  Eigen::Map<MatrixXd, Unaligned, Stride<Dynamic, Dynamic> > mapxVersion3(rawxFromMaps.data(), 2, 1, Stride<Dynamic, Dynamic>(0, 1)); //Floating point exception: 8

  // No maps involved

  x = A.triangularView<Eigen::Lower>().solve(b);

  std::cout<<"showing x (no maps)\n"<<x<<"\n"; //works

  mapxVersion1 << 0, 0;

  mapxVersion1 = A.triangularView<Eigen::Lower>().solve(b);

  std::cout<<"showing mapxVersion1 using b\n"<<mapxVersion1<<"\n"; //works

  mapxVersion1 << 0, 0;

  mapxVersion1 = A.triangularView<Eigen::Lower>().solve(mapbVersion1);

  std::cout<<"showing mapxVersion1 using mapbVersion1\n"<<mapxVersion1<<"\n"; //works

  mapxVersion1 << 0, 0;

  mapxVersion1 = A.triangularView<Eigen::Lower>().solve(mapbVersion2);

  std::cout<<"showing mapxVersion1 using mapbVersion2\n"<<mapxVersion1<<"\n"; //works

  mapxVersion2 << 0, 0;

  mapxVersion2 = A.triangularView<Eigen::Lower>().solve(b);

  std::cout<<"showing mapxVersion2 using b\n"<<mapxVersion2<<"\n"; // works

  mapxVersion2 << 0, 0;

  mapxVersion2 = A.triangularView<Eigen::Lower>().solve(mapbVersion1);

  std::cout<<"showing mapxVersion2 using mapbVersion1\n"<<mapxVersion2<<"\n"; // fails: gives 0,0

  mapxVersion2 << 0, 0;

  mapxVersion2 = A.triangularView<Eigen::Lower>().solve(mapbVersion2);

  std::cout<<"showing mapxVersion2 using mapbVersion2\n"<<mapxVersion2<<"\n"; // fails: gives 0,0

  // Any of these give a Floating point exception 8

  std::cout<<"There will be no further output except Floating point exception 8\n";

  mapxVersion3 << 0, 0;

  mapxVersion3 = A.triangularView<Eigen::Lower>().solve(b);

  std::cout<<"showing mapxVersion3 using b\n"<<mapxVersion3<<"\n";

  mapxVersion3 << 0, 0;

  mapxVersion3 = A.triangularView<Eigen::Lower>().solve(mapbVersion1);

  std::cout<<"showing mapxVersion3 using mapbVersion1\n"<<mapxVersion3<<"\n";

  

  mapxVersion3 << 0, 0;

  mapxVersion3 = A.triangularView<Eigen::Lower>().solve(mapbVersion2);

  std::cout<<"showing mapxVersion3 using mapbVersion2\n"<<mapxVersion3<<"\n";

  

}

I’m running on OS X 10.11.4 with Eigen 3.2.8 (eigen-eigen-07105f7124f9) installed a few days ago.

-Perry