now working, macos (m1 chip) compilation falls back to sequential LOL

CMakeLists.txt

@@ -1,18 +1,27 @@
-cmake_minimum_required(VERSION 3.10)
+cmake_minimum_required(VERSION 3.20)
 
 # Set the project name
 project(OpenQCD++)
 
 # Specify the C++ standard
-set(CMAKE_CXX_STANDARD 20)
-set(CMAKE_CXX_STANDARD_REQUIRED True)
+set(CMAKE_CXX_STANDARD 20)  # This is necessary for the std::execution policies.
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Specify the compilers
+set(CMAKE_C_COMPILER "clang")
+set(CMAKE_CXX_COMPILER "clang++")
+
+# Ensure the GCC standard library is used
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++20")
+
+
 
 # Add the executables
-add_executable(init_print_lattice app/init_print_lattice/main.cpp modules/lattice/LatticeField.cpp)
-add_executable(randomize app/randomize/main.cpp modules/lattice/LatticeField.cpp)
-add_executable(compute_plaquette app/compute_plaquette/main.cpp modules/lattice/LatticeField.cpp modules/gauge/GaugeField.cpp)
+add_executable(randomize app/parallel_computations/randomize.cpp modules/lattice/LatticeField.cpp)
+add_executable(apply_func app/parallel_computations/apply_func.cpp modules/lattice/LatticeField.cpp)
+add_executable(reduce_if app/parallel_computations/reduce_if.cpp modules/lattice/LatticeField.cpp)
 
 # Include directories
-target_include_directories(init_print_lattice PRIVATE include)
 target_include_directories(randomize PRIVATE include)
-target_include_directories(compute_plaquette PRIVATE include)
+target_include_directories(apply_func PRIVATE include)
+target_include_directories(reduce_if PRIVATE include)

app/parallel_computations/apply_func.cpp

@@ -0,0 +1,56 @@
+// This program performs some examples of the parallel computations that can be done with the
+// LatticeField class.
+// As examples, this performs and times:
+// - Parallel application of a function (an exp) to each element of the lattice field.
+
+#include "lattice/LatticeField.h"
+#include <complex>
+#include <chrono>
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <vector>
+#include <functional>
+
+int main() {
+    using namespace std::chrono;
+
+    std::ofstream csv_file("apply_func_times.csv");
+    csv_file << "Parallel_exp [ms], Sequential_exp [ms]\n";
+
+    for (int i = 0; i < 100; ++i) {
+        // Create a lattice field of size 8x8x8x8 with initial value 0.0
+        LatticeField<double> lattice(8, 8, 8, 8, 0.0);
+
+        // Randomize the lattice field
+        lattice.randomize();
+
+        // Apply exp function to each element of the lattice field in parallel
+        const auto start = high_resolution_clock::now();
+        lattice.func([](double& val) {
+            val = std::exp(val);
+        });
+        const auto end = high_resolution_clock::now();
+        auto timed_parallel_exp = duration_cast<microseconds>(end - start).count();
+
+        // Randomize the lattice field
+        lattice.randomize();
+
+        // Apply exp function to each element of the lattice field in parallel
+        const auto start_s = high_resolution_clock::now();
+        lattice.func_sequential([](double& val) {
+            val = std::exp(val);
+        });
+        const auto end_s = high_resolution_clock::now();
+        auto timed_seq_exp = duration_cast<microseconds>(end_s - start_s).count();
+
+
+        // Write the times to the CSV file
+        csv_file << timed_parallel_exp << ", " << timed_seq_exp << "\n";
+    }
+
+    csv_file.close();
+    std::cout << "Timing results have been written to apply_func_times.csv" << std::endl;
+
+    return 0;
+}

app/parallel_computations/randomize.cpp

@@ -0,0 +1,45 @@
+// This program performs some examples of the parallel computations that can be done with the
+// LatticeField class.
+// As examples, this performs and times:
+// - Randomization of a lattice field.
+
+#include "lattice/LatticeField.h"
+#include <complex>
+#include <chrono>
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <vector>
+#include <functional>
+
+int main() {
+    using namespace std::chrono;
+
+    std::ofstream csv_file("randomize_times.csv");
+    csv_file << "Randomization_par_unseq [ms], Randomization_seq [ms]\n";
+
+    for (int i = 0; i < 100; ++i) {
+        // Create a lattice field of size 4x4x4x4 with initial value 0.0
+        LatticeField<double> lattice(8, 8, 8, 8, 0.0);
+
+        // Parallel Randomize the lattice field
+        const auto start = high_resolution_clock::now();
+        lattice.randomize();
+        const auto end = high_resolution_clock::now();
+        auto timed_randomize_parallel = duration_cast<microseconds>(end - start).count();
+
+        // Randomize the lattice field
+        auto start_time = high_resolution_clock::now();
+        lattice.randomize_sequential();
+        auto end_time = high_resolution_clock::now();
+        auto timed_randomize_sequential = duration_cast<microseconds>(end_time - start_time).count();
+
+        // Write the times to the CSV file
+        csv_file << timed_randomize_parallel << "," << timed_randomize_sequential << "\n";
+    }
+
+    csv_file.close();
+    std::cout << "Timing results have been written to timing_results.csv" << std::endl;
+
+    return 0;
+}

app/parallel_computations/reduce_if.cpp

@@ -0,0 +1,60 @@
+// This program performs some examples of the parallel computations that can be done with the
+// LatticeField class.
+// As examples, this performs and times:
+// - Parallel reduction of the lattice field (values greater than 0.5 only).
+
+
+#include "lattice/LatticeField.h"
+#include <complex>
+#include <chrono>
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <vector>
+#include <functional>
+
+int main() {
+    using namespace std::chrono;
+
+    std::ofstream csv_file("reduce_if_times.csv");
+    csv_file << "Parallel_reduce_if [ms],Sequential_reduce_if [ms]\n";
+
+    for (int i = 0; i < 100; ++i) {
+        // Create a lattice field of size 8x8x8x8 with initial value 0.0
+        LatticeField<double> lattice(8, 8, 8, 8, 0.0);
+
+        // Randomize the lattice field
+        lattice.randomize();
+
+        // Perform parallel reduction of the lattice field (values greater than 0.5 only)
+        const auto start_parallel = high_resolution_clock::now();
+        double sum_parallel = lattice.reduce_if(
+            [](const double& val) { return val > 0.5; },
+            0.0,
+            [](const double& a, const double& b) { return a + b; }
+        );
+        const auto end_parallel = high_resolution_clock::now();
+        auto timed_parallel_reduce = duration_cast<microseconds>(end_parallel - start_parallel).count();
+
+        // Randomize the lattice field
+        lattice.randomize();
+
+        // Perform sequential reduction of the lattice field (values greater than 0.5 only)
+        const auto start_sequential = high_resolution_clock::now();
+        double sum_sequential = lattice.reduce_if_sequential(
+            [](const double& val) { return val > 0.5; },
+            0.0,
+            [](const double& a, const double& b) { return a + b; }
+        );
+        const auto end_sequential = high_resolution_clock::now();
+        auto timed_sequential_reduce = duration_cast<microseconds>(end_sequential - start_sequential).count();
+
+        // Write the times to the CSV file
+        csv_file << timed_parallel_reduce << "," << timed_sequential_reduce << "\n";
+    }
+
+    csv_file.close();
+    std::cout << "Timing results have been written to reduce_if_times.csv" << std::endl;
+
+    return 0;
+}

include/gauge/GaugeField.h

@@ -1,16 +1 @@
-#ifndef GAUGEFIELD_H
-#define GAUGEFIELD_H
-
-#include "lattice/LatticeField.h" // Include the header for LatticeField
-
-// Define the subclass GaugeField
-template <typename T>
-class GaugeField : public LatticeField<T> {
-public:
-    GaugeField(int dim1, int dim2, int dim3, int dim4, const T& default_value)
-        : LatticeField<T>(dim1, dim2, dim3, dim4, default_value) {};
-    ~GaugeField() = default;
-
-    // Additional methods specific to GaugeField can be added here
-    T compute_plaquette() const;
-};
+// TODO:

include/lattice/LatticeField.h

@@ -4,24 +4,45 @@
 #include <complex>
 #include <iostream>
 #include <vector>
-#include "mdarray/mdarray.h"
+#include <array>
+#include <functional>
 
+// TODO: For Gauge and Spinor field implementation, add an extra parameter of the representation.
+// I.e. this LatticeField is a container for the field values, of type T.
+// The container for Gauge fields would be a container of SU(N) matrices (templated type M), and
+// each matrix entry of type T (e.g., a std::complex<double>). Additionally, it will have an extra
+// set of internal indices for the dirac and color indices.
 template <typename T>
 class LatticeField {
 public:
     LatticeField(
-        const int dim0, const int dim1, const int dim2, const int dim3,
-        const T default_value
+        const size_t dim0, const size_t dim1, const size_t dim2, const size_t dim3,
+        const T init_value
     );
     ~LatticeField() = default;
 
-    void print(size_t i0, size_t i1, size_t i2, size_t i3) const;
-    void randomize();
+    T& operator()(size_t i0, size_t i1, size_t i2, size_t i3);
+    const T& operator()(size_t i0, size_t i1, size_t i2, size_t i3) const;
+    size_t get_flat_idx(size_t i0, size_t i1, size_t i2, size_t i3) const;
+
+    void randomize();  // TODO: add options to randomize with different distributions
+    void randomize_sequential();  // Example of a serial implementation.
+
+    // This will do a parallel application of fun on each element in field_values.
+    void func(std::function<void(T&)> fun);
+    void func_sequential(std::function<void(T&)> fun);
+
+    // This will perform a binary operator between all values in field_values after filtering by a
+    // logic predicate.
+    T reduce_if(std::function<bool(const T&)> predicate, T init, std::function<T(const T&, const T&)> binary_op);
+    T reduce_if_sequential(std::function<bool(const T&)> predicate, T init, std::function<T(const T&, const T&)> binary_op);
+
 
 private:
-    const int dim0, dim1, dim2, dim3;
-    const int volume;
-    std::array<T, volume> field_values;
+    const size_t dim0, dim1, dim2, dim3;
+    const size_t volume;
+    // Maybe in the future we can use a mdvector / mdarray if something like that gets into the std.
+    std::vector<T> field_values;
 };
 
 #endif // LATTICEFIELD_H

include/mdarray/mdarray.h

@@ -49,28 +49,28 @@ using mdarray = typename mdarray_impl<T, sizes...>::type;
 // Function template to fill the mdarray with a default value
 // FIXME: This is buggy.
 // template <typename T, std::size_t dim_i, std::size_t... dim_ixx>
-// void fill_mdarray_with_value(mdarray_impl<T, dim_i, dim_ixx...>& arr, const T& default_value) {
+// void fill_mdarray_with_value(mdarray_impl<T, dim_i, dim_ixx...>& arr, const T& init_value) {
 //     for (auto& sub_arr : arr) {
-//         fill_mdarray_with_value<T, dim_ixx...>(sub_arr, default_value); // Recursive call
+//         fill_mdarray_with_value<T, dim_ixx...>(sub_arr, init_value); // Recursive call
 //     }
 // }
 
 // Function to fill with default value on 4D array
 template <typename T, std::size_t dim_0, std::size_t dim_1, std::size_t dim_2, std::size_t dim_3>
-void fill_mdarray_with_value_4d(mdarray<T, dim_0, dim_1, dim_2, dim_3>& arr, const T& default_value) {
+void fill_mdarray_with_value_4d(mdarray<T, dim_0, dim_1, dim_2, dim_3>& arr, const T& init_value) {
     for (auto& sub_arr_3d : arr) {
         for (auto& sub_arr_2d : sub_arr_3d) {
             for (auto& sub_arr_1d : sub_arr_2d) {
-                std::fill(sub_arr_1d.begin(), sub_arr_1d.end(), default_value); // Fill innermost array
+                std::fill(sub_arr_1d.begin(), sub_arr_1d.end(), init_value); // Fill innermost array
             }
         }
     }
 }
 
 // Base case for the innermost array
 template <typename T, std::size_t dim_i>
-void fill_mdarray_with_value(std::array<T, dim_i>& arr, const T& default_value) {
-    std::fill(arr.begin(), arr.end(), default_value); // Fill innermost array
+void fill_mdarray_with_value(std::array<T, dim_i>& arr, const T& init_value) {
+    std::fill(arr.begin(), arr.end(), init_value); // Fill innermost array
 }
 
 
@@ -98,15 +98,15 @@ void fill_mdarray_with_iterator_4d(mdarray<T, dim_0, dim_1, dim_2, dim_3>& arr,
 
 // // Function template to fill the mdarray with a default value
 // template <typename T, std::size_t size, std::size_t... sizes>
-// void fill_mdarray_with_value(mdarray<T, size, sizes...>& arr, const T& default_value) {
+// void fill_mdarray_with_value(mdarray<T, size, sizes...>& arr, const T& init_value) {
 //     for (auto& sub_arr : arr) {
-//         fill_mdarray_with_value(sub_arr, default_value); // Recursive call
+//         fill_mdarray_with_value(sub_arr, init_value); // Recursive call
 //     }
 // }
 
 // // Base case for the innermost array
 // template <typename T, std::size_t size>
-// void fill_mdarray_with_value(std::array<T, size>& arr, const T& default_value) {
-//     std::fill(arr.begin(), arr.end(), default_value); // Fill innermost array
+// void fill_mdarray_with_value(std::array<T, size>& arr, const T& init_value) {
+//     std::fill(arr.begin(), arr.end(), init_value); // Fill innermost array
 // }
 

include/mdarray/test2.cpp

@@ -28,8 +28,8 @@ int main() {
     constexpr std::size_t dim_4 = 5;
 
     mdarray<int, dim_4> arr_1d;
-    int default_value = 33;
-    fill_mdarray_with_value(arr_1d, default_value);
+    int init_value = 33;
+    fill_mdarray_with_value(arr_1d, init_value);
     std::cout << arr_1d[0] << std::endl; // Should print 33
     return 0;
 }

modules/gauge/GaugeField.cpp

@@ -1,35 +1 @@
-#include "lattice/LatticeField.h" // Include the header for LatticeField
-#include "gauge/GaugeField.h"
-#include <numeric>           // For std::accumulate
-#include <vector>            // For std::vector
-
-// Define the subclass GaugeField
-template <typename T>  // In here, T would be a length 4 / 8 array of SU(3) matrices.
-class GaugeField : public LatticeField {
-public:
-    // Constructor
-    GaugeField() : LatticeField() {}
-
-    // Function to compute the plaquette
-    double compute_plaquette(const std::vector<double>& field_values) const{
-        // Perform a reduction sum over all the values in the field_values
-        return std::accumulate(field_values.begin(), field_values.end(), 0.0);
-    }
-};
-
-
-
-#include "GaugeField.h"
-
-template <typename T>
-GaugeField<T>::GaugeField(int dim1, int dim2, int dim3, int dim4)
-    : LatticeField<T>(dim1, dim2, dim3, dim4) {
-    // Additional initialization specific to GaugeField can be added here
-
-}
-
-// Explicit template instantiation
-template class GaugeField<float>;
-template class GaugeField<double>;
-template class GaugeField<std::complex<float>>;
-template class GaugeField<std::complex<double>>;
+// TODO: