Created attachment 826 [details]
Basic code example illustrating the bug
I'm generally struggling with using the triangularView or selfadjointView classes with any other matrix operation except a matrix product. In the Quick Reference guide, the stated operation of a scalar s1 * m1.adjoint().triangularView<xxx> doesn't compile for me (https://eigen.tuxfamily.org/dox/group__QuickRefPage.html#QuickRef_TriangularView). However, this bug report is not about these features, but only about a product expression of the form
m_id * (0.5*m1).template triangularView<UnitUpper>()
This code does compile and it is also used inside the test system at several points. For simplicity reasons, just assume m_id is the identity matrix. Now, using the UnitUpper (or UnitLower or UnitDiag) mode should result in only diagonal entries with value 1 in the returned matrix. However, the above code will return 0.5 as the diagonal values. This only happens when having an "outer" product expression and when multiplying a constant scalar to m2. In other words, the following lines will return the correct result:
m1 = (0.5 * m).template triangularView<UnitUpper>();
m2 = m_id * (m).template triangularView<UnitUpper>();
m3 = m_id * (m + m).template triangularView<UnitUpper>();
Also, forcing eval() before creating the triangularView will also give correct
m_id * (0.5 * m).eval().template triangularView<UnitUpper>();
As said above, code like this is present in the test system, for example in the product_trmm_21 test, line 54:
VERIFY_IS_APPROX( ge_xs.noalias() = (s1*mat.adjoint()).template triangularView<Mode>() * (s2*ge_left.transpose()), s1*triTr.conjugate() * (s2*ge_left.transpose()));
Note that 'Mode' is 'UnitUpper' in this test case. As you can see, the first argument fullfills the pattern shown above. Therefore, the matrix returned from triangularView will not have ones on the diagonal, but s1 instead. Adding an eval() call to this line will make the test fail. I do not know what the tests are calculating and testing, but regarding the UnitXXX modes it appears to me that the correct way should be to multiply s1 to the return value of the triangularView return matrix, and not to it's source matrix inside the expression. However, moving the s1 multiplication outside is not possible as stated at the beginning of this bug report.
I can confirm that this bug is present on the current development branch as well. You can find a simple code example illustrating the bug as attachment.
I confirm, what happens is that expressions as:
(scalar * A) * B
with A and B expressions formed by a matrix + trivial operations like block/triangularView/conjugate/etc are directly mapped to:
This is wrong with triangularView and a unit diagonal... I'll see how to prevent this optimization.
ah, and the respective unit test suffers from the same shortcoming:
VERIFY_IS_APPROX( ge_xs = (s1*mat).adjoint().template triangularView<Mode>() * ge_left.adjoint(), numext::conj(s1) * triTr.conjugate() * ge_left.adjoint());
which is why we did not catch it.
I fixed it by applying a correction to the result:
that's not super elegant, but that's the only way to avoid the introduction of a temporary without introducing heavy changes. I also believe that a such use case is very seldom.
-- GitLab Migration Automatic Message --
This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.
You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1517.