Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Tests/FermionMatrix/test_FermionMatrix_hermitian.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ COMPLEX_NSL_TEST_CASE( "FermionMatrix: MdaggerM hermitian", "[FermionMatrix, Mda
//Test if MMdagger is hermitian
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_FermionMatrix_MMdagger_hermitian(const NSL::size_t nt, LatticeType & Lattice, const std::string latticeName, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
typedef NSL::complex<NSL::RealTypeOf<Type>> ComplexType;
NSL::size_t nx = Lattice.sites();
NSL::Tensor<Type> phi(nt, nx);
phi.rand();
Expand Down Expand Up @@ -90,7 +90,7 @@ void test_FermionMatrix_MMdagger_hermitian(const NSL::size_t nt, LatticeType & L
//Test if MdaggerM is hermitian
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_FermionMatrix_MdaggerM_hermitian(const NSL::size_t nt, LatticeType & Lattice, const std::string latticeName, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
typedef NSL::complex<NSL::RealTypeOf<Type>> ComplexType;
NSL::size_t nx = Lattice.sites();
NSL::Tensor<Type> phi(nt, nx);
phi.rand();
Expand Down
22 changes: 7 additions & 15 deletions Tests/FermionMatrix/test_fermionMatrixHubbardDiag.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,6 @@ COMPLEX_NSL_TEST_CASE( "fermionMatrixHubbardDiag: logDetM_noninteracting", "[fer

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardDiag_M_dense(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Expand Down Expand Up @@ -164,7 +163,6 @@ void test_fermionMatrixHubbardDiag_M_dense(const NSL::size_t nt, LatticeType & L
//Test for the function Mdagger(psi)
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardDiag_Mdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Expand Down Expand Up @@ -206,7 +204,6 @@ void test_fermionMatrixHubbardDiag_Mdagger(const NSL::size_t nt, LatticeType & L
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardDiag_MMdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Expand All @@ -216,7 +213,7 @@ void test_fermionMatrixHubbardDiag_MMdagger(const NSL::size_t nt, LatticeType &

NSL::FermionMatrix::HubbardDiag M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.MMdagger(psi);
auto indirect = M.M(M.Mdagger(psi));
Expand All @@ -233,7 +230,6 @@ void test_fermionMatrixHubbardDiag_MMdagger(const NSL::size_t nt, LatticeType &
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardDiag_MdaggerM(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

//hardcoding the calculation done in the method MdaggerM of fermionMatrixHubbardDiag class
Expand All @@ -244,7 +240,7 @@ void test_fermionMatrixHubbardDiag_MdaggerM(const NSL::size_t nt, LatticeType &

NSL::FermionMatrix::HubbardDiag M(Lattice,nt,beta);
M.populate(phi,NSL::Hubbard::Particle);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.MdaggerM(psi);
auto indirect = M.Mdagger(M.M(psi));
Expand All @@ -265,7 +261,6 @@ void test_logDetM_time_shift_invariance(const NSL::size_t nt, LatticeType & Latt
// We should find that shifting phi in time doesn't change the determinant.
int slices_to_shift_by=4;

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
size_t nx = Lattice.sites();
NSL::Tensor<Type> phi(nt, nx), phiShift(nt, nx);
phi.rand();
Expand Down Expand Up @@ -301,8 +296,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c

// We should find that by shifting any element of phi by 2π the real part of the determinant doesn't change.

typedef typename NSL::RT_extractor<Type>::value_type RealType;
typedef NSL::complex<RealType> ComplexType;
NSL::size_t nx = Lattice.sites();
Type delta = beta/nt;

Expand All @@ -321,7 +314,7 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
random.rand();
NSL::Tensor<Type> orbits = static_cast<NSL::Tensor<int>>(10*random);

RealType two_pi = 2*std::numbers::pi_v<RealType>;
NSL::RealTypeOf<Type> two_pi = 2*std::numbers::pi_v<NSL::RealTypeOf<Type>>;
phiShift = phi + two_pi * orbits;

NSL::FermionMatrix::HubbardDiag<Type,LatticeType> M (Lattice,nt ,beta);
Expand All @@ -331,9 +324,9 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c

Type result = M.logDetM();
Type result_shift = Mshift.logDetM();
RealType diff_imag_mod_two_pi = std::remainder(
static_cast<RealType>(NSL::imag(result - result_shift)),
static_cast<RealType>(two_pi)
NSL::RealTypeOf<Type> diff_imag_mod_two_pi = std::remainder(
static_cast<NSL::RealTypeOf<Type>>(NSL::imag(result - result_shift)),
static_cast<NSL::RealTypeOf<Type>>(two_pi)
);


Expand All @@ -345,7 +338,7 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c

//comparing only the real parts
REQUIRE(almost_equal(result_shift.real(),result.real(),std::numeric_limits<Type>::digits10-1));
REQUIRE(almost_equal(static_cast<RealType>(0),
REQUIRE(almost_equal(static_cast<NSL::RealTypeOf<Type>>(0),
diff_imag_mod_two_pi,
std::numeric_limits<Type>::digits10-3));

Expand All @@ -360,7 +353,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_logDetM_noninteracting(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx), sausage = NSL::Matrix::Identity<Type>(nx);
Expand Down
57 changes: 22 additions & 35 deletions Tests/FermionMatrix/test_fermionMatrixHubbardExp.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,6 @@ COMPLEX_NSL_TEST_CASE( "fermionMatrixHubbardExp: logDetM_uniform_timeslices", "[
//Test for the function M(psi)
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_M(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();
//hardcoding the calculation done in the method M of fermionMatrixHubbardExp class
NSL::Tensor<Type> phi(nt, nx);
Expand All @@ -153,7 +152,7 @@ void test_fermionMatrixHubbardExp_M(const NSL::size_t nt, LatticeType & Lattice,
Type delta = beta/nt;
NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I ={0,1};
Type I ={0,1};

// apply kronecker delta
//NSL::Tensor<Type> psiShift = NSL::LinAlg::shift(psi,1);
Expand Down Expand Up @@ -183,7 +182,6 @@ void test_fermionMatrixHubbardExp_M(const NSL::size_t nt, LatticeType & Lattice,

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_M_dense(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Expand Down Expand Up @@ -219,42 +217,41 @@ void test_fermionMatrixHubbardExp_M_dense(const NSL::size_t nt, LatticeType & La
//Test for the function Mdagger(psi)
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();
INFO("nt: "+NSL::to_string(nt)+" nx: "+NSL::to_string(nx));

NSL::Tensor<ComplexType> phi(nt, nx);
NSL::Tensor<ComplexType> psi(nt, nx);
NSL::Tensor<Type> phi(nt, nx);
NSL::Tensor<Type> psi(nt, nx);
phi.rand();
psi.rand();

// To simplify this test one can force the field
//phi = phi.real() + ComplexType(0,0); // phi to be real
//psi = psi.real() + ComplexType(0,0); // psi to be real
//psi = psi.imag() * ComplexType(0,1); // psi to be imaginary
//phi = phi.real() + Type(0,0); // phi to be real
//psi = psi.real() + Type(0,0); // psi to be real
//psi = psi.imag() * Type(0,1); // psi to be imaginary
NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);


// First let's check Mdagger against the dagger of the dense implementation of M.

// Construct a dense representation of M† from a dense representation of M
NSL::Tensor<ComplexType> Mdense_dagger = M.M_dense(nt).transpose(0,2).transpose(1,3).conj();
NSL::Tensor<Type> Mdense_dagger = M.M_dense(nt).transpose(0,2).transpose(1,3).conj();

// We can also apply Mdagger to the identity matrix in order to get a dense Mdagger.
// Follow the M_dense implementation:

NSL::Tensor<ComplexType> dense(nt, nx, nt, nx);
NSL::Tensor<Type> dense(nt, nx, nt, nx);

// Construct the identity matrix.
NSL::Tensor<ComplexType> identity(nt, nx, nt, nx);
NSL::Tensor<Type> identity(nt, nx, nt, nx);
for(int t = 0; t < nt; t++){
identity(t,NSL::Slice(), t, NSL::Slice()) = NSL::Matrix::Identity<ComplexType>(nx);
identity(t,NSL::Slice(), t, NSL::Slice()) = NSL::Matrix::Identity<Type>(nx);
}

// Ensure it's really the identity in the mat-vec sense.
// Apply the identity to psi via obvious mat-vec
NSL::Tensor<ComplexType> Ipsi(nt, nx);
NSL::Tensor<Type> Ipsi(nt, nx);
for(int t=0; t < nt; t++){
for(int x=0; x < nx; x++){
for(int i=0; i< nt; i++){
Expand All @@ -275,7 +272,7 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La
}

// So, we can compare M†.I to (M.I)†
REQUIRE( almost_equal(Mdense_dagger-dense, ComplexType(0,0)).all() );
REQUIRE( almost_equal(Mdense_dagger-dense, Type(0,0)).all() );
// This REQUIREment looks funny; why not just check that the two tensors are almost_equal directly, as in
// REQUIRE( almost_equal(Mdense_dagger, dense).all() );
// TODO: almost_equal of +0.0000... and -0.0000... evaluates to False and that's extremely misleading
Expand All @@ -284,7 +281,7 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La


// Finally, compare two ways of computing M†ψ
NSL::Tensor<ComplexType>Mdense_dagger_psi(nt, nx), M_dagger_psi(nt, nx);
NSL::Tensor<Type>Mdense_dagger_psi(nt, nx), M_dagger_psi(nt, nx);
// by doing the obvious mat-vec,
for(int t=0; t < nt; t++){
for(int x=0; x < nx; x++){
Expand All @@ -308,7 +305,6 @@ void test_fermionMatrixHubbardExp_Mdagger(const NSL::size_t nt, LatticeType & La
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_MdaggerM(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

//hardcoding the calculation done in the method MdaggerM of fermionMatrixHubbardExp class
Expand All @@ -319,7 +315,7 @@ void test_fermionMatrixHubbardExp_MdaggerM(const NSL::size_t nt, LatticeType & L

NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.MdaggerM(psi);
auto indirect = M.Mdagger(M.M(psi));
Expand All @@ -335,7 +331,6 @@ void test_fermionMatrixHubbardExp_MdaggerM(const NSL::size_t nt, LatticeType & L

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_MMdagger(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Expand All @@ -345,7 +340,7 @@ void test_fermionMatrixHubbardExp_MMdagger(const NSL::size_t nt, LatticeType & L

NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.MMdagger(psi);
auto indirect = M.M(M.Mdagger(psi));
Expand All @@ -361,7 +356,6 @@ void test_fermionMatrixHubbardExp_MMdagger(const NSL::size_t nt, LatticeType & L

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_M_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::size_t Nbatch = 10;
Expand All @@ -373,7 +367,7 @@ void test_fermionMatrixHubbardExp_M_batched(const NSL::size_t nt, LatticeType &

NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.M(psi);

Expand All @@ -389,7 +383,6 @@ void test_fermionMatrixHubbardExp_M_batched(const NSL::size_t nt, LatticeType &

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_Mdagger_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::size_t Nbatch = 10;
Expand All @@ -401,7 +394,7 @@ void test_fermionMatrixHubbardExp_Mdagger_batched(const NSL::size_t nt, LatticeT

NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.Mdagger(psi);

Expand All @@ -417,7 +410,6 @@ void test_fermionMatrixHubbardExp_Mdagger_batched(const NSL::size_t nt, LatticeT

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_MMdagger_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::size_t Nbatch = 10;
Expand All @@ -429,7 +421,7 @@ void test_fermionMatrixHubbardExp_MMdagger_batched(const NSL::size_t nt, Lattice

NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.MMdagger(psi);

Expand All @@ -445,7 +437,6 @@ void test_fermionMatrixHubbardExp_MMdagger_batched(const NSL::size_t nt, Lattice

template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_fermionMatrixHubbardExp_MdaggerM_batched(const NSL::size_t nt, LatticeType & Lattice, const Type & beta){
typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::size_t Nbatch = 10;
Expand All @@ -457,7 +448,7 @@ void test_fermionMatrixHubbardExp_MdaggerM_batched(const NSL::size_t nt, Lattice

NSL::FermionMatrix::HubbardExp M(Lattice,nt,beta);
M.populate(phi);
ComplexType I={0,1};
Type I={0,1};

auto direct = M.MdaggerM(psi);

Expand All @@ -477,7 +468,6 @@ void test_logDetM_time_shift_invariance(const NSL::size_t nt, LatticeType & Latt
// We should find that shifting phi in time doesn't change the determinant.
int slices_to_shift_by=4;

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
size_t nx = Lattice.sites();
NSL::Tensor<Type> phi(nt, nx), phiShift(nt, nx);
phi.rand();
Expand Down Expand Up @@ -512,7 +502,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c

// We should find that by shifting any element of phi by 2π the determinant doesn't change.

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();
Type delta = beta/nt;

Expand Down Expand Up @@ -552,7 +541,6 @@ void test_logDetM_phi_plus_two_pi(const NSL::size_t nt, LatticeType & Lattice, c
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_logDetM_noninteracting(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Expand Down Expand Up @@ -583,12 +571,11 @@ void test_logDetM_noninteracting(const NSL::size_t nt, LatticeType & Lattice, co
template<NSL::Concept::isNumber Type, NSL::Concept::isDerived<NSL::Lattice::SpatialLattice<Type>> LatticeType>
void test_logDetM_uniform_timeslices(const NSL::size_t nt, LatticeType & Lattice, const Type & beta) {

typedef NSL::complex<typename NSL::RT_extractor<Type>::value_type> ComplexType;
NSL::size_t nx = Lattice.sites();

NSL::Tensor<Type> phi(nt, nx);
Type delta = beta/nt;
ComplexType I ={0,1};
Type I ={0,1};

// When phi on a given timeslice is the same on every spatial site
NSL::Tensor<Type> tmp (nt); tmp.rand();
Expand All @@ -598,9 +585,9 @@ void test_logDetM_uniform_timeslices(const NSL::size_t nt, LatticeType & Lattice
// exp(i phi(t)) matrix is proportional to the identity matrix and can be
// treated like a scalar.
// When EVERY timeslice is like that we gather all the scalars together
NSL::Tensor<ComplexType> sum(1);
NSL::Tensor<Type> sum(1);
sum(0) = I*phi( NSL::Slice(), 0).sum();
ComplexType expsum = NSL::LinAlg::exp(sum)(0);
Type expsum = NSL::LinAlg::exp(sum)(0);

// to get
// logdet M = logdet( 1 + exp(sum(phi(t))) exp(kappa_tilde * Nt))
Expand Down
Loading