Ruge–Stüben solver for complex-valued matrix

[learning-chip]

smoothed_aggregation_solver supports complex-valued matrix (without converting to equivalent real form), as shown in this demo, but ruge_stuben_solver seems to take real-valued matrix only:

import pyamg  # v4.2.3
A = pyamg.gallery.gauge_laplacian(100)
ml = pyamg.smoothed_aggregation_solver(A)  # works
pyamg.ruge_stuben_solver(A)  # TypeError

Error message goes like:

TypeError: rs_direct_interpolation_pass2(): incompatible function arguments. The following argument types are supported:
    1. (n_nodes: int, Ap: numpy.ndarray[numpy.int32], Aj: numpy.ndarray[numpy.int32], Ax: numpy.ndarray[numpy.float32], Sp: numpy.ndarray[numpy.int32], Sj: numpy.ndarray[numpy.int32], Sx: numpy.ndarray[numpy.float32], splitting: numpy.ndarray[numpy.int32], Bp: numpy.ndarray[numpy.int32], Bj: numpy.ndarray[numpy.int32], Bx: numpy.ndarray[numpy.float32]) -> None
    2. (n_nodes: int, Ap: numpy.ndarray[numpy.int32], Aj: numpy.ndarray[numpy.int32], Ax: numpy.ndarray[numpy.float64], Sp: numpy.ndarray[numpy.int32], Sj: numpy.ndarray[numpy.int32], Sx: numpy.ndarray[numpy.float64], splitting: numpy.ndarray[numpy.int32], Bp: numpy.ndarray[numpy.int32], Bj: numpy.ndarray[numpy.int32], Bx: numpy.ndarray[numpy.float64]) -> None

This seems a software limitation, instead of a mathematical requirement? The RS-AMG interpolation formula should work for complex matrix naturally, based on the references:

• The paper Algebraic Multigrid Solvers for Complex-Valued Matrices, Section 3.2 "we use the natural complex extension of (2.5) to define interpolation"
• The paper A multigrid-based shifted Laplacian preconditioner for a fourth-order Helmholtz discretization, Section 3.2.2. AMG interpolation, only the strength-of-connection is slightly changed (use real part rather than modulus), while the interpolation formula looks exactly the same as real-valued case.

Are there any non-trivial modifications required, in order to apply RS-type AMG for complex matrices?

[learning-chip]

This rs_direct_interpolation_pass2 function is templated and should take complex arrays as well? So the TypeError comes from the pybind11 interface?

pyamg/pyamg/amg_core/ruge_stuben.h

Lines 721 to 733 in 1fd1ef6

	template<class I, class T>
	void rs_direct_interpolation_pass2(const I n_nodes,
	const I Ap[], const int Ap_size,
	const I Aj[], const int Aj_size,
	const T Ax[], const int Ax_size,
	const I Sp[], const int Sp_size,
	const I Sj[], const int Sj_size,
	const T Sx[], const int Sx_size,
	const I splitting[], const int splitting_size,
	const I Bp[], const int Bp_size,
	I Bj[], const int Bj_size,
	T Bx[], const int Bx_size)
	{

[learning-chip]

A breakdown:

import pyamg
A = pyamg.gallery.gauge_laplacian(100)
S = pyamg.strength.classical_strength_of_connection(A, norm='abs')  # works, although norm='min' requires real-number
splitting = pyamg.classical.classical.split.RS(S)  # works, since S is real
pyamg.classical.interpolate.direct_interpolation(A, S, splitting)  # TypeError at rs_direct_interpolation_pass2()

[learning-chip]

rs_direct_interpolation_pass2 contains comparisons to zero (Ax[jj] < 0) which is undefined for complex numbers. Does it make sense to use real part for comparison while keeping the rest of complex arithmetic unchanged?

pyamg/pyamg/amg_core/ruge_stuben.h

Lines 755 to 764 in ece7eb3

	} else {
	if (Ax[jj] < 0)
	sum_all_neg += Ax[jj];
	else
	sum_all_pos += Ax[jj];
	}
	}
	T alpha = sum_all_neg / sum_strong_neg;
	T beta = sum_all_pos / sum_strong_pos;

[bensworth]

Sorry for the delay. I see no problem generalizing to complex. Note that this is not RS interpolation, this is a simpler version, but the same principle applies, that the pointwise formula can be extended more or less directly to complex values. As Scott points out in his paper, the efficacy will depend on the problem. If you're interested in making a PR similar to the Julia package though please do. For strength I think we would want to use absolute value rather than comparison w/ real part, but this option already exists for real numbers so should not be hard to generalize to complex. @lukeolson thoughts?