diff --git a/Sheet_5/bsp_5_3/Makefile b/Sheet_5/bsp_5_3/Makefile new file mode 100644 index 0000000..5b84561 --- /dev/null +++ b/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 diff --git a/Sheet_5/bsp_5_3/bsp_5_3_lib.cpp b/Sheet_5/bsp_5_3/bsp_5_3_lib.cpp new file mode 100644 index 0000000..4e19ef7 --- /dev/null +++ b/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 +#ifdef _OPENMP +#include +#endif + +using namespace std; + +int single_goldbach(const int &k) +{ + vector 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 = 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 count_goldbach(const int &n) +{ + vector count_vec((n-4)/2+1,0); + vector primes = get_primes(n); + for(size_t k=0; k < primes.size() && primes[k]<=n/2; ++k) + { + for(size_t i=k; i count_goldbach_par(const int &n) +{ + vector count_vec((n-4)/2+1,0); + vector 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 +#include +#include + +/** \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 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 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 +std::vector &operator+=(std::vector &a, std::vector 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 : omp_out += omp_in) \ + initializer (omp_priv=omp_orig) diff --git a/Sheet_5/bsp_5_3/main.cpp b/Sheet_5/bsp_5_3/main.cpp new file mode 100644 index 0000000..ab46935 --- /dev/null +++ b/Sheet_5/bsp_5_3/main.cpp @@ -0,0 +1,46 @@ +#include "bsp_5_3_lib.h" +#include "mayer_primes.h" +#include +#include +#include +// 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 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 nvec{static_cast(1e4), static_cast(1e5), static_cast(4*1e5), static_cast(1e6), static_cast(2*1e6)}; // Vektor für n + for(size_t k=0; k 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 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; +} diff --git a/Sheet_5/bsp_5_3/mayer_primes.h b/Sheet_5/bsp_5_3/mayer_primes.h new file mode 100644 index 0000000..6c6dc23 --- /dev/null +++ b/Sheet_5/bsp_5_3/mayer_primes.h @@ -0,0 +1,73 @@ +#pragma once + +#include //memset +#include +//using namespace std; + +/** \brief Determines all prime numbers in interval [2, @p max]. + * + * The sieve of Eratosthenes is used. + * + * The implementation originates from Florian Mayer. + * + * \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 +std::vector get_primes(T max) +{ + std::vector 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 primes; +// primes = get_primes(10000000); +// // return 0; +// // Print out result. +// vector::iterator it; +// for(it=primes.begin(); it < primes.end(); it++) +// cout << *it << " "; +// +// cout << endl; +// return 0; +//}