Ok, even though it is discutable whether it make sense to allow an outer-stride smaller than the inner-size for matrices (e.g., BLAS does not allow that), I've still worked around the issue in both 3.2 and devel branches.
Everything should be fine now, thanks for the detailed report.
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:
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:
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
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