Use Hermitian transpose for complex-valued constraint

[conpierce8]

When using complex values, it is advantageous to pre-multiply the system by the Hermitian transpose of the prolongation matrix rather than the standard transpose, as this preserves the Hermitian character of the original matrix, i.e. applying the prolongation matrix as $$K^H A K = K^H b$$ yields a Hermitian system matrix $K^H A K$ if $A$ is Hermitian, whereas applying the prolongation matrix as $$K^T A K = K^T b$$ yields a non-Hermitian matrix if $K$ has complex values. (Note $K^H = \overline{K^T}$, where $\overline{K}$ is the complex conjugate of $K$.)

This should probably be the default behavior for complex numbers, since destroying the Hermitian character of $A$ is not likely to be advantageous to users. This behavior could, however, be enabled/disabled by the user, e.g. by adding an option to MultiPointConstraint (which would be ignored for real-valued builds):

class MultiPointConstraint():
    ...
    def __init__(self, V: _fem.FunctionSpace, use_hermitian_transpose: bool = True):
        ...

[conpierce8]

@jorgensd I took a crack at implementing this: https://github.com/conpierce8/dolfinx_mpc/tree/hermitian. I've only implemented it for assemble_matrix so far; maybe I'll take a look at assemble_vector soon. I've lightly tested my implementation with some small meshes and so far it matches the matrix I obtain using a manually-populated prolongation matrix.

[bleyerj]

Hi, I confirm that I encountered a similar issue for my application. Preserving Hermitian character of the reduced matrix is definitely needed in practice. I will test your implementation soon and let you know.

[jorgensd]

I've got no issues with changing the default behavior to be the Hermitian, as it makes more sense for complex numbers, and doesn't change the behavior in real mode. Shouldn't there be a conj here: https://github.com/jorgensd/dolfinx_mpc/compare/main...conpierce8:dolfinx_mpc:hermitian?expand=1#diff-1a2d251f84ec0019705af3746d92fae95ff774bee9dd824aa859df5ce7923396R173 ?

[conpierce8]

I will test your implementation soon and let you know.

@bleyerj: Great, thanks!

Shouldn't there be a conj here: https://github.com/jorgensd/dolfinx_mpc/compare/main...conpierce8:dolfinx_mpc:hermitian?expand=1#diff-1a2d251f84ec0019705af3746d92fae95ff774bee9dd824aa859df5ce7923396R173 ?

@jorgensd: I think it's okay as-is...? I wrote two implementations apiece for modify_mpc_cell_rows and modify_mpc_cell_masters, one for double and one for std::complex<double>. You've linked to the double version, where conj would have no effect. (conj is present on line 205 in the complex version.) I separated the two implementations because I thought it would improve performance in the real case if I avoided calling conj during assembly. Am I correct in my thinking, or would the change in performance be minimal?

[jorgensd]

I guess it would be better to use templating and a if constexpr clause to avoid code duplication. See for instance: https://github.com/FEniCS/dolfinx/blob/b5dc1daa76958bb5d5110a1b70bbdf862c9f31e0/cpp/dolfinx/io/xdmf_function.cpp#L109 and https://github.com/FEniCS/dolfinx/blob/627059a4dcc694bed249b3a6f0bddd837357fa45/cpp/dolfinx/fem/assemble_vector_impl.h#L138-L170

[conpierce8]

I guess it would be better to use templating and a if constexpr clause to avoid code duplication.

:flushed: Definitely. I didn't know about constexprs. C++ rookie mistake, I guess. This will simplify the implementation greatly; I no longer see any need to create separate methods for modify_mpc_cell_rows and modify_mpc_cell_masters. I'll just add if constexpr at the appropriate places, namely here and here.

This should probably be the default behavior for complex numbers, since destroying the Hermitian character of A is not likely to be advantageous to users. This behavior could, however, be enabled/disabled by the user, e.g. by adding an option to MultiPointConstraint (which would be ignored for real-valued builds):

class MultiPointConstraint():
    ...
    def __init__(self, V: _fem.FunctionSpace, use_hermitian_transpose: bool = True):
        ...

@jorgensd @bleyerj: Do either of you see a need to implement such an option? Would it be better to place that option in MultiPointConstraint (as I suggested above) or in assemble_matrix? E.g.:

def assemble_matrix(form: _fem.FormMetaClass,
                    constraint: Union[MultiPointConstraint,
                                      Sequence[MultiPointConstraint]],
                    bcs: Sequence[_fem.DirichletBCMetaClass] = [],
                    diagval: _PETSc.ScalarType = 1,
                    A: Optional[_PETSc.Mat] = None,
                    premultiply_hermitian: bool = True) -> _PETSc.Mat:
    ...

[jorgensd]

I think we could force the code to always use the hermitian, and not leave the option for $K^TAK$ multiplication. I do not see any use-cases where it would be favorable over $K^HAK$. If either of you have any comments on this please let me know.

If we choose to keep the option, it should probably be on the assemble_matrix/vector level, even if this introduces the possibility of people creating inconsistent systems, i.e. $K^HAK = K^TB$.

[bleyerj]

I agree, I don't see any use-case where the transpose would be preferred either. Let's keep this simple.

[conpierce8]

Sounds logical. If someone has a use case, they can make a feature request. constexpr implementation is pushed; pull request opened.