OMP version comparable runtime to mpi

PoissonSolver/2D-omp/Makefile

@@ -0,0 +1,62 @@
+#=======================================================================================
+# Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+# All rights reserved.
+# Use of this source code is governed by a MIT-style
+# license that can be found in the LICENSE file.
+#=======================================================================================
+
+#CONFIGURE BUILD SYSTEM
+TARGET	   = exe-$(TAG)
+BUILD_DIR  = ./$(TAG)
+SRC_DIR    = ./src
+MAKE_DIR   = ./
+Q         ?= @
+
+#DO NOT EDIT BELOW
+include $(MAKE_DIR)/config.mk
+include $(MAKE_DIR)/include_$(TAG).mk
+INCLUDES  += -I$(SRC_DIR)/includes -I$(BUILD_DIR)
+
+VPATH     = $(SRC_DIR)
+ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s,$(wildcard $(SRC_DIR)/*.c))
+OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o,$(wildcard $(SRC_DIR)/*.c))
+CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
+
+${TARGET}: $(BUILD_DIR) $(OBJ)
+	$(info ===>  LINKING  $(TARGET))
+	$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
+
+$(BUILD_DIR)/%.o:  %.c $(MAKE_DIR)/include_$(TAG).mk
+	$(info ===>  COMPILE  $@)
+	$(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
+	$(Q)$(GCC) $(CPPFLAGS) -MT $(@:.d=.o) -MM  $< > $(BUILD_DIR)/$*.d
+
+$(BUILD_DIR)/%.s:  %.c
+	$(info ===>  GENERATE ASM  $@)
+	$(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@
+
+.PHONY: clean distclean tags info asm
+
+clean:
+	$(info ===>  CLEAN)
+	@rm -rf $(BUILD_DIR)
+	@rm -f tags
+
+distclean: clean
+	$(info ===>  DIST CLEAN)
+	@rm -f $(TARGET)
+
+info:
+	$(info $(CFLAGS))
+	$(Q)$(CC) $(VERSION)
+
+asm:  $(BUILD_DIR) $(ASM)
+
+tags:
+	$(info ===>  GENERATE TAGS)
+	$(Q)ctags -R
+
+$(BUILD_DIR):
+	@mkdir $(BUILD_DIR)
+
+-include $(OBJ:.o=.d)

PoissonSolver/2D-omp/README.md

@@ -0,0 +1,48 @@
+# C source skeleton
+
+## Build
+
+1. Configure the toolchain and additional options in `config.mk`:
+```
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+ENABLE_OPENMP ?= false
+
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER
+```
+
+The verbosity options enable detailed output about affinity settings, allocation sizes and timer resolution.
+
+
+2. Build with:
+```
+make
+```
+
+You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set.
+Intermediate build results are located in the `<TOOLCHAIN>` directory.
+
+To output the executed commands use:
+```
+make Q=
+```
+
+3. Clean up with:
+```
+make clean
+```
+to clean intermediate build results.
+
+```
+make distclean
+```
+to clean intermediate build results and binary.
+
+4. (Optional) Generate assembler:
+```
+make asm
+```
+The assembler files will also be located in the `<TOOLCHAIN>` directory.

PoissonSolver/2D-omp/animate.plot

@@ -0,0 +1,15 @@
+set term png size 1024,768 enhanced font ,12
+set datafile separator whitespace
+set grid
+set hidden3d
+set xrange [0:40]
+set yrange [0:40]
+set zrange [-2:2]
+
+input(n) = sprintf("p-%d.dat", n)
+output(n) = sprintf("%03d.png", n)
+
+do for [i=1:50] {
+    set output output(i)
+  splot input(i) matrix using 1:2:3 with lines
+}

PoissonSolver/2D-omp/config.mk

@@ -0,0 +1,9 @@
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+
+#Feature options
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER
+ENABLE_OPENMP += true

PoissonSolver/2D-omp/include_CLANG.mk

@@ -0,0 +1,18 @@
+CC   = clang
+GCC  = cc
+LINKER = $(CC)
+
+ifeq ($(ENABLE_OPENMP),true)
+OPENMP   = -fopenmp
+#OPENMP   = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
+LIBS     = # -lomp
+endif
+
+VERSION  = --version
+CFLAGS   = -Ofast -std=c99 $(OPENMP)
+#CFLAGS   = -Ofast -fnt-store=aggressive  -std=c99 $(OPENMP) #AMD CLANG
+LFLAGS   = $(OPENMP) -lm
+DEFINES  =  -D_GNU_SOURCE
+DEFINES  += -DANIMATE
+# DEFINES  += -DDEBUG
+INCLUDES =

PoissonSolver/2D-omp/include_GCC.mk

@@ -0,0 +1,14 @@
+CC   = gcc
+GCC  = gcc
+LINKER = $(CC)
+
+ifeq ($(ENABLE_OPENMP),true)
+OPENMP   = -fopenmp
+endif
+
+VERSION  = --version
+CFLAGS   = -Ofast -ffreestanding -std=c99 $(OPENMP)
+LFLAGS   = $(OPENMP)
+DEFINES  = -D_GNU_SOURCE
+INCLUDES =
+LIBS     = -lm

PoissonSolver/2D-omp/include_ICC.mk

@@ -0,0 +1,14 @@
+CC   = icc
+GCC  = gcc
+LINKER = $(CC)
+
+ifeq ($(ENABLE_OPENMP),true)
+OPENMP   = -qopenmp
+endif
+
+VERSION  = --version
+CFLAGS   =  -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP)
+LFLAGS   = $(OPENMP)
+DEFINES  = -D_GNU_SOURCE
+INCLUDES =
+LIBS     =

PoissonSolver/2D-omp/poisson.par

@@ -0,0 +1,22 @@
+# Problem specific Data:
+# ---------------------
+
+name poisson
+
+# Geometry Data:
+# -------------
+
+xlength    1.0		# domain size in x-direction
+ylength    1.0		# domain size in y-direction
+imax       6000		# number of interior cells in x-direction
+jmax       6000		# number of interior cells in y-direction
+
+# Pressure Iteration Data:
+# -----------------------
+
+itermax  100000		# maximal number of pressure iteration in one time step
+eps      0.000001	# stopping tolerance for pressure iteration
+rho      0.99999		# relaxation parameter for SOR iteration
+omg      1.2		# relaxation parameter for SOR iteration
+
+#===============================================================================

PoissonSolver/2D-omp/src/allocate.c

@@ -0,0 +1,37 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+
+void* allocate (int alignment, size_t bytesize)
+{
+    int errorCode;
+    void* ptr;
+
+    errorCode =  posix_memalign(&ptr, alignment, bytesize);
+
+    if (errorCode) {
+        if (errorCode == EINVAL) {
+            fprintf(stderr,
+                    "Error: Alignment parameter is not a power of two\n");
+            exit(EXIT_FAILURE);
+        }
+        if (errorCode == ENOMEM) {
+            fprintf(stderr,
+                    "Error: Insufficient memory to fulfill the request\n");
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    if (ptr == NULL) {
+        fprintf(stderr, "Error: posix_memalign failed!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return ptr;
+}

PoissonSolver/2D-omp/src/allocate.h

@@ -0,0 +1,11 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#ifndef __ALLOCATE_H_
+#define __ALLOCATE_H_
+#include <stdlib.h>
+
+extern void* allocate(int alignment, size_t bytesize);
+
+#endif

PoissonSolver/2D-omp/src/likwid-marker.h

@@ -0,0 +1,53 @@
+/*
+ * =======================================================================================
+ *
+ *      Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *      Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *
+ *      Permission is hereby granted, free of charge, to any person obtaining a copy
+ *      of this software and associated documentation files (the "Software"), to deal
+ *      in the Software without restriction, including without limitation the rights
+ *      to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ *      copies of the Software, and to permit persons to whom the Software is
+ *      furnished to do so, subject to the following conditions:
+ *
+ *      The above copyright notice and this permission notice shall be included in all
+ *      copies or substantial portions of the Software.
+ *
+ *      THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ *      IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ *      FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ *      AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ *      LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ *      OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ *      SOFTWARE.
+ *
+ * =======================================================================================
+ */
+#ifndef LIKWID_MARKERS_H
+#define LIKWID_MARKERS_H
+
+#ifdef LIKWID_PERFMON
+#include <likwid.h>
+#define LIKWID_MARKER_INIT likwid_markerInit()
+#define LIKWID_MARKER_THREADINIT likwid_markerThreadInit()
+#define LIKWID_MARKER_SWITCH likwid_markerNextGroup()
+#define LIKWID_MARKER_REGISTER(regionTag) likwid_markerRegisterRegion(regionTag)
+#define LIKWID_MARKER_START(regionTag) likwid_markerStartRegion(regionTag)
+#define LIKWID_MARKER_STOP(regionTag) likwid_markerStopRegion(regionTag)
+#define LIKWID_MARKER_CLOSE likwid_markerClose()
+#define LIKWID_MARKER_RESET(regionTag) likwid_markerResetRegion(regionTag)
+#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count) likwid_markerGetRegion(regionTag, nevents, events, time, count)
+#else  /* LIKWID_PERFMON */
+#define LIKWID_MARKER_INIT
+#define LIKWID_MARKER_THREADINIT
+#define LIKWID_MARKER_SWITCH
+#define LIKWID_MARKER_REGISTER(regionTag)
+#define LIKWID_MARKER_START(regionTag)
+#define LIKWID_MARKER_STOP(regionTag)
+#define LIKWID_MARKER_CLOSE
+#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)
+#define LIKWID_MARKER_RESET(regionTag)
+#endif /* LIKWID_PERFMON */
+
+#endif /*LIKWID_MARKERS_H*/

PoissonSolver/2D-omp/src/main.c

@@ -0,0 +1,76 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.ke
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "likwid-marker.h"
+#include "parameter.h"
+#include "solver.h"
+#include "timing.h"
+#include "omp.h"
+
+#define LIKWID_PROFILE(tag, call)                                                        \
+    startTime = getTimeStamp();                                                          \
+    LIKWID_MARKER_START(#tag);                                                           \
+    call(&solver);                                                                       \
+    LIKWID_MARKER_STOP(#tag);                                                            \
+    endTime = getTimeStamp();
+
+enum VARIANT { SOR = 1, RB, RBA };
+
+int main(int argc, char** argv)
+{
+    int volatile dummy = 0;
+    int variant = RB;
+    double startTime, endTime;
+    Parameter params;
+    Solver solver;
+    initParameter(&params);
+    LIKWID_MARKER_INIT;
+    #pragma omp parallel
+    {
+        if(dummy==1 || omp_get_thread_num()==0)
+        printf("OMP_THREADS_DETECTED: %d\n",omp_get_num_threads());
+    }
+    if (argc < 2) {
+        printf("Usage: %s <configFile>\n", argv[0]);
+        exit(EXIT_SUCCESS);
+    }
+
+    readParameter(&params, argv[1]);
+    // printParameter(&params);
+    if (argc == 3) {
+        variant = atoi(argv[2]);
+    }
+    if (argc == 4) {
+        sscanf("%lf", argv[3], &params.omg);
+    }
+
+    initSolver(&solver, &params, 2);
+    writeResult(&solver, "p-0.dat");
+
+    switch (variant) {
+    case SOR:
+        printf("Plain SOR\n");
+        fflush(stdout);
+        LIKWID_PROFILE("SOR", solve);
+        break;
+    case RB:
+        printf("Red-black SOR\n");
+        fflush(stdout);
+        LIKWID_PROFILE("RB", solveRB);
+        break;
+    case RBA:
+        printf("Red-black SOR with acceleration\n");
+        fflush(stdout);
+        LIKWID_PROFILE("RBA", solveRBA);
+        break;
+    }
+    printf(" %.2fs\n", endTime - startTime);
+    writeResult(&solver, "p-final.dat");
+
+    LIKWID_MARKER_CLOSE;
+    return EXIT_SUCCESS;
+}

PoissonSolver/2D-omp/src/parameter.c

@@ -0,0 +1,79 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+//---
+#include "parameter.h"
+#include "util.h"
+#define MAXLINE 4096
+
+void initParameter(Parameter* param)
+{
+    param->xlength = 1.0;
+    param->ylength = 1.0;
+    param->imax    = 100;
+    param->jmax    = 100;
+    param->itermax = 1000;
+    param->eps     = 0.0001;
+    param->omg     = 1.8;
+    param->rho     = 0.99;
+}
+
+void readParameter(Parameter* param, const char* filename)
+{
+    FILE* fp = fopen(filename, "r");
+    char line[MAXLINE];
+    int i;
+
+    if (!fp) {
+        fprintf(stderr, "Could not open parameter file: %s\n", filename);
+        exit(EXIT_FAILURE);
+    }
+
+    while (!feof(fp)) {
+        line[0] = '\0';
+        fgets(line, MAXLINE, fp);
+        for (i = 0; line[i] != '\0' && line[i] != '#'; i++)
+            ;
+        line[i] = '\0';
+
+        char* tok = strtok(line, " ");
+        char* val = strtok(NULL, " ");
+
+#define PARSE_PARAM(p, f)                                                                \
+    if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) {                         \
+        param->p = f(val);                                                               \
+    }
+#define PARSE_STRING(p) PARSE_PARAM(p, strdup)
+#define PARSE_INT(p)    PARSE_PARAM(p, atoi)
+#define PARSE_REAL(p)   PARSE_PARAM(p, atof)
+
+        if (tok != NULL && val != NULL) {
+            PARSE_REAL(xlength);
+            PARSE_REAL(ylength);
+            PARSE_INT(imax);
+            PARSE_INT(jmax);
+            PARSE_INT(itermax);
+            PARSE_REAL(eps);
+            PARSE_REAL(omg);
+            PARSE_REAL(rho);
+        }
+    }
+
+    fclose(fp);
+}
+
+void printParameter(Parameter* param)
+{
+    printf("Parameters:\n");
+    printf("Geometry data:\n");
+    printf("\tDomain box size (x, y): %e, %e\n", param->xlength, param->ylength);
+    printf("\tCells (x, y): %d, %d\n", param->imax, param->jmax);
+    printf("Iterative solver parameters:\n");
+    printf("\tMax iterations: %d\n", param->itermax);
+    printf("\tepsilon (stopping tolerance) : %e\n", param->eps);
+    printf("\tomega (SOR relaxation): %e\n", param->omg);
+}

PoissonSolver/2D-omp/src/parameter.h

@@ -0,0 +1,18 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#ifndef __PARAMETER_H_
+#define __PARAMETER_H_
+
+typedef struct {
+    double xlength, ylength;
+    int imax, jmax;
+    int itermax;
+    double eps, omg, rho, gamma;
+} Parameter;
+
+void initParameter(Parameter*);
+void readParameter(Parameter*, const char*);
+void printParameter(Parameter*);
+#endif

PoissonSolver/2D-omp/src/solver.c

@@ -0,0 +1,276 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+
+#include "allocate.h"
+#include "parameter.h"
+#include "solver.h"
+
+#define PI        3.14159265358979323846
+#define P(i, j)   p[(j) * (imax + 2) + (i)]
+#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
+
+void initSolver(Solver* solver, Parameter* params, int problem)
+{
+    solver->imax    = params->imax;
+    solver->jmax    = params->jmax;
+    solver->dx      = params->xlength / params->imax;
+    solver->dy      = params->ylength / params->jmax;
+    solver->eps     = params->eps;
+    solver->omega   = params->omg;
+    solver->rho     = params->rho;
+    solver->itermax = params->itermax;
+
+    int imax        = solver->imax;
+    int jmax        = solver->jmax;
+    size_t bytesize = (imax + 2) * (jmax + 2) * sizeof(double);
+    solver->p       = allocate(64, bytesize);
+    solver->rhs     = allocate(64, bytesize);
+
+    double dx   = solver->dx;
+    double dy   = solver->dy;
+    double* p   = solver->p;
+    double* rhs = solver->rhs;
+    #pragma omp parallel for collapse(2)
+    for (int j = 0; j < jmax + 2; j++) {
+        for (int i = 0; i < imax + 2; i++) {
+            P(i, j) = sin(2.0 * PI * i * dx * 2.0) + sin(2.0 * PI * j * dy * 2.0);
+        }
+    }
+
+    if (problem == 2) {
+        #pragma omp parallel for collapse(2)
+        for (int j = 0; j < jmax + 2; j++) {
+            for (int i = 0; i < imax + 2; i++) {
+                RHS(i, j) = sin(2.0 * PI * i * dx);
+            }
+        }
+    } else {
+        #pragma omp parallel for collapse(2)
+        for (int j = 0; j < jmax + 2; j++) {
+            for (int i = 0; i < imax + 2; i++) {
+                RHS(i, j) = 0.0;
+            }
+        }
+    }
+}
+
+void solve(Solver* solver)
+{
+    int imax      = solver->imax;
+    int jmax      = solver->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->dx * solver->dx;
+    double dy2    = solver->dy * solver->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double* p     = solver->p;
+    double* rhs   = solver->rhs;
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+    char filename[20];
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+
+                double r = RHS(i, j) -
+                           ((P(i - 1, j) - 2.0 * P(i, j) + P(i + 1, j)) * idx2 +
+                               (P(i, j - 1) - 2.0 * P(i, j) + P(i, j + 1)) * idy2);
+
+                P(i, j) -= (factor * r);
+                res += (r * r);
+            }
+        }
+
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0)        = P(i, 1);
+            P(i, jmax + 1) = P(i, jmax);
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j)        = P(1, j);
+            P(imax + 1, j) = P(imax, j);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+#ifdef ANIMATE
+        sprintf(filename, "p-%d.dat", it);
+        writeResult(solver, filename);
+#endif
+        it++;
+    }
+
+    printf("%d, %f\n", it, solver->omega);
+}
+
+void solveRB(Solver* solver)
+{
+    int imax      = solver->imax;
+    int jmax      = solver->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->dx * solver->dx;
+    double dy2    = solver->dy * solver->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double* p     = solver->p;
+    double* rhs   = solver->rhs;
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+    int pass, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        jsw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+            #pragma omp parallel for firstprivate(isw)
+            for (int j = 1; j < jmax + 1; j++) {
+                for (int i = isw; i < imax + 1; i += 2) {
+
+                    double r = RHS(i, j) -
+                               ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
+                                   (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
+
+                    P(i, j) -= (factor * r);
+                    res += (r * r);
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+        #pragma omp parallel for 
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0)        = P(i, 1);
+            P(i, jmax + 1) = P(i, jmax);
+        }
+        #pragma omp parallel for
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j)        = P(1, j);
+            P(imax + 1, j) = P(imax, j);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+// #ifdef ANIMATE
+//         sprintf(filename, "p-%d.dat", it);
+//         writeResult(solver, filename);
+// #endif
+//         it++;
+    }
+
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+    printf("%d, %f\n", it, solver->omega);
+}
+
+void solveRBA(Solver* solver)
+{
+    int imax      = solver->imax;
+    int jmax      = solver->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->dx * solver->dx;
+    double dy2    = solver->dy * solver->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double rho    = solver->rho;
+    double* p     = solver->p;
+    double* rhs   = solver->rhs;
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+    int pass, jsw, isw;
+    double omega = 1.0;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        jsw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+
+            for (int j = 1; j < jmax + 1; j++) {
+                for (int i = isw; i < imax + 1; i += 2) {
+
+                    double r = RHS(i, j) -
+                               ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
+                                   (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
+
+                    P(i, j) -= (omega * factor * r);
+                    res += (r * r);
+                }
+                isw = 3 - isw;
+            }
+            jsw   = 3 - jsw;
+            omega = (it == 0 && pass == 0 ? 1.0 / (1.0 - 0.5 * rho * rho)
+                                          : 1.0 / (1.0 - 0.25 * rho * rho * omega));
+        }
+
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0)        = P(i, 1);
+            P(i, jmax + 1) = P(i, jmax);
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j)        = P(1, j);
+            P(imax + 1, j) = P(imax, j);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e Omega: %e\n", it, res, omega);
+#endif
+#ifdef ANIMATE
+        sprintf(filename, "p-%d.dat", it);
+        writeResult(solver, filename);
+#endif
+        it++;
+    }
+
+    // printf("Final omega: %f\n", omega);
+    // printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+    printf("%d, %f\n", it, omega);
+}
+
+void writeResult(Solver* solver, char* filename)
+{
+    int imax  = solver->imax;
+    int jmax  = solver->jmax;
+    double* p = solver->p;
+
+    FILE* fp;
+    fp = fopen(filename, "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    for (int j = 0; j < jmax + 2; j++) {
+        for (int i = 0; i < imax + 2; i++) {
+            fprintf(fp, "%f ", P(i, j));
+        }
+        fprintf(fp, "\n");
+    }
+
+    fclose(fp);
+}

PoissonSolver/2D-omp/src/solver.h

@@ -0,0 +1,22 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#ifndef __SOLVER_H_
+#define __SOLVER_H_
+#include "parameter.h"
+
+typedef struct {
+    double dx, dy;
+    int imax, jmax;
+    double *p, *rhs;
+    double eps, omega, rho;
+    int itermax;
+} Solver;
+
+extern void initSolver(Solver*, Parameter*, int problem);
+extern void writeResult(Solver*, char*);
+extern void solve(Solver*);
+extern void solveRB(Solver*);
+extern void solveRBA(Solver*);
+#endif

PoissonSolver/2D-omp/src/timing.c

@@ -0,0 +1,22 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#include <stdlib.h>
+#include <time.h>
+
+double getTimeStamp()
+{
+    struct timespec ts;
+    clock_gettime(CLOCK_MONOTONIC, &ts);
+    return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
+}
+
+double getTimeResolution()
+{
+    struct timespec ts;
+    clock_getres(CLOCK_MONOTONIC, &ts);
+    return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
+}
+
+double getTimeStamp_() { return getTimeStamp(); }

PoissonSolver/2D-omp/src/timing.h

@@ -0,0 +1,11 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#ifndef __TIMING_H_
+#define __TIMING_H_
+
+extern double getTimeStamp();
+extern double getTimeResolution();
+
+#endif // __TIMING_H_

PoissonSolver/2D-omp/src/util.h

@@ -0,0 +1,20 @@
+/* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.    */
+#ifndef __UTIL_H_
+#define __UTIL_H_
+#define HLINE                                                                            \
+    "----------------------------------------------------------------------------\n"
+
+#ifndef MIN
+#define MIN(x, y) ((x) < (y) ? (x) : (y))
+#endif
+#ifndef MAX
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef ABS
+#define ABS(a) ((a) >= 0 ? (a) : -(a))
+#endif
+
+#endif // __UTIL_H_

PoissonSolver/2D-omp/surface.plot

@@ -0,0 +1,7 @@
+set terminal png size 1024,768 enhanced font ,12
+set output 'p.png'
+set datafile separator whitespace
+
+set grid
+set hidden3d
+splot 'p.dat' matrix using 1:2:3 with lines