ILU preconditioner error

[ZhuonanLin]

Hi I am trying to solve a complex number system with the ILU0 preconditioner (single level), the declaration is

typedef amgcl::backend::builtin< std::complex<data_type> > Backend;
	typedef amgcl::make_solver<
		amgcl::relaxation::as_preconditioner<Backend, relaxation::ilu0>,
		amgcl::solver::fgmres<Backend>
	> Solver;

this->solve = new Solver(std::tie(mat_size, preconditioner_ia, 
                                 preconditioner_ja, preconditioner_a), prm);
std::tie(iters, error_out) = (*(this->solve))(*this, x, y);

the preconditioner_ia, preconditioner_ja, preconditioner_a are the row, column, and value array of the crs format preconditioner matrix P. I make sure that diagonal elements of P is non-zero. When solving the system, in some case this throws the error
`
Zero pivot in ILU

No diagonal value in system matrix
`

I tracked the error message in the amgcl source codes, the error message comes from the codes in ilu0.hpp:

for(ptrdiff_t j = row_beg; j < row_end; ++j) {
                ptrdiff_t c = A.col[j];

               // Exit if diagonal is reached
                if (c >= i) {
                    precondition(c == i, "No diagonal value in system matrix");
                    precondition(!math::is_zero((*D)[i]), "Zero pivot in ILU");

                    (*D)[i] = math::inverse((*D)[i]);
                    break;
                }
             ...

 }

It seems to me that in the preconditioner matrix P, it cannot contains element that has column number c greater than the row number i . Am I correct? And what should I do to solve this issue? Should I modify the P to remove the element with c>i? Thanks in advance for any help!

[ddemidov]

This error usually means what it says: the system matrix A has zero somewhere on the diagonal. The ILU implementation in amgcl does not have pivoting, so it can not cope with this situation.

[ZhuonanLin]

Thanks for the reply. For the system matrix A I use a matrix-vector product instead of a matrix directly, since my system matrix A is dense. Is there any method to deal with it?

[ddemidov]

How are you using amgcl then? If you set up ILU0 for a matrix P that is approximating your system matrix A, then your P is the "system matrix" for amgcl.

[ZhuonanLin]

I overwrite the spmv_impl and residual_impl method

namespace amgcl {
	namespace backend {
		template <class Alpha, class DATA, class Vector1, class Beta, class Vector2>
		struct spmv_impl<Alpha, Linear_Solver<DATA>, Vector1, Beta, Vector2> {
			static void apply(Alpha alpha, const Linear_Solver<DATA> &A, const Vector1 &x, Beta beta, Vector2 &y)
			{
				vector<arcomplex<DATA>> t(y.size());
				A.linear_solver_lhs(&t[0], &x[0]);
				
				//// This should implement operation y = beta * y + alpha * A * x
				//// you can probably adjust implementation of your `lhs_gmres_solver_complex()` function,
				//// but right now you would need a temporary vector to do this:
				//std::vector<std::complex<T>> t(y.size());
				for (int i = 0; i < y.size(); ++i) {
					y[i] = beta * y[i] + alpha * t[i];
				}
			}
		};


		template <class DATA, class Vector1, class Vector2, class Vector3>
		struct residual_impl<
			Linear_Solver<DATA>, Vector1, Vector2, Vector3
		>
		{
			static void apply(
				const Vector1 &rhs,
				const LLG_Linear_Solver<DATA> &A,
				const Vector2 &x,
				Vector3       &res
			)
			{
				A.linear_solver_lhs(&res[0], &x[0]);
				for (int p = 0; p < x.size(); p++)
				{
					res[p] = rhs[p] - res[p];
				}
			}
		};
	} // namespace backend
} //namespace amgcl

where Linear_solver is my own class, and Linear_solver.linear_solver_lhs calculates the matrix-vector product Ax. And then use single level ILU0 for the solver, as described above. (I will write it again below):

typedef amgcl::backend::builtin< std::complex<data_type> > Backend;
	typedef amgcl::make_solver<
		amgcl::relaxation::as_preconditioner<Backend, relaxation::ilu0>,
		amgcl::solver::fgmres<Backend>
	> Solver;

this->solve = new Solver(std::tie(mat_size, preconditioner_ia, 
                                 preconditioner_ja, preconditioner_a), prm);
std::tie(iters, error_out) = (*(this->solve))(*this, x, y);

[ddemidov]

Ok, I see (sorry for not getting it from the first time). So it looks like your preconditioner matrix (std::tie(mat_size, preconditioner_ia, preconditioner_ja, preconditioner_a) ) has zero somewhere on a diagonal.

[ZhuonanLin]

Thanks for your reply. I have double-checked my preconditioner matrix P, it does have a problem on the diagonal matrix arrangement, and after reorganizing the elements it works now.

However, when I switched to the ILUt preconditioner, i.e. I use

typedef amgcl::backend::builtin< std::complex<data_type> > Backend;
	typedef amgcl::make_solver<
		amgcl::relaxation::as_preconditioner<Backend, relaxation::ilut>,
		amgcl::solver::fgmres<Backend>
	> Solver;

all other codes remain same. In this case, the preconditioner part doesn't fail but the solver fails, it quits in the first iteration with error = nan.

I also trited ILUk factorization, and it works. Both ILUt and ILUk are used with default prm. Do you have any quick idea why the ILUt fails? Thanks~

[ddemidov]

The ILUt is somewhat tricky, as it often requires manual tuning of its tau and p parameters for a specific problem (and sometimes it is hard to find the working ones). Try reducing tau from its default value of 1e-2 to 1e-4, 1e-6, etc. And increasing p to 5, 10, etc. See https://amgcl.readthedocs.io/en/latest/components/relaxation.html#ilut