This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1741 - Map<>.noalias()=A*B gives wrong result
Summary: Map<>.noalias()=A*B gives wrong result
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - matrix products (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Wrong Result
Assignee: Nobody
URL:
Whiteboard:
Keywords: test-needed, wrong-result
Depends on:
Blocks: 3.4
  Show dependency treegraph
 
Reported: 2019-08-16 15:45 UTC by Patrick Peltzer
Modified: 2019-12-04 18:45 UTC (History)
2 users (show)



Attachments
Minimal code example to reproduce the issue (688 bytes, text/x-c++src)
2019-08-16 15:45 UTC, Patrick Peltzer
no flags Details

Description Patrick Peltzer 2019-08-16 15:45:12 UTC
Created attachment 949 [details]
Minimal code example to reproduce the issue

In an alias free multiplication with a Map<> type to assign to, using noalias() can lead to wrong results. Example:

C_Map = A * B;          // Correct
C.setZero();
C_Map.noalias() = A*B;  // Incorrect

See the attached code example on how to reproduce. I've noticed this for larger matrices (C_Map >= 8x8), where only half of the rows are modified and the other half stays untouched. Typical output for the minimal code example:

C_Map
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8

C_Map.noalias()
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

I've reproduced the issue on version 3.3.7 and the current 3.3 and the default branch.
Comment 1 Christoph Hertzberg 2019-08-16 16:29:29 UTC
The information about the inner stride appears to get lost when calling Gemm::run at this line:

https://bitbucket.org/eigen/eigen/src/8071cda5714d7f4/Eigen/src/Core/products/GeneralMatrixMatrix.h#lines-228

The easiest fix I can think of, is to fall back to lazyproduct::eval_dynamic in 
   generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemmProduct>
       ::evalTo, addTo, subTo
if the inner stride of the destination is not 1.

More complicated solutions would require to store packets in a strided way (cf Bug 1737)
Comment 2 Gael Guennebaud 2019-09-10 14:27:36 UTC
We already had all the mechanisms to properly fix this issue:
https://bitbucket.org/eigen/eigen/commits/9037b27e44a5/

However, we also need to fix other products the same way, so I'm not closing it yet.
Comment 3 Gael Guennebaud 2019-09-10 14:36:29 UTC
3.3 backport: https://bitbucket.org/eigen/eigen/commits/cb30738c6762/
Comment 4 Gael Guennebaud 2019-09-10 21:57:26 UTC
https://bitbucket.org/eigen/eigen/commits/509c59338ca4/
Summary:     Bug 1741: fix SelfAdjointView::rankUpdate and product to triangular part for destination with non-trivial inner stride

More to come later...
Comment 5 Gael Guennebaud 2019-09-11 13:15:36 UTC
last chunk of fixes:

https://bitbucket.org/eigen/eigen/commits/df36abf6b401/
Summary:     Bug 1741: fix self-adjoint*matrix, triangular*matrix, and triangular^1*matrix with a destination having a non-trivial inner-stride

and 3.3 backports:

https://bitbucket.org/eigen/eigen/commits/f7f34d725b8b/
https://bitbucket.org/eigen/eigen/commits/4da1c565b22e/
https://bitbucket.org/eigen/eigen/commits/a5c3d12011a0/
Comment 6 Nobody 2019-12-04 18:45:37 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/1741.

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