diff --git a/meas/Makefile.in b/meas/Makefile.in
index 1d04ce4e8..fa7333375 100644
--- a/meas/Makefile.in
+++ b/meas/Makefile.in
@@ -37,7 +37,8 @@ libmeas_TARGETS = measurements \
pion_norm \
polyakov_loop \
measure_clover_field_strength_observables \
- gradient_flow
+ gradient_flow \
+ eigenvalues
libmeas_OBJECTS = $(addsuffix .o, ${libmeas_TARGETS})
diff --git a/meas/eigenvalues.c b/meas/eigenvalues.c
new file mode 100644
index 000000000..34d8a204c
--- /dev/null
+++ b/meas/eigenvalues.c
@@ -0,0 +1,101 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2024 Bartosz Kostrzewa
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with tmLQCD. If not, see .
+ *
+ ************************************************************************/
+
+#ifdef HAVE_CONFIG_H
+# include
+#endif
+#ifdef TM_USE_OMP
+# include
+#endif
+
+#include
+#include
+
+#include "global.h"
+#include "geometry_eo.h"
+#include "fatal_error.h"
+#include "measurements.h"
+#include "operator.h"
+#include "gettime.h"
+#include "tm_debug_printf.h"
+#include "misc_types.h"
+
+#ifdef TM_USE_QUDA
+#include "quda_interface.h"
+#endif
+
+void eigenvalues_measurement(const int traj, const int id, const int ieo) {
+ tm_stopwatch_push(&g_timers, __func__, "");
+ init_operators();
+
+ if(no_operators < 1){
+ tm_debug_printf(0, 0, "Error: no operators defined in input file, cannot perform eigenvalues online measurement!\n");
+#ifdef TM_USE_MPI
+ MPI_Finalize();
+#endif
+ exit(1);
+ }
+
+ eig_param_t eig = measurement_list[id].eig;
+
+ double * evals = calloc(eig.n_evals, sizeof(double));
+
+ for( int op_id = 0; op_id < no_operators; op_id++ ){
+ operator * optr = &operator_list[op_id];
+
+ if( !(optr->type == TMWILSON || optr->type == WILSON ||
+ optr->type == CLOVER || optr->type == DBTMWILSON ||
+ optr->type == DBCLOVER) ){
+ tm_debug_printf(0, 0, "Error: only operator types WILSON, TMWILSON, CLOVER, DBTMWILSON and DBCLOVER are supported.\n");
+#ifdef TM_USE_MPI
+ MPI_Finalize();
+#endif
+ exit(1);
+ }
+
+ op_backup_restore_globals(TM_BACKUP_GLOBALS);
+ op_set_globals(op_id);
+
+ if( measurement_list[id].external_library == QUDA_LIB ){
+#ifdef TM_USE_QUDA
+ eigsolveQuda(evals, eig.n_evals, eig.tol, 1, 0, eig.max_iter, eig.maxmin,
+ optr->eps_sq, optr->maxiter, eig.polydeg, eig.amin, eig.amax, eig.n_kr,
+ optr->solver, optr->solver, optr->rel_prec, ieo,
+ optr->sloppy_precision, optr->sloppy_precision, optr->compression_type,
+ optr->type == CLOVER || optr->type == TMWILSON || optr->type == WILSON);
+#else
+ tm_debug_printf(0, 0, "Error: Attempted to use QUDA eigensolver but this build was not configured for QUDA usage.\n");
+#ifdef TM_USE_MPI
+ MPI_Finalize();
+#endif
+ exit(1);
+#endif
+ }
+
+ op_backup_restore_globals(TM_RESTORE_GLOBALS);
+
+ } // loop over operators
+
+ free(evals);
+ tm_stopwatch_pop(&g_timers, 0, 1, "");
+ return;
+}
+
diff --git a/meas/eigenvalues.h b/meas/eigenvalues.h
new file mode 100644
index 000000000..fbf24284c
--- /dev/null
+++ b/meas/eigenvalues.h
@@ -0,0 +1,32 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2024 Bartosz Kostrzewa
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with tmLQCD. If not, see .
+ ***********************************************************************/
+
+#ifndef _MEASURE_EIGENVALUES_H
+#define _MEASURE_EIGENVALUES_H
+
+#include "su3.h"
+
+/* implements the eigenvalue online measurement
+ * writes into eigenvalues.mXY.IJKLMN
+ * where XY corresponds to the operator index and IJKLMN refers to the trajectory counter */
+
+void eigenvalues_measurement(const int traj, const int id, const int ieo);
+
+#endif
diff --git a/meas/measurements.c b/meas/measurements.c
index b18d263f9..70862cfac 100644
--- a/meas/measurements.c
+++ b/meas/measurements.c
@@ -30,6 +30,7 @@
#include "default_input_values.h"
#include "read_input.h"
+#include "eigenvalues.h"
#include "pion_norm.h"
#include "correlators.h"
#include "polyakov_loop.h"
@@ -80,6 +81,10 @@ int init_measurements(){
measurement_list[i].measurefunc = &gradient_flow_measurement;
}
+ if(measurement_list[i].type == EIGENVALUES) {
+ measurement_list[i].measurefunc = &eigenvalues_measurement;
+ }
+
measurement_list[i].id = i;
}
return(0);
diff --git a/meas/measurements.h b/meas/measurements.h
index 0f5d6d6ab..810c38d2a 100644
--- a/meas/measurements.h
+++ b/meas/measurements.h
@@ -22,6 +22,8 @@
#ifndef _MEASUREMENTS_H
#define _MEASUREMENTS_H
+#include "misc_types.h"
+
#define max_no_measurements 20
/* Give the measurement types an unambiguous ID*/
@@ -30,9 +32,23 @@ enum MEAS_TYPE {
PIONNORM,
POLYAKOV,
ORIENTED_PLAQUETTES,
- GRADIENT_FLOW
+ GRADIENT_FLOW,
+ EIGENVALUES
};
+typedef struct {
+ int n_evals;
+ int blksize;
+ int blkwise;
+ int max_iter;
+ double tol;
+ int maxmin;
+ int polydeg;
+ double amin;
+ double amax;
+ int n_kr;
+} eig_param_t;
+
typedef struct {
enum MEAS_TYPE type;
int initialised;
@@ -70,6 +86,8 @@ typedef struct {
void (*measurefunc) (const int traj, const int id, const int ieo);
ExternalLibrary external_library;
+
+ eig_param_t eig;
} measurement;
diff --git a/read_input.l b/read_input.l
index 3b39524dc..24f8cbefe 100644
--- a/read_input.l
+++ b/read_input.l
@@ -553,6 +553,7 @@ static inline double fltlist_next_token(int * const list_end){
%x PLOOP
%x ORIENTEDPLAQUETTESMEAS
%x GRADIENTFLOWMEAS
+%x EIGENVALUESMEAS
%x REWEIGH
%x REWSAMPLES
@@ -3058,6 +3059,10 @@ static inline double fltlist_next_token(int * const list_end){
meas->type = GRADIENT_FLOW;
strcpy(meas->name, "GRADIENTFLOW");
}
+ else if(strcmp(yytext, "EIGENVALUES")==0) {
+ meas->type = EIGENVALUES;
+ strcpy(meas->name, "EIGENVALUES");
+ }
else {
fprintf(stderr, "Unknown measurement type %s in line %d\n", yytext, line_of_file);
exit(1);
@@ -3086,9 +3091,10 @@ static inline double fltlist_next_token(int * const list_end){
else if(meas->type == POLYAKOV) BEGIN(PLOOP);
else if(meas->type == ORIENTED_PLAQUETTES) BEGIN(ORIENTEDPLAQUETTESMEAS);
else if(meas->type == GRADIENT_FLOW) BEGIN(GRADIENTFLOWMEAS);
+ else if(meas->type == EIGENVALUES) BEGIN(EIGENVALUESMEAS);
}
-{
+{
^EndMeasurement{SPC}* {
if(myverbose) printf("Measurement with id %d parsed in line %d\n\n", meas->id, line_of_file);
BEGIN(0);
@@ -3142,6 +3148,65 @@ static inline double fltlist_next_token(int * const list_end){
}
}
+{
+ {SPC}*Tolerance{EQL}{FLT} {
+ sscanf(yytext, " %[2a-zA-Z] = %lf", name, &c);
+ meas->eig.tol = c;
+ if(myverbose) printf(" Eigenvalue measurement tolerance set to %lf line %d, measurement id=%d\n", meas->eig.tol, line_of_file, meas->id);
+ }
+ {SPC}*NoEigenvalues{EQL}{DIGIT}+ {
+ sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+ meas->eig.n_evals = a;
+ if(myverbose) printf(" Number of eigenvalues to be determined set to %d, line %d measurement id=%d\n", a, line_of_file, meas->id);
+ }
+ {SPC}*BlockSize{EQL}{DIGIT}+ {
+ sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+ meas->eig.blksize = a;
+ if(myverbose) printf(" Block size for eigensolver set to %d, line %d measurement id=%d\n", a, line_of_file, meas->id);
+ }
+ {SPC}*Blockwise{EQL}yes {
+ meas->eig.blkwise = 1;
+ if(myverbose) printf(" Blockwise eigenvalue determination, line %d, measurement id=%d\n", line_of_file, meas->id);
+ }
+ {SPC}*Blockwise{EQL}no {
+ meas->eig.blkwise = 0;
+ if(myverbose) printf(" NO blockwise eigenvalue deterination, line %d, measurement id=%d\n", line_of_file, meas->id);
+ }
+ {SPC}*LargestEigenvalues{EQL}yes {
+ meas->eig.maxmin = 1;
+ if(myverbose) printf(" Largest eigenvalues to be determined, line %d, measurement id=%d\n", line_of_file, meas->id);
+ }
+ {SPC}*LargestEigenvalues{EQL}no {
+ meas->eig.maxmin = 0;
+ if(myverbose) printf(" Smallest eigenvalues to be determined, line %d, measurement id=%d\n", line_of_file, meas->id);
+ }
+ {SPC}*MaxIterations{EQL}{DIGIT}+ {
+ sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+ meas->eig.max_iter = a;
+ if(myverbose) printf(" Maximum number of iterations for eigensolver set to %d, line %d measurement id=%d\n", a, line_of_file, meas->id);
+ }
+ {SPC}*PolynomialDegree{EQL}{DIGIT}+ {
+ sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+ meas->eig.polydeg = a;
+ if(myverbose) printf(" Degree of polynomial for eigensolver set to %d, line %d measurement id=%d\n", a, line_of_file, meas->id);
+ }
+ {SPC}*PolyMin{EQL}{FLT} {
+ sscanf(yytext, " %[2a-zA-Z] = %lf", name, &c);
+ meas->eig.amin = c;
+ if(myverbose) printf(" Minimum eigenvalue to exclude using polynomial acceleration in eigensolver set to %e line %d, measurement id=%d\n", c, line_of_file, meas->id);
+ }
+ {SPC}*PolyMax{EQL}{FLT} {
+ sscanf(yytext, " %[2a-zA-Z] = %lf", name, &c);
+ meas->eig.amax = c;
+ if(myverbose) printf(" Maximum eigenvalue to exclude using polynomial acceleration in eigensolver set to %e line %d, measurement id=%d\n", c, line_of_file, meas->id);
+ }
+ {SPC}*KrylovSubspaceSize{EQL}{DIGIT}+ {
+ sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+ meas->eig.n_kr = a;
+ if(myverbose) printf(" Krylov subspace size for polynomial acceleration in eigensolver set to %d, line %d measurement id=%d\n", a, line_of_file, meas->id);
+ }
+}
+
{
{SPC}*Direction{EQL}[03] {
sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
@@ -3691,7 +3756,7 @@ static inline double fltlist_next_token(int * const list_end){
BEGIN(comment_caller);
}
-{SPC}*\n {
+{SPC}*\n {
line_of_file++;
}
<*>{SPC}*\n {
diff --git a/sample-input/sample-eigenvalues-measurement.input b/sample-input/sample-eigenvalues-measurement.input
new file mode 100644
index 000000000..80965e097
--- /dev/null
+++ b/sample-input/sample-eigenvalues-measurement.input
@@ -0,0 +1,50 @@
+# example input file for eigenvalue measurements using "offline_measurement"
+# requires 2 8^4 gauge configuration conf.0000 and conf.0002
+
+L=8
+T=8
+
+DebugLevel = 5
+ompnumthreads=4
+
+InitialStoreCounter = 0
+Measurements = 2
+# measurements will be carried out in nsave steps
+# e.g. for conf.0000 and conf.0002 in this case
+nsave=2
+2kappamu = 0.05
+kappa = 0.177
+BCAngleT = 1
+GaugeConfigInputFile = conf
+UseEvenOdd = yes
+
+# the eigenvalues measurement requires at least ONE operator to be defined
+# if multiple operators are defined, the measurement will loop over them
+# and produce numbered output files
+BeginMeasurement EIGENVALUES
+ Frequency = 1
+ NoEigenvalues = 10
+ LargestEigenvalues = yes
+ BlockSize = 1
+ Blockwise = no
+ MaxIterations = 10000
+ Tolerance = 1e-14
+ PolynomialDegree = 0
+ PolyMin = 0.0002
+ PolyMax = 0.003
+ KrylovSubspaceSize = 100
+ UseExternalLibrary = quda
+EndMeasurement
+
+# note: setting the solver to CGMMS will result in the CGMMS inversion taking place
+# because the solver is not properly decoupled form the rest of the code
+BeginOperator TMWILSON
+ 2kappaMu = 0.05
+ kappa = 0.177
+ UseEvenOdd = yes
+ Solver = CG
+ SolverPrecision = 1e-14
+ MaxSolverIterations = 1000
+ AddDownPropagator = no
+EndOperator
+