How to use ilu preconditioner correctly?

[zjvskobe]

I want to utilize Ginkgo's restarted GMRES with an ilu preconditioner as a solver for sparse linear systems. But there are no detailed instructions provided in the tutorial. I implemented the GMRES+ilu preconditioner with reference to the example, but the result of the solution is nan. Here is my code snippet:

using ValueType =double;
using IndexType = int;

// create matrix and vector
auto A = gko::share(gko::matrix::Coo<ValueType , IndexType >::create(exec, gko::dim<2>{size,size}, nEntries));
auto x = gko::matrix::Dense<ValueType >::create(exec, gko::dim<2>(size,1));
auto b = gko::matrix::Dense<ValueType >::create(exec, gko::dim<2>(size,1));
// fill data
double * xdata = x->get_values();
double * bdata = b->get_values();
double *val = A->get_values();
int *r = A->get_row_idxs();
int *c = A->get_col_idx();
for(int i = 0; i < size; i++)
{
     /* fill x,and b*/
}
for(int i = 0; i < nEntries; i++)
{
     /* fill r, c and val*/
}
/// next is from the example:
auto par_ilu_fact = gko::factorization::ParIlu<ValueType, IndexType>::build().on(exec);
auto par_ilu = gko::share(par_ilu_fact->generate(A));
auto ilu_pre_factory = gko::preconditioner::Ilu<gko::solver::LowerTrs<ValueType, IndexType>,
                                 gko::solver::UpperTrs<ValueType, IndexType>,
                                 false>::build().on(exec);
auto ilu_preconditioner = gko::share(ilu_pre_factory->generate(par_ilu));
ValueType reduction_factor = 1e-2;
auto ilu_gmres_factory =gmres::build().with_criteria(gko::stop::Iteration::build().with_max_iters(1000u).on(exec),
                gko::stop::RelativeResidualNorm<ValueType>::build().with_reduction_factor(reduction_factor).on(exec))
            .with_krylov_dim(30u)
            .with_generated_preconditioner(ilu_preconditioner)
            .on(exec);
auto ilu_gmres = ilu_gmres_factory->generate(A);

// Solve system
ilu_gmres->apply(gko::lend(b), gko::lend(x));

The assembly of the matrix and the right-hand-side vector should be no problem because when switching to the Jacobi preconditioner I get the correct result:

using bj = gko::preconditioner::Jacobi<ValueType, IndexType>;
auto ilu_gmres_factory =gmres::build().with_criteria(gko::stop::Iteration::build().with_max_iters(1000u).on(exec),
                gko::stop::RelativeResidualNorm<ValueType>::build().with_reduction_factor(reduction_factor).on(exec))
            .with_krylov_dim(30u)
            .with_preconditioner(bj::build().with_max_block_size(1u).on(exec))
            .on(exec);
auto ilu_gmres = ilu_gmres_factory->generate(A);

// Solve system
ilu_gmres->apply(gko::lend(b), gko::lend(x));

What are the problems with my way of using GMRES with the ilu-preconditioner?

[upsj]

In case you are using Reference, CUDA or HIP, you could try switching from ParIlu (which incurs a certain approximation error), to Ilu (which is exact). Can you give some background on your matrix? There are certain cases where ParIlu may generate a rather poor preconditioner. We need to improve our handling of diverging solvers anyways.

[zjvskobe]

My application background implements a CFD solver that solves the Navier-Stokes equations. And it uses an unstructured grid, so the matrix may be ill-conditioned. Could you please tell me how to use the correct ilu preconditioner? The relevant documentation and examples are not particularly detailed.

[upsj]

That is true, our documentation needs some expansion, we will focus on that before the next release. It should be sufficient to just replace ParIlu<...> by Ilu<...> if you included ginkgo.hpp directly, including ginkgo/core/factorization/ilu.hpp otherwise.

[zjvskobe]

I replaced ParILU with ILU and it still doesn't work. I have a doubt if there is any effect of using the matrix of COO type? I opened "ginkgo/core/factorization/ilu.hpp" and I found that the matrix type is fixed to Csr:

template <typename ValueType = gko::default_precision,
          typename IndexType = gko::int32>
class Ilu : public Composition<ValueType> {
public:
    using value_type = ValueType;
    using index_type = IndexType;
    using matrix_type = matrix::Csr<ValueType, IndexType>;

[upsj]

You are using it correctly, the matrix automatically gets converted into Csr internally to set up the preconditioner. Would it be possible to share the matrix? You can run gko::write_binary(output_stream, A.get()) to write the matrix into a std::ofstream. There are a couple ways that a GMRES solver can break down, I'd like to investigate further what is happening in your case :) Just to be sure: Are you making sure that the matrix data you pass to your Coo matrix is sorted/grouped by row index? Otherwise, weird things might happen.

[zjvskobe]

Matrices A and b (suffixed with .dat) are provided in the attachment. As you said, my Coo matrix is most likely not sorted/grouped by row index, since it's filled with 5x5 submatrix of blocks. So I also provide the row, column, and value information (A.txt) of the Coo matrix in the attachment for your reference.
matrix.zip

[upsj]

I guess we need to push #770 forward a bit to be able to figure out such issues quickly - passing unsorted data to Coo may cause all kinds of weird things to happen, depending on the executor. Did you ever try running your code with the Reference executor? Does it fail there as well? Reference doesn't have any issues with unordered Coo data, so we can use it for comparison

Also, I can't seem to access the uploaded file (Not Found), could you try attaching it again?

[zjvskobe]

The executed is created by gko::ReferenceExecutor::create() so it is exactly the Reference executor. Now let me try to upload the file again.
matrix.zip

[greole]

Could you elaborate on how the matrix is generated, ie. which CFD method is behind it? Could you try a different solver - preconditioner combination like BiCGStab + Jacobi, just to see if it also produces nans. What I find interesting is that there are off-diagonal elements within the blocks, which are much larger compared to main diagonal entries, however, if that is the cause of the convergence issues is unclear.

Additionally, A.dat and A.txt seem to contain different values, or at least different ordering A[0,1] = -0.00314816 in A.txt, vs -0.23123 in A.dat. The value -0.00314816 appears much later A.dat, so could it be that the ordering is broken?

Furthermore, looking while the sparsity pattern of A.txt looks ok, the one presented in A.dat looks broken. Can you double-check that the conversion is correct.

[greole]

I'll close this issue now, as it seems stale. Feel free to reopen if needed.