Product brakes UnitX Modes in nested triangularView/selfadjointView expressions of a scalar_constant_op * Matrix
Submitted by Patrick Peltzer
Assigned to Nobody
Link to original bugzilla bug (#1517)
Version: 3.3 (current stable)
Operating system: Linux
Description
Created attachment 826
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
results:
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() = (s1mat.adjoint()).template triangularView<Mode>() * (s2ge_left.transpose()), s1triTr.conjugate() * (s2ge_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.
Attachment 826, "Basic code example illustrating the bug":
triangularView_bug.cpp