sheet 5

Sheet_5/bsp_5_3/Makefile

@@ -0,0 +1,30 @@
+#
+# use GNU-Compiler tools
+COMPILER=GCC_
+# alternatively from the shell
+# export COMPILER=GCC_
+# or, alternatively from the shell
+# make COMPILER=GCC_
+
+# use Intel compilers
+#COMPILER=ICC_
+
+# use PGI compilers
+# COMPILER=PGI_
+
+
+SOURCES = main.cpp bsp_5_3_lib.cpp
+OBJECTS = $(SOURCES:.cpp=.o)
+
+PROGRAM	= main.${COMPILER}
+
+# uncomment the next to lines for debugging and detailed performance analysis
+CXXFLAGS += -g
+LINKFLAGS += -g
+# do not use -pg with PGI compilers
+
+ifndef COMPILER
+  COMPILER=GCC_
+endif
+
+include ../${COMPILER}default.mk

Sheet_5/bsp_5_3/bsp_5_3_lib.cpp

@@ -0,0 +1,101 @@
+#include "bsp_5_3_lib.h"
+#include "mayer_primes.h"
+#include "algorithm"
+#include <iostream>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+using namespace std;
+
+int single_goldbach(const int &k)
+{
+    vector<int> primes = get_primes(k);
+    int amount = 0;
+    for(size_t it = 0; primes[it]<=k/2.0; ++it) //für Primzahl größer als k/2 haben wir bereits Zerlegung gezählt: 3+7 = 7+3
+    {
+       for(size_t j = it; j<primes.size(); ++j)
+       {
+           if(primes[it] + primes[j] == k)
+           {
+               amount += 1;
+           }
+       }
+    }
+
+    return amount;
+}
+
+
+int single_goldbach_par(const int &k)
+{
+    vector<int> primes = get_primes(k);
+    int amount = 0;
+    size_t it_max = 0;
+    // to get a barrier for the loop; omp needs a fix upper bound
+	while (it_max < primes.size() && primes[it_max] <= k/2)
+    {
+		++it_max;
+	}
+    
+#pragma omp parallel for default(none) shared(primes, k, it_max) reduction(+:amount)
+    for(size_t it = 0; it<it_max; ++it) //für Primzahl größer als k/2 haben wir bereits Zerlegung gezählt: 3+7 = 7+3
+    {
+       for(size_t j = it; j<primes.size(); ++j)
+       {
+           if(primes[it] + primes[j] == k)
+           {
+               amount += 1;
+           }
+       }
+    }
+
+    return amount;
+}
+
+
+vector<int> count_goldbach(const int &n)
+{
+    vector<int> count_vec((n-4)/2+1,0);
+    vector<int> primes = get_primes(n);
+    for(size_t k=0; k < primes.size() && primes[k]<=n/2; ++k)
+    {
+        for(size_t i=k; i<primes.size(); ++i)
+        {
+            int sum = primes[k] + primes[i];
+            if(sum<=n)
+            {
+                count_vec[(sum-4)/2] += 1;
+            }
+        }
+    }
+
+    return count_vec;
+}
+
+
+vector<int> count_goldbach_par(const int &n)
+{
+    vector<int> count_vec((n-4)/2+1,0);
+    vector<int> primes = get_primes(n);
+    size_t it_max = 0;
+    while (it_max < primes.size() && primes[it_max] <= n/2)
+    {
+		++it_max;
+	}
+    
+#pragma omp parallel for schedule(dynamic) shared(primes, n, it_max) reduction(VecAdd:count_vec)
+    for(size_t k=0; k<it_max; ++k)
+    {
+        for(size_t i=k; i<primes.size(); ++i)
+        {
+            int sum = primes[k] + primes[i];
+            if(sum<=n)
+            {
+                count_vec[(sum-4)/2] += 1;
+            }
+        }
+    }
+
+    return count_vec;
+}

Sheet_5/bsp_5_3/bsp_5_3_lib.h

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <cassert>
+#include <omp.h>
+#include <vector>
+
+/** \brief Zaehlt fuer eine gegebene Zahl @k die Anzahl der moeglichen Zerlegungen als Summe zweier Primzahlen
+ *
+ * \param[in]   k   fuer diese Zahl wird die Anzahl der Zerlegungen berechnet
+ * \return  Anzahl der Zerlegungen
+ *
+ */
+int single_goldbach(const int &k);
+
+
+/** \brief Zaehlt fuer eine gegebene Zahl @k die Anzahl der moeglichen Zerlegungen als Summe zweier Primzahlen [parallel]
+ *
+ * \param[in]   k   fuer diese Zahl wird die Anzahl der Zerlegungen berechnet
+ * \return  Anzahl der Zerlegungen
+ *
+ */
+int single_goldbach_par(const int &k);
+
+
+/** \brief Zaehlt die Anzahl der Dekomposition für alle geraden Zahlen in [4,n]
+ *
+ * \param[in]   n   obere Intervallgrenze (in diesem Bereich werden die Dekompositionen berechnet
+ * \return Vektor mit den Anzahl der Dekompositionen (ad Adressierung: x[0] = (n=4), x[1] = (n=6), x[2] = (n=8), Umrechnung *2 + 4
+ *
+ */
+std::vector<int> count_goldbach(const int &n);
+
+
+/** \brief Zaehlt die Anzahl der Dekomposition für alle geraden Zahlen in [4,n] - [parallel]
+ *
+ * \param[in]   n   obere Intervallgrenze (in diesem Bereich werden die Dekompositionen berechnet
+ * \return Vektor mit den Anzahl der Dekompositionen (ad Adressierung: x[0] = (n=4), x[1] = (n=6), x[2] = (n=8), Umrechnung *2 + 4
+ *
+ */
+std::vector<int> count_goldbach_par(const int &n);
+
+
+/** 	 Vector @p b adds its elements to vector @p a .
+	@param[in] a	vector
+	@param[in] b	vector
+	@return 	     a+=b componentwise
+*/
+template<class T>
+std::vector<T> &operator+=(std::vector<T> &a, std::vector<T> const &b)
+{
+    assert(a.size()==b.size());
+    for (size_t k = 0; k < a.size(); ++k) {
+        a[k] += b[k];
+    }
+    return a;
+}
+
+#pragma omp declare reduction(VecAdd : std::vector<int>  : omp_out += omp_in) \
+ initializer (omp_priv=omp_orig)

Sheet_5/bsp_5_3/main.cpp

@@ -0,0 +1,46 @@
+#include "bsp_5_3_lib.h"
+#include "mayer_primes.h"
+#include <iostream>
+#include <algorithm>
+#include <chrono>
+// BSP 5_3 Goldbach conjunction
+using namespace std;
+
+int main()
+{
+	omp_set_num_threads(8);
+	
+	cout << "\nChecking for correct result\n";
+	cout << "694 has " << single_goldbach_par(694) << " decompositions" << endl;
+	
+    // Auswertung für n=100000 bzw herausfinden der Zahl, welche die meisten Dekompositionen hat
+    vector<int> v = count_goldbach_par(1e5);
+    auto ip = max_element(v.begin(), v.end());
+    cout << "number in [4,1e5] with most decompositions: " << distance(v.begin(), ip)*2+4 << " has " << *ip << " decompositions" << endl;
+
+    cout << "\nTiming for parallel" << endl;
+    vector<int> nvec{static_cast<int>(1e4), static_cast<int>(1e5), static_cast<int>(4*1e5), static_cast<int>(1e6), static_cast<int>(2*1e6)}; // Vektor für n
+    for(size_t k=0; k<nvec.size(); ++k)
+    {
+        auto timestart = omp_get_wtime();
+        vector<int> vall = count_goldbach_par(nvec[k]);
+        auto time = omp_get_wtime() - timestart;
+        cout << "n = " << nvec[k] << "\ttime in s: " << time << endl;
+    }
+    
+    cout << "\nTiming for serial" << endl;
+    for(size_t k=0; k<nvec.size(); ++k)
+    {
+        auto timestart = omp_get_wtime();
+        vector<int> vall = count_goldbach(nvec[k]);
+        auto time = omp_get_wtime() - timestart;
+        cout << "n = " << nvec[k] << "\ttime in s: " << time << endl;
+    }
+    
+    // for large n the parallel version is slower than the sequential
+    // the reason is the reduction(VecAdd) which is slow especially for large vectors
+
+//    cout << single_goldbach(694) << endl;
+//    cout << count_goldbach(10000)[690/2] << endl;
+    return 0;
+}

Sheet_5/bsp_5_3/mayer_primes.h

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <cstring> //memset
+#include <vector>
+//using namespace std;
+
+/** \brief Determines all prime numbers in interval [2, @p max].
+ *
+ *  The sieve of Eratosthenes is used.
+ *
+ *  The implementation originates from <a href="http://code.activestate.com/recipes/576559-fast-prime-generator/">Florian Mayer</a>.
+ *
+ * \param[in]   max end of interval for the prime number search.
+ * \return  vector of prime numbers @f$2,3,5, ..., p<=max @f$.
+ *
+ * \copyright
+ * Copyright (c) 2008 Florian Mayer          (adapted by Gundolf Haase 2018)
+ *
+ * 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.
+ *
+ */
+template <class T>
+std::vector<T> get_primes(T max)
+{
+    std::vector<T> primes;
+    char *sieve;
+    sieve = new char[max / 8 + 1];
+    // Fill sieve with 1
+    memset(sieve, 0xFF, (max / 8 + 1) * sizeof(char));
+    for (T x = 2; x <= max; x++)
+    {
+        if (sieve[x / 8] & (0x01 << (x % 8))) {
+            primes.push_back(x);
+            // Is prime. Mark multiplicates.
+            for (T j = 2 * x; j <= max; j += x)
+            {
+                sieve[j / 8] &= ~(0x01 << (j % 8));
+			}
+        }
+	}
+    delete[] sieve;
+    return primes;
+}
+
+//---------------------------------------------------------------
+//int main()        // by  Florian Mayer
+//{g++ -O3  -std=c++14 -fopenmp main.cpp && ./a.out
+//    vector<unsigned long> primes;
+//    primes = get_primes(10000000);
+//    // return 0;
+//    // Print out result.
+//    vector<unsigned long>::iterator it;
+//    for(it=primes.begin(); it < primes.end(); it++)
+//        cout << *it << " ";
+//
+//    cout << endl;
+//    return 0;
+//}