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

On Wed, Apr 13, 2016 at 6:49 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Hi,

I've fixed the shortcoming with the detection of identical input/output (the problem was that extract_data returned 0 for both).

Then, before looking more deeply into the remaining issue, let me first mention that outer-stride==0 for one dimensional objects might not always work in the future, as discussed in the following entry:

http://eigen.tuxfamily.org/bz/show_bug.cgi?id=416

Also, much better use compile-time vectors for 1D matrices -> this will enable much faster evaluation paths.

best,
gael.

On Tue, Apr 12, 2016 at 6:18 PM, Perry de Valpine <pdevalpine@xxxxxxxxxxxx> wrote:
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