This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 356 - Solve() in TriangularSolver.h depending on the storage format
Summary: Solve() in TriangularSolver.h depending on the storage format
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Sparse (show other bugs)
Version: unspecified
Hardware: All All
: --- Unknown
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2011-09-30 11:15 UTC by Donald McGillavry
Modified: 2019-12-04 11:09 UTC (History)
1 user (show)



Attachments

Description Donald McGillavry 2011-09-30 11:15:15 UTC
I wanted to use the triangularView() for sparse matrix and found that it does not work as i need.

Eigen::SparseMatrix<double> m_matrix(m,m);
// Filling matrix: 
// 3 3 0 0 0
// 5 3 3 0 0
// 5 5 3 3 0
// 0 5 5 3 3
// 0 0 5 5 3

// The following functions return incorrect results:
m_matrix.triangularView<Eigen::UnitLower>();
m_matrix.triangularView<Eigen::UnitUpper>();

// UnitLower		Should be
// 3 3 0 0 0		1 0 0 0 0
// 0 3 3 0 0		5 1 0 0 0
// 0 0 3 3 0		5 5 1 0 0
// 0 0 0 3 3		0 5 5 1 0
// 0 0 0 0 3		0 0 5 5 1
// UnitUpper		Should be
// 3 3 0 0 0		1 3 0 0 0
// 0 3 3 0 0		0 1 3 0 0
// 0 0 3 3 0		0 0 1 3 0
// 0 0 0 3 3		0 0 0 1 3
// 0 0 0 0 3		0 0 0 0 1

Please implement this, I myself could not understand how to do it :(
Comment 1 Donald McGillavry 2011-10-11 11:21:07 UTC
I just realized that there is another error.
Function solve () returns different values ​​depending on the storage format of sparse matrices.
For example:
Eigen::SparseMatrix<fp_type,Eigen::RowMajor> Ar(m,m);
Eigen::SparseMatrix<fp_type,Eigen::ColMajor> Ac(m,m);
Eigen::Matrix<double,Eigen::Dynamic,1> b(m);
b.setConstant(5.0);
// Filling matrix Ar and Ac
// 2.0 0.0 3.0
// 0.0 2.0 0.0
// 0.5 0.0 0.5

std::cout<<Ac.triangularView<UnitLower>().solve(b);
//-2.5 5 1.25

std::cout<<Ar.triangularView<UnitLower>().solve(b);
// 5 5 2.5

I think the result should not depend on the format. Maybe I'm wrong?
Comment 2 Donald McGillavry 2011-10-13 08:38:43 UTC
Eigen/src/Sparse/TriangularSolver.h
line 161:
other.coeffRef(i,col) /= lhs.innerVector(i).lastCoeff();

If i==0 lhs.innerVector(i).lastCoeff() returns a value m_matrix (m-1, 0), but must return a value m_matrix (0,0). Because of this, get the wrong value the first element of the solution vector.

This is reflected in
Eigen::SparseMatrix<double,Eigen::ColMajor> m_matrix(m,m);
std::cout<<m_matrix.triangularView<Eigen::Upper>().solve(b);
Comment 3 Donald McGillavry 2011-10-14 08:35:39 UTC
(In reply to comment #2)

It is better to leave the function lhs.coeff(i,i), because the error is shown for each column. This is understandable, because the need to divide by the diagonal element and not the latter. If you want to quickly find the diagonal elements, it is best to store references to them.
Comment 4 Gael Guennebaud 2011-12-03 19:56:52 UTC
With sparse matrices it was assumed that the underlying matrix only holds elements inside the respective triangular part because there is no memory advantage in storing other elements. So what's your use case for that?
Comment 5 Gael Guennebaud 2011-12-03 23:51:32 UTC
The issue with the triangular solver is solved in the change set 0696a3beb7cb using a reverse iterator.
Comment 6 Gael Guennebaud 2011-12-04 14:41:14 UTC
https://bitbucket.org/eigen/eigen/changeset/275a8666a5ce/
changeset:   275a8666a5ce
date:        2011-12-04 14:39:24
summary:     fix bug 356: fix TriangularView::InnerIterator for unit diagonals
Comment 7 Nobody 2019-12-04 11:09:49 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/356.

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