Use symmetric block matrix for complex coarse solve

[zhanggy2004]

Thanks for developing and maintaining this nice project.

For "complex coarse solve", instead of solving

$$A_1\begin{pmatrix}x_r \\ x_i \end{pmatrix}:=\begin{pmatrix} A_r & -A_i \\ A_i & A_r \end{pmatrix} \begin{pmatrix}x_r \\ x_i \end{pmatrix}=\begin{pmatrix}y_r \\ y_i \end{pmatrix},$$

we can solve

$$A_2\begin{pmatrix}x_r \\ -x_i \end{pmatrix}:=\begin{pmatrix} A_r & A_i \\ A_i & -A_r \end{pmatrix} \begin{pmatrix}x_r \\ -x_i \end{pmatrix}=\begin{pmatrix}y_r \\ y_i \end{pmatrix}.$$

The idea is similar to mfem::ComplexOperator::Convention::BLOCK_SYMMETRIC.

When $A_r$ and $A_i$ are both symmetric, $A_2$ will also be symmetric, which can be utitlized by MUMPS solver for more efficient solve.

I realized that $A_r$ and $A_i$ are unsymmetric when there is periodic boundary (?), but am not sure if there are other exceptions. Please let me know your suggestions for the correct/better approach. Thanks.

[simlapointe]

This looks good! I tested it on a few different cases with Floquet periodic BCs (cylinder, nail) and without (cpw, single transmon, adapter waveguide) and it works as expected in all cases. I checked that mumps, strumpack, and superlu all converge in 1-2 iterations when ComplexCoarseSolve == true.

I think Floquet periodicity is the only symmetry-breaking term we have at the moment, so the wave_vector == std::array<double, 3>{0.0, 0.0, 0.0} check is ok (though I don't love == 0.0). However, if we ever encounter other symmetry-breaking scenarios, it might be better to define a single variable such as iodata.solver.linear.is_symmetric that gets set accordingly after reading the config. @hughcars, what do you think?

[hughcars]

This looks good! I tested it on a few different cases with Floquet periodic BCs (cylinder, nail) and without (cpw, single transmon, adapter waveguide) and it works as expected in all cases. I checked that mumps, strumpack, and superlu all converge in 1-2 iterations when ComplexCoarseSolve == true.

I think Floquet periodicity is the only symmetry-breaking term we have at the moment, so the wave_vector == std::array<double, 3>{0.0, 0.0, 0.0} check is ok (though I don't love == 0.0). However, if we ever encounter other symmetry-breaking scenarios, it might be better to define a single variable such as iodata.solver.linear.is_symmetric that gets set accordingly after reading the config. @hughcars, what do you think?

I think it's ok for now, but ultimately we want this kind of information extract as part of a validator that ingests a configuration object. Are there similar symmetry flags for strumpack and superlu we might exploit?

[simlapointe]

I think it's ok for now, but ultimately we want this kind of information extract as part of a validator that ingests a configuration object. Are there similar symmetry flags for strumpack and superlu we might exploit?

Not that I'm aware of. SuperLU has SetSymmetricPattern and STRUMPACKRowLocMatrix has a sym_sparse argument which are about the sparsity pattern and we already always set them to true, but nothing about the matrix data being symmetric. I think only mumps cares about that.