23 template <
typename Matrix,
typename Float>
29 for (
int i=0; i<in.
size(); i++) {
31 for (
int j=0; j<in.
size(); j++) {
32 if (fabs(
in(i,j).real() - inv(j,i).real()) > tol ||
33 fabs(
in(i,j).imag() + inv(j,i).imag()) >
tol)
return false;
41 for (
int i=0; i<in.
size(); i++) {
42 if (fabs(identity(i,i).real() - static_cast<Float>(1.0)) > tol ||
43 fabs(identity(i,i).imag()) >
tol)
46 for (
int j=0; j<in.
size(); j++) {
48 if (fabs(identity(i,j).real()) > tol || fabs(identity(i,j).imag()) > tol ||
49 fabs(identity(j,i).real()) > tol || fabs(identity(j,i).imag()) >
tol )
65 template <
typename Matrix>
68 for (
int i=0; i<in.
size(); i++) {
69 for (
int j=0; j<in.
size(); j++) {
70 printf(
"TESTR: %+.13le %+.13le %+.13le\n",
71 in(i,j).real(), inv(j,i).real(), fabs(
in(i,j).real() - inv(j,i).real()));
72 printf(
"TESTI: %+.13le %+.13le %+.13le\n",
73 in(i,j).imag(), inv(j,i).imag(), fabs(
in(i,j).imag() + inv(j,i).imag()));
86 template <
typename Float>
89 constexpr Float negative_third = -1.0/3.0;
90 constexpr Float negative_sixth = -1.0/6.0;
94 constexpr
int max_iter = 100;
97 out = 0.5*(out +
conj(inv));
99 }
while (!
checkUnitary(inv, out, tol) && ++i < max_iter);
104 Float angle =
arg(det);
106 complex<Float> cTemp;
107 sincos(negative_third * angle, &cTemp.y, &cTemp.x);
109 in = (mod*cTemp)*out;
__host__ __device__ void polarSu3(Matrix< complex< Float >, 3 > &in, Float tol)
Project the input matrix on the SU(3) group. First unitarize the matrix and then project onto the spe...
__device__ __host__ constexpr int size() const
__host__ __device__ ValueType norm(const complex< ValueType > &z)
Returns the magnitude of z squared.
__host__ __device__ bool checkUnitary(const Matrix &inv, const Matrix &in, const Float tol)
Check the unitarity of the input matrix to a given tolerance.
__host__ __device__ void checkUnitaryPrint(const Matrix &inv, const Matrix &in)
Print out deviation for each component (used for debugging only).
__host__ __device__ ValueType pow(ValueType x, ExponentType e)
__device__ __host__ Matrix< T, 3 > inverse(const Matrix< T, 3 > &u)
cpuColorSpinorField * out
__host__ __device__ ValueType arg(const complex< ValueType > &z)
Returns the phase angle of z.
static int mod(int a, int b)
__device__ __host__ T getDeterminant(const Mat< T, 3 > &a)
__host__ __device__ ValueType conj(ValueType x)