14bool report_pass(std::string message,
double eps,
double limit) {
16 static int is_terminal = -1;
17 static std::string ok_string, fail_string;
19 if (is_terminal < 0) {
20 is_terminal = isatty(fileno(stdout));
23 ok_string =
"\x1B[32m --- \033[0m ";
24 fail_string =
"\x1B[31m *** \033[0m ";
27 fail_string =
" *** ";
32 hila::out0 << ok_string << message <<
" passed" << std::endl;
35 hila::out0 << fail_string << message <<
" FAILED: eps " << eps <<
" limit " << limit
43void check_reductions() {
49 f[
ALL] = expi(2 * M_PI * X.coordinate(e_x) / lattice.size(e_x));
52 onsites(
ALL) sum += f[X];
54 sum /= lattice.volume();
55 report_pass(
"Complex reduction value " + hila::prettyprint(sum),
abs(sum), 1e-4);
60 rv[X.coordinate(e_x)] += f[X];
64 for (
int i = 0; i < lattice.size(e_x); i++) {
66 expi(2 * M_PI * i / lattice.size(e_x)) - rv[i] / (lattice.volume() / lattice.size(e_x));
68 report_pass(
"Vector reduction, sum " + hila::prettyprint(sum),
abs(sum), 1e-4);
75void test_site_access() {
81 for (
int i = 0; i < 3; i++) {
87 report_pass(
"Setting and reading a value at " + hila::prettyprint(c.
transpose()),
106 report_pass(
"Maxloc is " + hila::prettyprint(loc.
transpose()), (c - loc).norm(), 1e-8);
107 report_pass(
"Max value " + hila::prettyprint(v), v - 2, 1e-9);
114 report_pass(
"Minloc is " + hila::prettyprint(loc.
transpose()), (c - loc).norm(), 1e-8);
115 report_pass(
"Min value " + hila::prettyprint(v), v + 1, 1e-9);
124 constexpr int n_loops = 100;
128 double fsum = 0, fsqr = 0;
129 for (
int i=0; i<n_loops; i++) {
131 double s = 0, s2 = 0;
137 fsum += s/lattice.volume();
138 fsqr += s2/lattice.volume();
144 report_pass(
"Gaussian random average (6 sigma limit) " + hila::prettyprint(fsum),
145 abs(fsum), 6/sqrt(((
double)n_loops)*lattice.volume()));
147 report_pass(
"Gaussian random width^2 " + hila::prettyprint(fsqr),
148 fsqr - 1, 6/sqrt(((
double)n_loops)*lattice.volume()));
157void test_set_elements_and_select() {
162 std::vector<CoordinateVector> cvec;
163 std::vector<Complex<double>> vals;
168 for (
int i = 0; i <= 50; i++) {
187 f.set_elements(vals, cvec);
192 for (
int i = 0; i < vals.size(); i++) {
196 report_pass(
"Field set_elements and get_elements with " + hila::prettyprint(vals.size()) +
203 SiteValueSelect<Complex<double>> sv;
212 report_pass(
"SiteSelect size " + hila::prettyprint(s.size()), s.size() - cvec.size(), 1e-3);
213 report_pass(
"SiteValueSelect size " + hila::prettyprint(sv.size()), sv.size() - cvec.size(),
218 for (
auto &c : cvec) {
220 for (
int i = 0; !found && i < s.size(); i++) {
221 if (s.coordinates(i) == c)
228 report_pass(
"SiteSelect content", 1 - ok, 1e-6);
231 for (
int k = 0; k < cvec.size(); k++) {
234 for (
int i = 0; !found && i < sv.size(); i++) {
235 if (sv.coordinates(i) == cvec[k] && sv.value(i) == vals[k])
242 report_pass(
"SiteValueSelect content", 1 - ok, 1e-6);
250void test_subvolumes() {
253 for (
int i = 0; i < 20; i++) {
257 if (si.coordinates() != c)
261 report_pass(
"SiteIndex", ok ==
false, 1e-2);
269 size_t vol = lattice.volume();
276 vol /= lattice.size(d);
279 report_pass(hila::prettyprint(NDIM - (
int)d - 1) +
"-dimensional slice size " +
280 hila::prettyprint(vol),
281 slice.size() - vol, 1e-3);
298 for (
auto s : slice) {
303 if (add && c[d2] < 0) {
305 if (mycoord[d2] < lattice.size(d2)) {
313 if (pass && mycoord != s.coordinates()) {
315 << hila::prettyprint(mycoord.
transpose()) <<
" is "
316 << hila::prettyprint(s.coordinates().transpose()) <<
" slice is "
317 << hila::prettyprint(c.
transpose()) <<
'\n';
323 report_pass(
"slice content", pass ==
false, 1e-2);
344 p2[{0, 0, 0}] = lattice.volume();
348 double eps = squarenorm_relative(p, p2);
350 report_pass(
"FFT constant field", eps, 1e-13 * sqrt(lattice.volume()));
360 sum += (f[X] - lattice.volume()).
squarenorm();
364 eps = fabs(sum / tnorm);
366 report_pass(
"FFT inverse transform", eps, 1e-10);
369 for (
int iter = 0; iter < 5; iter++) {
379 double d = kv.
dot(X.coordinates());
386 p2[kx] = lattice.volume();
388 eps = squarenorm_relative(p, p2);
390 report_pass(
"FFT of wave vector " + hila::prettyprint(kx.
transpose()), eps,
391 1e-13 * sqrt(lattice.volume()));
401 p = f.FFT(fft_direction::back) / lattice.volume();
402 eps = squarenorm_relative(r, p);
404 report_pass(
"FFT real to complex", eps, 1e-13 * sqrt(lattice.volume()));
406 auto r2 = f.FFT_complex_to_real(fft_direction::back) / lattice.volume();
407 eps = squarenorm_relative(r, r2);
409 report_pass(
"FFT complex to real", eps, 1e-13 * sqrt(lattice.volume()));
418 f = p.FFT(fft_direction::back) / sqrt(lattice.volume());
422 report_pass(
"Norm of field = " + hila::prettyprint(nf) +
" and FFT = " + hila::prettyprint(np),
423 (nf - np) / nf, 1e-10);
426 b.
k_max(M_PI * sqrt(3.0));
436 report_pass(
"Norm of binned FFT = " + hila::prettyprint(s), (s - np) / np, 1e-10);
441void spectraldensity_test() {
449 b.
k_max(M_PI * sqrt(3.0));
453 for (
int iter = 0; iter < 3; iter++) {
460 auto absk = kv.
norm();
469 for (
int i = 0; i < b.
bins(); i++) {
470 if (b.
bin_min(i) <= absk && b.bin_max(i) > absk) {
478 report_pass(
"Binning test at vector " + hila::prettyprint(kx.
transpose()), sum, 1e-10);
483 double d = kv.
dot(X.coordinates());
490 for (
int i = 0; i < b.
bins(); i++) {
491 if (b.
bin_min(i) <= absk && b.bin_max(i) > absk) {
493 sum += fabs(fsd[i] / pow(lattice.volume(), 2) - 1);
499 report_pass(
"Spectral density test with above vector ", sum, 1e-10);
505void test_field_slices() {
516void test_matrix_algebra() {
523 M.gaussian_random(2.0);
529 auto H = M[X] * M[X].dagger();
531 auto r = H.eigen_hermitean();
532 delta[X] = (H - r.eigenvectors * r.eigenvalues * r.eigenvectors.
dagger()).norm();
535 auto max_delta = delta.
max();
537 report_pass(
"Eigenvalue analysis with " + hila::prettyprint(myMatrix::rows()) +
"x" +
538 hila::prettyprint(myMatrix::columns()) +
" Hermitean matrix",
545 delta[X] = (M[X] - r.U * r.singularvalues * r.V.
dagger()).norm();
548 max_delta = delta.
max();
550 report_pass(
"SVD with " + hila::prettyprint(myMatrix::rows()) +
"x" +
551 hila::prettyprint(myMatrix::columns()) +
" Complex matrix",
559 auto r = M[X].svd_pivot(hila::sort::ascending);
560 delta[X] = (M[X] - r.U * r.singularvalues * r.V.
dagger()).norm();
563 max_delta = delta.
max();
565 report_pass(
"Fully pivoted SVD with " + hila::prettyprint(myMatrix::rows()) +
"x" +
566 hila::prettyprint(myMatrix::columns()) +
" Complex matrix",
575int main(
int argc,
char **argv) {
582 long seed = par.get(
"random seed");
588 lattice.
setup(lsize);
605 test_set_elements_and_select();
611 spectraldensity_test();
613 test_matrix_algebra();
hila::arithmetic_type< T > squarenorm(const Array< n, m, T > &rhs)
Return square norm of Array.
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Field< A > dagger() const
Returns dagger or Hermitian conjugate of Field depending on how it is defined for Field type T.
T max(Parity par=ALL) const
Find maximum value from Field.
Field< Complex< hila::arithmetic_type< T > > > FFT_real_to_complex(fft_direction fdir=fft_direction::forward) const
Field< T > conj() const
Returns field with all elements conjugated depending how conjugate is defined for type.
double squarenorm() const
Squarenorm.
std::vector< T > get_elements(const std::vector< CoordinateVector > &coord_list, bool broadcast=false) const
Get a list of elements and store them into a vector.
T min(Parity par=ALL) const
Find minimum value from Field.
std::vector< T > get_slice(const CoordinateVector &c, bool broadcast=false) const
Get a slice (subvolume)
hila::arithmetic_type< T > norm() const
Calculate vector norm - sqrt of squarenorm.
const Mtype & fill(const S rhs)
Matrix fill.
Rtype transpose() const
Transpose of matrix.
R dot(const Matrix< p, q, S > &rhs) const
Dot product.
Matrix class which defines matrix operations.
Running index for locating sites on the lattice.
std::vector< double > spectraldensity(const Field< T > &f)
k_binning & k_max(double km)
Max value of k.
std::vector< T > bin_k_field(const Field< T > &f)
Generic k-space field binner routine - returns field element type vector.
k_binning & bins(int n)
Set number of bins in histogram.
void setup(const CoordinateVector &siz)
General lattice setup.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
#define foralldir(d)
Macro to loop over (all) Direction(s)
constexpr Parity ALL
bit pattern: 011
constexpr T sqr(const T &arg)
Define convenience function sqr(), returning square of argument.
Vector< 4, double > convert_to_k(const CoordinateVector &cv)
Convert momentum space CoordinateVector to wave number k, where -pi/2 < k_i <= pi/2.
void FFT_field(const Field< T > &input, Field< T > &result, const CoordinateVector &directions, fft_direction fftdir=fft_direction::forward)
int myrank()
rank of this node
double random()
Real valued uniform random number generator.
double gaussrand()
Gaussian random generation routine.
std::ostream out0
This writes output only from main process (node 0)
void initialize(int argc, char **argv)
Read in command line arguments. Initialise default stream and MPI communication.
void seed_random(uint64_t seed, bool device_rng=true)
Seed random generators with 64-bit unsigned value. On MPI shuffles the seed so that different MPI ran...
T broadcast(T &var, int rank=0)
Broadcast the value of var to all MPI ranks from rank (default=0).
void finishrun()
Normal, controlled exit - all nodes must call this. Prints timing information and information about c...