[eigen] solve_complex_linear_system

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]


Hi Eigen experts,
I am part of a developer team to evaluate Eigen's capabilities to solve a system of complex linear equation (such as A * x = b).
Our test case are small self-adjoint matrices such as (here as matlab code):
A=[-1 1-2i 0; 1+2i 1 -i; 0 i 2];
b=[1 2+i 3];
Eigen's solution is:
x_eigen=[0.9245 - 0.7358i 0.6792 + 0.6226i 0.7358 - 1.0755i];
Matlab's solution is correct though: 
x_matlab=[0.4545 + 0.6364i 0.5455 - 0.4545i 1.7273 + 0.2727i];

We are using Eigen's Conjugate Gradient or BiConjugate Gradient algorithm to solve this, e.g.:

template <class ColumnMatrixType, template <typename> class SolverType>
class SolveLinearSystemAlgorithmEigenCGImpl
{
public:
SolveLinearSystemAlgorithmEigenCGImpl(SharedPointer<ColumnMatrixType> rhs, double tolerance, int maxIterations) :
tolerance_(tolerance), maxIterations_(maxIterations), rhs_(rhs) {}

using SolutionType = ColumnMatrixType;

template <class MatrixType>
typename ColumnMatrixType::EigenBase solveWithEigen(const MatrixType& lhs)
{
SolverType<typename MatrixType::EigenBase> solver;
solver.compute(lhs);

if (solver.info() != Eigen::Success)
BOOST_THROW_EXCEPTION(AlgorithmInputException()
<< LinearAlgebraErrorMessage("Eigen solver initialization was unsuccessful")
<< EigenComputationInfo(solver.info()));

solver.setTolerance(tolerance_);
solver.setMaxIterations(maxIterations_);
auto solution = solver.solve(*rhs_).eval();
tolerance_ = solver.error();
maxIterations_ = solver.iterations();
return solution;
}

double tolerance_;
int maxIterations_;
private:
SharedPointer<ColumnMatrixType> rhs_;
}; 

The algorithm runs (e.g., 5000 iterations and target error of 0.00001) quickly but outputs a wrong solution (see above).
We checked obvious things (like that the inputs and outputs gets messed up somehow) indicating that we did a mistake - but that seems not to be the case.
Has anyone experienced that as well?
Is there anything we need to do to make it work (e.g., specific algorithm parameters, specific matrix properties we missed to check for)?
Are there limits what Eigen can do in that regard?

Thanks so much!
Best,
Moritz


Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/