Any way to improve the convergence of ParIlut?

[FlyingCat2023]

I have a matrix which is generated from CutFEM. The matrix is symmetric but may be not positive definite. I try to solve Ax = b by using the ParIlut preconditioner but fail to get a convergent solution. I found that it was because ILU was not converged in both 2D and 3D problems. I wonder if there is a solution to ensure convergence of ParIlut?

[pratikvn]

I believe ILU should be more robust than ParILU/ParILUT. If ILU does not converge, then I am not sure if there is a way to make ParILU/ParILUT converge.

[FlyingCat2023]

Hello, I am late again. Unfortunately, The off-diagonal value are not < 1 by scaling them with the diagonal elements. I found some diagonal elements are very small as compared to the off-diagonal elements around it.
Here are the results.
parilu 2 2 2
Reading ../data/a_1_0.mtx
1;1689553197;1638503874;0;0;
2;3315644637;3236100911;0;0;
parilu 2 2 2
Reading ../data/a_1_0.mtx
1;155339553;192793444;6.52954e-08;6.52954e-08;
2;275994243;296592474;6.52954e-08;6.52954e-08;
parilut 2 2 2
Reading ../data/a_1_0.mtx
1;3837950080;3671037208;6.52954e-08;6.52954e-08;
2;21346279645;22399902701;0.125565;0.125565;
parilut 4 2 2
Reading ../data/a_1_0.mtx
1;3826550646;3740850592;6.52954e-08;6.52954e-08;
2;22383866145;22235150445;0.125565;0.125565;
3;49282028763;51056535073;135.172;513.352;
4;79245004217;79275728756;1325.01;1331.49;

[nbeams]

I found some diagonal elements are very small as compared to the off-diagonal elements around it.

This will definitely cause problems for ParILU(T). You could try a reordering meant to move large values to the diagonal. Ginkgo has the Mc64 reordering algorithm in its experimental namespace: https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1experimental_1_1reorder_1_1Mc64.html

An example of how you could use it:

      // Create the Mc64 object
       auto perm_comp = gko::experimental::reorder::Mc64<ValueType,IndexType>::build()
                   .on(exec)->generate(orig_mtx);
       // Get row/column operators
       auto row_perm_scale = gko::as<gko::matrix::ScaledPermutation<ValueType,IndexType>>(
                      perm_comp->get_operators()[0]);
       auto col_perm_scale = gko::as<gko::matrix::ScaledPermutation<ValueType,IndexType>>(
                      perm_comp->get_operators()[1]);
       // Make a new matrix with Mc64 applied
       mtx_op = gko::share(orig_mtx->scale_permute(row_perm_scale, col_perm_scale));

Which Executor were you using for the results you added? Previously you said "OpenMPExecutor leads to convergence at the first iteration. However, for further iterations, it converges to some extent, but not meeting the tolerance requirements, e.g., 1e-5," which is not what these results show. I also find it very odd for the residual to be 0 after 1 iteration in the first test...

[FlyingCat2023]

Sorry for the confusion.
The first example is using the Reference Executor. The others are using OMP Executor.

[FlyingCat2023]

Test 2 converge to tolerance by using ParILU. However, Test 3 and 4 does not converges by using ParILUT, starting from iteration 3, as can be seen in test 4.

[nbeams]

But the residuals for 1 or 2 iterations are the same in test 2. It's very odd for ParILU to converge after only 1 iteration.
Regardless, if the matrix is not diagonally dominant, then we can't expect that it will converge quickly. We could get lucky -- but we can't guarantee that it should be converging. A reordering may help, or else ParILU(T) may not work for your matrix.