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.

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.

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

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.