Re: [eigen] Sparse non-SPD symmetric solving

[Alberto Luaces <alberto.luaces@xxxxxxxxx>]

Thank you, Gaël.

I was wondering why the SparseLU method was more than three times slower
than SimplicialLDLT for my problem.  Now that you explained
selfadjointView to me, I also realize that SparseLU does not take
symmetry as an advantage, and therefore it is slower, in addition of
Cholesky-based methods being faster in a general basis.

My problem is symmetric, but not generally SPD.  In this test it is, so
I could make a performance comparison between those methods.

Did I read someone talking about a future symmetric SparseLU
implementation in Eigen or is it just my imagination?

Regards,

Alberto

Gael Guennebaud <gael.guennebaud@xxxxxxxxx> writes:

> Hi,
>
> yes, since SparseLU is for general square problems, you need to build
> a full version explicitly using A.selfadjointView<Upper>(). Note that
> Upper is a template parameter, not an argument of the method.
>
> On the other hand, solvers specialized for self-adjoint problems
> (e.g., Simplicial*, ConjugateGradient) read only the triangular part
> as specified by their 2nd template argument. For them, calling
> sefladjointView is thus not required, and even counter productive.
> Here is an example:
>
> ConjugateGradient<SparseMatrix<double>, Upper> cg;
> cg.compute(A);
>
> cheers,
> Gaël
>
> On Fri, Jun 5, 2015 at 10:46 PM, Alberto Luaces
> <alberto.luaces@xxxxxxxxx> wrote:
>
>     Hi,
>     
>     I want to solve a symmetric, non-definite positive sparse linear
>     system. My input is a matrix whose values are stored into the
>     upper
>     part.
>     
>     I wonder if the following is the best way to solve the system. I
>     am
>     concerned about my handling of the symmetric matrix: should I use
>     the
>     adjoint view or try to build an adjoint matrix from the beginning?
>     
>     Eigen::SparseMatrix<double> A;
>     A.setFromTriplets(...);
>     
>     Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
>     
>     solver.isSymmetric(true);
>     solver.compute(A.selfadjointView(Eigen::Upper)).solve(b);
>     
>     Thanks,
>     
>     --
>     Alberto
>     
>     
>