diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc index 4560714d..e56c918a 100644 --- a/tests/sp2n/Test_project_on_Sp.cc +++ b/tests/sp2n/Test_project_on_Sp.cc @@ -107,62 +107,67 @@ void test_group_projections(T U) { template bool has_correct_algebra_block_structure(const T& U) { + // this only checks for the anti-hermitian part of the algebra const int nsp = Nc / 2; Complex i(0., 1.); std::cout << GridLogMessage << "Checking the structure is " << std::endl; - std::cout << GridLogMessage << "U = ( W X ) " << std::endl; - std::cout << GridLogMessage << " ( X^* -W^* ) " << std::endl; - std::cout << GridLogMessage << std::endl; - for (int c1 = 0; c1 < nsp; c1++) // check on W - { - for (int c2 = 0; c2 < nsp; c2++) { - auto W = PeekIndex(U, c1, c2); - auto Wstar = PeekIndex(U, c1 + nsp, c2 + nsp); - auto Ww = conjugate(Wstar); - auto amizero = sum(W + Ww); - auto amizeroo = TensorRemove(amizero); - assert(amizeroo.real() < 10e-6); - amizeroo *= i; - assert(amizeroo.real() < 10e-6); + std::cout << GridLogMessage << "U = ( W X ) " << std::endl; + std::cout << GridLogMessage << " ( -X^* W^* ) " << std::endl; + for (int c1 = 0; c1 < nsp; c1++) //check on W + { + for (int c2 = 0; c2 < nsp; c2++) + { + auto W = PeekIndex(U,c1,c2); + auto Wstar = PeekIndex(U,c1+nsp,c2+nsp); + auto Ww = conjugate( Wstar ); + auto amizero = sum(W - Ww); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } + } + for (int c1 = 0; c1 < nsp ; c1++) + { + for (int c2 = 0; c2 < nsp; c2++) + { + auto X = PeekIndex(U,c1,c2+nsp); + auto minusXstar = PeekIndex(U,c1+nsp,c2); + auto minusXx = conjugate(minusXstar); + auto amizero = sum (X + minusXx); + auto amizeroo = TensorRemove(amizero); + assert( amizeroo.real() < 10e-6 ); + amizeroo *= i; + assert( amizeroo.real() < 10e-6 ); + } } - } - for (int c1 = 0; c1 < nsp; c1++) { - for (int c2 = 0; c2 < nsp; c2++) { - auto X = PeekIndex(U, c1, c2 + nsp); - auto minusXstar = PeekIndex(U, c1 + nsp, c2); - auto minusXx = conjugate(minusXstar); - auto amizero = sum(X - minusXx); - auto amizeroo = TensorRemove(amizero); - assert(amizeroo.real() < 10e-6); - amizeroo *= i; - assert(amizeroo.real() < 10e-6); - } - } return true; } template -bool is_element_of_sp2n_algebra(T U) { - // does explicitly take a copy in order to not spoil the matrix for further - // use +bool is_element_of_sp2n_algebra(const T& U) { LatticeColourMatrixD aux(U.Grid()); LatticeColourMatrixD identity(U.Grid()); identity = 1.0; - - aux = U - adj(U); - std::cout << GridLogMessage << "T - Tda = " << norm2(aux) << std::endl; - aux = U + adj(U); - std::cout << GridLogMessage << "T + Tda = " << norm2(aux) << std::endl; - - std::cout << GridLogMessage << "Check that Omega U Omega = conj(U)" - << std::endl; - LatticeColourMatrixD Omega(U.Grid()); Sp::Omega(Omega); - aux = Omega * U * Omega - conjugate(U); - std::cout << GridLogMessage << "Omega U Omega - conj(U) = " << norm2(aux) + + std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl; + assert(norm2(U) > 1e-8); + + aux = U - adj(U); + std::cout << GridLogMessage << "T - Tda = " << norm2(aux) + << " (not supposed to vanish)" << std::endl; + + aux = U + adj(U); + std::cout << GridLogMessage << "T + Tda = " << norm2(aux) + << " (supposed to vanish)" << std::endl; + assert(norm2(aux) - 1 < 1e-8); + + std::cout << GridLogMessage << "Check that Omega T Omega + conj(T) = 0 " << std::endl; + aux = Omega * U * Omega + conjugate(U); assert(norm2(aux) < 1e-8); return has_correct_algebra_block_structure(U);