amgcl

amgcl copied to clipboard

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!

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.

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?

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.

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);

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.

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~

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