Reply to topic

Iterative solvers and ILU pre conditioner

kinetic
Registered Member
Posts
15
Karma
0
Hi, I am trying to use an iterative solver but am running into a few problems. Maybe on my end. Specifically I would like to use the IncompleteLUT class as a preconditioner. But I cannot perform the factorisation in my case

Consider the following pseudo code
Code: Select all
int main() {

  int unx(10000), uny(10000), unz(10000);
  Eigen::SparseMatrix<Complex>    C(unx+uny+unz , unx+uny+unz);
  fillC(C); // C is symmetric non-Hermitian only Upper part is stored. In fact it is all real, except the diagonal.

  // Compute ILUT factorisation
  Eigen::IncompleteLUT< Complex > ILUTC;
  ILUTC.setFillfactor(7);
  ILUTC.compute( C.selfadjointView<Eigen::Upper>() );

}


However this fails with a lot of errors, ending with
/home/tirons/src/eigen/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h:265:69: error: 'const class Eigen::SparseSelfAdjointView<Eigen::SparseMatrix<std::complex<double>, 0, int>, 2u>' has no member named 'nonZeros'

On a somewhat related note. Trying to use the BiCGSTab iterative solver on my symmetric, non-Hermitian matrix does not seem possible without storing both parts.
Code: Select all
BiCGSTAB< SparseMatrix<Complex> > solver;
solver.setTolerance(tol);
A = solver.compute( C.selfadjointView<Eigen::Upper>() ).solve(Se);


Also fails to compile with errors

/home/tirons/src/eigen/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h:103:12: note: no known conversion for argument 1 from 'Eigen::SparseSelfAdjointView<Eigen::SparseMatrix<std::complex<double>, 0, int>, 2u>' to 'const MatrixType& {aka const Eigen::SparseMatrix<std::complex<double>, 0, int>&}'

I can make a copy so that it will be a const
Code: Select all
 SparseMatrix<Complex>  Csym;
        Csym = C.selfadjointView<Eigen::Upper>();
        BiCGSTAB< SparseMatrix<Complex> > solver;
        solver.setTolerance(tol);
        //A = solver.compute( C.selfadjointView<Eigen::Upper>() ).solve(Se);
        A = solver.compute( Csym ).solve(Se);

And this compiles fine. But I'd like to avoid this copy for obvious reasons.

Am I doing something wrong? Is there a way to perform an ILU decomposition on a Symmetric matrix. I can't do a Cholesky due to non-Hermitian nature. I would also like to limit the fill in. This is just for use as a preconditioner.

Thanks for any insight.
kinetic
Registered Member
Posts
15
Karma
0
I should mention I can get around the ILUT decomposition in the same way

Code: Select all
         
        SparseMatrix<Complex>  Csym;
        Csym = C.selfadjointView<Eigen::Upper>();
        Eigen::IncompleteLUT< Complex > ILUTC;
        ILUTC.setFillfactor(7);
        //ILUTC.compute( C.selfadjointView<Eigen::Upper>() );
        ILUTC.compute( Csym );


Works fine. But I'd like to avoid the deep copy and 2X replication.
Dee33
Registered Member
Posts
54
Karma
0
OS
Yes I'm afraid ILUT requires the full matrix.
Note that you don't have to create separately the solver and the preconditioner.
The best way is like
Code: Select all
BiCGSTAB< SparseMatrix<Complex>,IncompleteLUT<double> > solver;
solver.preconditioner().setFillfactor(7);  //Get the reference of the preconditioner and set properties
solver.setTolerance(tol);
SparseMatrix<Complex>  Csym;
Csym = C.selfadjointView<Eigen::Upper>();
solver.compute(Csym) // Compute the ILUT factorization
kinetic
Registered Member
Posts
15
Karma
0
Thanks for the reply.

As I dig into the source code for CongugateGradient and BiCGSTAB I notice that ConjugateGradient actually calls the underlying internal conjugate_gradient with a call like
Code: Select all
      internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
                                   Base::m_preconditioner, m_iterations, m_error);


But I don't really understand the mp_matrix->template selfadjointView<UpLo>() syntax. Is there some way that I could do something similar with BiCGSTAB or IncompleteLUT? Because how ever this works, it allow the selfadjointView to be passed to a function that is expecting a const reference. Which is exactly the problem I am having.

I realize that my use case may be sort of isolated, but it seems like there is no underlying reason why this should not work other than the author of the code thinking (reasonably) that no one using BiCGSTAB would be dealing with a symmetric matrix.

Thanks in advance.

 
Reply to topic

Bookmarks



Who is online

Registered users: ab4bd, Baidu [Spider], BeS, Bing [Bot], Exabot [Bot], GeekQuack, ggael, ghevan, Google [Bot], google01103, igorm, jensreuterberg, kainz.a, maarten, mmistretta, Nuc!eoN, orbmiser, plfiorini, scummos, searchfgold6789, SecretCode, Steve T, TheraHedwig, Uri_Herrera, vblazquez, Yahoo [Bot]