diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0a084dd7900afebacfaf4ca2c585f1fb05ac14ad..4056a443f56a3efe6b49d8b696f16cfb1f12fed3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,9 +14,12 @@ include_directories(${GSL_INCLUDE_DIRS})
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIRS})
+find_package(OpenMP REQUIRED)
add_subdirectory(src)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
set_target_properties(wline PROPERTIES
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin" )
diff --git a/src/ipglasma.cpp b/src/ipglasma.cpp
index 3a3d232572652fb83cc3d78cfd1792e3f42309d5..b356be9614a0bb3331e48afa2dd78d88c21aeecc 100644
--- a/src/ipglasma.cpp
+++ b/src/ipglasma.cpp
@@ -475,8 +475,8 @@ int IPGlasma::WilsonLineCoordinate(int xind, int yind)
}
-// attempt to overload amplitude with one wilson line and two indices
-double IPGlasma::Amplitude(double xpom, WilsonLine quark, int q2[2])
+// compute complex amplitude with one wilson line and two indices
+std::complex<double> IPGlasma::ComplexAmplitude(double xpom, WilsonLine quark, int q2[2])
{
//TODO: make work with int?
//ApplyPeriodicBoundaryConditions(q1);
@@ -503,27 +503,22 @@ double IPGlasma::Amplitude(double xpom, WilsonLine quark, int q2[2])
cout << antiquark << endl;
exit(1);
}
- std::complex<double > amp = 1.0 - 1.0/NC * prod.Trace();
+ std::complex<double> amp = 1.0 - 1.0/NC * prod.Trace();
double result = amp.real();
- if (result < 0) return 0;
- return result;
- if (result > 1)
- return 1;
- if (result < 0)
- return 0;
-
- if (isnan(result))
+ //if (result < 0) return 0;
+ //if (result > 1) return 1;
+ //if (isnan(result))
{
cerr << "Wilson line trance NaN, antiquark coords " << q2[0] << ", " << q2[1] << "and quark wline" << quark << endl;
exit(1);
}
- return result;
+ return amp;
}
-// attempt to overload amplitude when given coordinates are integers i.e. indices
-double IPGlasma::Amplitude(double xpom, int q1[2], int q2[2])
+// overload amplitude when given coordinates are integers i.e. indices
+std::complex<double> IPGlasma::ComplexAmplitude(double xpom, int q1[2], int q2[2])
{
//TODO: make work with int?
//ApplyPeriodicBoundaryConditions(q1);
@@ -586,9 +581,9 @@ double IPGlasma::DistanceToOrigin(int q[2]){
//TODO: check
-double IPGlasma::SumAmplitudes(int w)
+std::complex<double> IPGlasma::SumAmplitudes(int w)
{
- double sum = 0.0;
+ std::complex<double> sum = 0.0;
int L = xcoords.size();
cout << "# Summing with lattice size L = " << L << endl;
WilsonLine quark;
@@ -609,7 +604,7 @@ double IPGlasma::SumAmplitudes(int w)
continue;
}
- sum += Amplitude(0, quark, q2);
+ sum += ComplexAmplitude(0, quark, q2);
}
}
}
@@ -621,7 +616,7 @@ double IPGlasma::SumAmplitudes(int w)
}
//TODO: check
-double IPGlasma::SumAmplitudes()
+std::complex<double> IPGlasma::SumAmplitudes()
{
int L = xcoords.size();
cout << "# Summing with lattice size L = " << L << endl;
@@ -663,9 +658,12 @@ double IPGlasma::SumAmplitudes()
std::complex<double> sum = 0.0;
+ //double sum = 0.0;
+
- //TODO: does pragma omp work?
- //#pragma omp parallel for reduction(+:sum)
+ #pragma omp declare reduction(+:std::complex<double>:omp_out+=omp_in)
+ //TODO: check that this sums correctly
+ #pragma omp parallel for reduction(+:sum)
//error: user defined reduction not found for ‘sum’
for(int i =0; i < L; i += 1){
@@ -693,7 +691,7 @@ double IPGlasma::SumAmplitudes()
double WF_d = WFGammaqqbarLongitudinal(r, z0, Q, d_quark);
double WF_c = WFGammaqqbarLongitudinal(r, z0, Q, c_quark);
- sum += (WF_u + WF_d + WF_c) * Amplitude(0, quark, q2)
+ sum += (WF_u + WF_d + WF_c) * ComplexAmplitude(0, quark, q2)
* exp( -i * ( std::inner_product(std::begin(vecL), std::end(vecL), std::begin(diffX), 0)
+ std::inner_product(std::begin(vecDelta), std::end(vecDelta), std::begin(sumZX), 0) ) );
@@ -709,8 +707,8 @@ double IPGlasma::SumAmplitudes()
}
}
}
- cout << "# total amplitude squared = " << std::norm(sum) << endl;
- return std::norm(sum);
+ cout << "# total amplitude = " << sum << endl;
+ return sum;
}
diff --git a/src/ipglasma.hpp b/src/ipglasma.hpp
index 354240c0e505b8cebe230abf63fd1eb1619f66e2..1da59da3353e0f64580bfa894f6c773bc1c48000 100644
--- a/src/ipglasma.hpp
+++ b/src/ipglasma.hpp
@@ -57,10 +57,10 @@ public:
//Pyry added:
double DistanceToOrigin(int q[2]);
- double Amplitude(double xpom, int q1[2], int q2[2]);
- double Amplitude(double xpom, WilsonLine quark, int q2[2]);
- double SumAmplitudes();
- double SumAmplitudes(int w);
+ std::complex<double> ComplexAmplitude(double xpom, int q1[2], int q2[2]);
+ std::complex<double> ComplexAmplitude(double xpom, WilsonLine quark, int q2[2]);
+ std::complex<double> SumAmplitudes();
+ std::complex<double> SumAmplitudes(int w);
struct quark_flavor{
double mass;
diff --git a/src/main.cpp b/src/main.cpp
index 46f1543b349b7d4a2ed0c96765977a9fa8153225..064ae291aa1d43ae4bf0dcef7c220fe7f8e21ea2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -58,6 +58,14 @@ int printDipoleAmplitudes(string fname, IPGlasma ipglasma)
int main(int argc, char* argv[])
{
+ //TODO if num_threads > 1 cout << "using omp";
+
+ if (omp_get_max_threads() > 1){
+ cout << "Running with OpenMP." << " Using " << omp_get_max_threads() << " threads." << std::endl;
+ }
+ else { cout << "Running without OpenMP." << endl;}
+
+
if (argc < 3 || argc > 4){
cerr << "Usage: " << argv[0] << " filename 1(integrate)/0(print amplitudes) [task_id]" << endl;
}
@@ -78,7 +86,7 @@ int main(int argc, char* argv[])
else{
- double sum = ipglasma.SumAmplitudes();
+ std::complex<double> sum = ipglasma.SumAmplitudes();
// Save results
//std::filesystem::create_directory("output_"+fname);
ofstream outfile("da_sum_ +" + fname + ".txt");