This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1517 - Product brakes UnitX Modes in nested triangularView/selfadjointView expressions of a scalar_constant_op * Matrix
Summary: Product brakes UnitX Modes in nested triangularView/selfadjointView expressio...
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - expression templates (show other bugs)
Version: 3.3 (current stable)
Hardware: x86 - 64-bit Linux
: Normal Wrong Result
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2018-02-08 23:49 UTC by Patrick Peltzer
Modified: 2019-12-04 17:28 UTC (History)
3 users (show)



Attachments
Basic code example illustrating the bug (1.06 KB, text/x-c++src)
2018-02-08 23:49 UTC, Patrick Peltzer
no flags Details

Description Patrick Peltzer 2018-02-08 23:49:11 UTC
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 
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() = (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.
Comment 1 Gael Guennebaud 2018-02-09 10:33:19 UTC
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:

  adequate_kernel_routine(A,B,scalar)

implementing scalar*(A*B).

This is wrong with triangularView and a unit diagonal... I'll see how to prevent this optimization.
Comment 2 Gael Guennebaud 2018-02-09 10:36:35 UTC
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.
Comment 3 Gael Guennebaud 2018-02-09 16:54:53 UTC
I fixed it by applying a correction to the result:

https://bitbucket.org/eigen/eigen/commits/253cc1138131/
https://bitbucket.org/eigen/eigen/commits/a546d43bdd4f/  (3.3)

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.
Comment 4 Nobody 2019-12-04 17:28:15 UTC
-- 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.

Note You need to log in before you can comment on or make changes to this bug.