Reply to topic

UmfPack/SuperLU wrong result if LHS is real()/imag()

phmarti
Registered Member
Posts
5
Karma
0
OS
Hi,

Using the UmfPackSupport or SuperLUSupport module from Eigen 3.2.4. The code below compiles but does not provide the correct solution. Both modules access the output memory directly and not only "real()" or "imag()". SparseLU gives the correct solution. If this is the expected behaviour is it possible to make the compilation fail as that leads to very hard to find bugs?

Code: Select all
#include <Eigen/UmfPackSupport>
#include <Eigen/SparseLU>
#include <Eigen/OrderingMethods>
#include <Eigen/SuperLUSupport>
#include <vector>
#include <iostream>
#include <complex>

int main()
{
   Eigen::UmfPackLU<Eigen::SparseMatrix<double> > solver;
//   Eigen::SuperLU<Eigen::SparseMatrix<double> > solver;
//   Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<Eigen::SparseMatrix<double>::Index> > solver;

   Eigen::SparseMatrix<double>   mat(5,5);
   typedef Eigen::Triplet<double> T;
   std::vector<T> triplets;

   triplets.push_back(T(0,0, 0.1667));
   triplets.push_back(T(1,1, 0.5000));
   triplets.push_back(T(0,2,-1.6000));
   triplets.push_back(T(2,2, 0.6000));
   triplets.push_back(T(1,3,-1.6667));
   triplets.push_back(T(3,3, 0.6667));
   triplets.push_back(T(0,4, 1.0000));
   triplets.push_back(T(2,4,-1.7143));
   triplets.push_back(T(4,4, 0.7143));

   mat.setFromTriplets(triplets.begin(), triplets.end());

   solver.compute(mat);
   Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> lhs(5);
   Eigen::Matrix<double, Eigen::Dynamic, 1> rhs(5);

   lhs.real().setConstant(-42);
   lhs.imag().setConstant(-42);

   rhs.setConstant(1);
   lhs.real() = solver.solve(rhs);
   rhs.setConstant(2);
   lhs.imag() = solver.solve(rhs);

   std::cerr << lhs.transpose() << std::endl;
}


regards,
Philippe
User avatar ggael
Moderator
Posts
3447
Karma
19
OS
Fixed in devel branch:

https://bitbucket.org/eigen/eigen/commits/97797a4c10d4/
Changeset: 97797a4c10d4
User: ggael
Date: 2015-02-03 22:46:05+00:00
Summary: Use Ref<> to ensure that both x and b in Ax=b are compatible with Umfpack/SuperLU expectations

We also need to fix Cholmod and SPQR.

 
Reply to topic

Bookmarks



Who is online

Registered users: Bing [Bot], Google [Bot]