Feature/overlap

[SaltyChiang]

Implemented the overlap chiral fermion in QUDA.

Because the overlap fermion breaks the locality, we cannot build an even-odd preconditioned Dirac matrix, so only the unpreconditioned version is implemented.

The DiracOverlap::Dslash is defined as below:

$$ D\equiv D_\mathrm{overlap}=\frac{1}{2}\left[1+\gamma_5\mathrm{sign}\left(\gamma_5M_\mathrm{Wilson}\right)\right] $$

Take the chiral projector $P^\pm=\frac{1\pm\gamma_5}{2}$. Under the DeGrand-Rossi basis, we get the upper two spinor rows with $P^+$, and the lower two spinor rows with $P^-$ if we apply the projector to a spinor.

We can easily get:

$$ D^\dagger D=DD^\dagger=\frac{1}{2}\left(D^\dagger+D\right),\text{ and then } D^\dagger DP^\pm=P^\pm DP^\pm $$

$D^\dagger DP^\pm=P^\pm DP^\pm$ allows us to separate vectors into two parts $v=P^+v+P^-v=v^++v^-$, and $D^\dagger D$ will return a pure upper(lower) vector if the input vector is pure upper(lower). Thus, all solvers can work on a single chirality instead of a full spinor. This is why we added the new MdagMChiral interface in Dirac, along with the new DiracMdagMChiral class. We should always use a MdagMChiral type solver when dealing with overlap fermions.

Then we define the chiral fermion matrix

$$ M_\mathrm{chiral}=m+\frac{D}{1-D} $$

We define DiracOverlap::M and DiracOverlap::MdagM as

$$ M=m+\left(1-m\right)D,\quad M^\dagger M=m^2+\left(1-m^2\right)D^\dagger D $$

Though it's a bit weird, it will simplify the expression of the chiral fermion matrix inverse

$$ \begin{aligned} M_\mathrm{chiral}^{-1}=\left(1-D\right)M^{-1}=\frac{1}{1-m}\left[M^{-1}-1\right] \\ \left(M_\mathrm{chiral}^\dagger M_\mathrm{chiral}\right)^{-1}=\left(1-D^\dagger D\right)\left(M^\dagger M\right)^{-1}=\frac{1}{1-m^2}\left[\left(M^\dagger M\right)^{-1}-1\right] \end{aligned} $$

These steps are implemented in the reconstruct function.

However, some issues arise when we are dealing with the multishift solver. We have to set the mass of a DiracOverlap to 0 and set offsets to $\frac{m^2}{1-m^2}$, and apply $\frac{1}{1-m^2}$ as the postprocessing. We also need to call construct in the multishift solving process.

[maddyscientist]

Thanks for this amazing contribution @SaltyChiang. I've enabled CI for the PR, but see quite a failures, we can review this properly once you've resolved this.

[maddyscientist]

What do these changes to the Eigen solver do?

[SaltyChiang]

I'm sorry, this code should be reverted. When I want to match QUDA's eigensolver with the textbook, I noticed that there are some differences. It shows that QUDA's eigensolver takes lambda1 = 0.0 and the following sigma1 and d2 can be simplified. It shouldn't matter which lambda1 is selected.

[weinbe2]

This is an awesome contribution, @SaltyChiang ! I'm looking forward to looking into this more closely.

One interface comment specific to how overlap only admits a full-parity form. The Kahler-Dirac preconditioned staggered operator also only admits a full parity form, and the convention I adopted for the Dirac* class was to leave Dslash and DslashXpay undefined (error-throwing): https://github.com/lattice/quda/blob/develop/lib/dirac_improved_staggered_kd.cpp#L37 , and only leave M, MdagM, etc, defined.

The reason I settled on this convention was because in all other cases Dslash + DslashXpay only act on single-parity fields (and a lot of code that calls those routines assumes as such), and I didn't want to deviate from that. (Whether or not that's a good convention is another story.) I didn't run into any fundamental issues from doing this, and I think in some cases it simplified things.

Unless you have a very specific reason to define Dslash and DslashXpay along with M, etc, I'd recommend matching what I did for the KD operators.

[SaltyChiang]

@maddyscientist I hope compilers will be happy with the update. Also, I reverted some code in eigensolve, now those parameters are the same as before.
@weinbe2 Thank you for the explanation. I disabled Dslash and DslashXpay in DiracOverlap.

[havogt]

cscs-ci run

[SaltyChiang]

@havogt Thank you for triggering the CI. I mistakenly changed a parameter and now it's reverted.

[maddyscientist]

cscs-ci run

[SaltyChiang]

@maddyscientist Those failed tests are caused by the modification in the multishift solver. I reverted the offsets-related part and all tests passed on my local machine.