From 05a1cc7f52c40dcac613c5321b0e516ce33b6907 Mon Sep 17 00:00:00 2001
From: "Pyry Runko (Puck)" <pyjorunk@student.jyu.fi>
Date: Wed, 26 Feb 2025 09:09:07 +0200
Subject: [PATCH] Changed double Amplitude to std::complex<double>
 ComplexAmplitude, contains old fundouble Amplitude to std::complex<double>
 ComplexAmplitude, contains old functions to be deletedd

---
 CMakeLists.txt   |  3 +++
 src/ipglasma.cpp | 44 +++++++++++++++++++++-----------------------
 src/ipglasma.hpp |  8 ++++----
 src/main.cpp     | 10 +++++++++-
 4 files changed, 37 insertions(+), 28 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0a084dd..4056a44 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 3a3d232..b356be9 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 354240c..1da59da 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 46f1543..064ae29 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");
-- 
GitLab