From 96026b91d83fe79bba92e15ba94c9ab5181bf181 Mon Sep 17 00:00:00 2001
From: "lisa.pizzo" <lisa.pizzo@noreply.imsc.uni-graz.at>
Date: Tue, 16 Dec 2025 11:03:17 +0100
Subject: [PATCH] Main File

---
 Sheet5/Ex5_Second_Attempt/main.cpp | 130 +++++++++++++++++++++++++++++
 1 file changed, 130 insertions(+)
 create mode 100644 Sheet5/Ex5_Second_Attempt/main.cpp

diff --git a/Sheet5/Ex5_Second_Attempt/main.cpp b/Sheet5/Ex5_Second_Attempt/main.cpp
new file mode 100644
index 0000000..4162e91
--- /dev/null
+++ b/Sheet5/Ex5_Second_Attempt/main.cpp
@@ -0,0 +1,130 @@
+// Exercise 5 – Parallelization of Exercise 3
+
+#include "geom.h"
+#include "getmatrix.h" //in here it is defined the parallel function 
+#include "jacsolve.h" //in here it is defined the parallel function 
+#include "userset.h"
+#include "vdop.h"
+
+#include <chrono>
+#include <cmath>
+#include <iostream>
+#include <omp.h>
+
+using namespace std;
+using namespace std::chrono;
+
+int main(int, char **)
+{
+    // ------------------------------------------------------------
+    // OpenMP info
+    // ------------------------------------------------------------
+    #pragma omp parallel
+    {
+        #pragma omp single
+        cout << "Using " << omp_get_num_threads()
+             << " OpenMP threads\n\n";
+    }
+
+    // ------------------------------------------------------------
+    // Domain decomposition (MPI-style layout kept, MPI=1)
+    // ------------------------------------------------------------
+    const int numprocs = 1;
+    const int myrank   = 0;
+
+    if (myrank == 0)
+        cout << "There are " << numprocs << " processes running.\n\n";
+
+    const int procx = static_cast<int>(sqrt(numprocs + 0.0));
+    const int procy = procx;
+
+    if (procx * procy != numprocs)
+    {
+        cout << "Wrong number of processors!\n";
+        return 1;
+    }
+
+    // ------------------------------------------------------------
+    // Problem size
+    // ------------------------------------------------------------
+    int nx = 100;
+    int ny = 100;
+    int NXglob = nx * procx;
+    int NYglob = ny * procy;
+
+    cout << "Intervals: " << NXglob << " x " << NYglob << "\n";
+
+// ##################### STL ###########################################
+        {
+        Mesh_2d_3_square const mesh(nx, ny);
+        //mesh.Debug();
+
+        CRS_Matrix SK(mesh);                  // sparse FEM matrix
+        //SK.Debug();
+        vector<double> u(SK.Nrows(), 0.0);    // solution
+        vector<double> f(SK.Nrows(), 0.0);    // RHS
+
+        // Parallel FEM matrix assembly
+        auto t0 = high_resolution_clock::now();
+
+        SK.CalculateLaplace_omp(f);   // PARALLEL assembly
+        //SK.Debug();
+
+        auto t1 = high_resolution_clock::now();
+        cout << "FEM assembly time (parallel): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+
+        // Initial condition
+        mesh.SetValues(u, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
+
+        // Apply Dirichlet BC
+        SK.ApplyDirichletBC(u, f);
+        //SK.Debug();
+
+        // Parallel Jacobi solver
+        t0 = high_resolution_clock::now();
+
+        JacobiSolve_omp(SK, f, u);    // PARALLEL Jacobi
+
+        t1 = high_resolution_clock::now();
+        cout << "JacobiSolve time (parallel): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+    }
+
+    // ============================================================
+    {
+        Mesh_2d_3_matlab const mesh("square_100.txt");
+
+        CRS_Matrix SK(mesh);
+        //SK.Debug();
+        vector<double> u(SK.Nrows(), 0.0);
+        vector<double> f(SK.Nrows(), 0.0);
+
+        auto t0 = high_resolution_clock::now();
+
+        SK.CalculateLaplace_omp(f);   // PARALLEL assembly
+        //SK.Debug();
+
+        auto t1 = high_resolution_clock::now();
+        cout << "FEM assembly time (parallel, matlab mesh): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+
+        mesh.SetValues(u, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
+        SK.ApplyDirichletBC(u, f);
+        //SK.Debug();
+
+        t0 = high_resolution_clock::now();
+
+        JacobiSolve_omp(SK, f, u);    // PARALLEL Jacobi
+
+        t1 = high_resolution_clock::now();
+        cout << "JacobiSolve time (parallel, matlab mesh): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+    }
+
+    return 0;
+}