140 lines
3.5 KiB
C++
140 lines
3.5 KiB
C++
#include "mylib.h"
|
|
#include <cassert> // assert()
|
|
#include <cmath>
|
|
#include <iostream>
|
|
#include <functional> // multiplies<>{}
|
|
#include <list>
|
|
#include <numeric> // iota()
|
|
#ifdef _OPENMP
|
|
#include <omp.h>
|
|
#endif
|
|
#include <vector>
|
|
using namespace std;
|
|
|
|
double scalar(vector<double> const &x, vector<double> const &y)
|
|
{
|
|
assert(x.size() == y.size()); // switch off via compile flag: -DNDEBUG
|
|
size_t const N = x.size();
|
|
double sum = 0.0;
|
|
// int check; // for using omp_in_parallel
|
|
//#pragma omp parallel for default(none) shared(x,y,N,check) reduction(+:sum)
|
|
#pragma omp parallel for default(none) shared(x,y,N) reduction(+:sum)
|
|
for (size_t i = 0; i < N; ++i)
|
|
{
|
|
sum += x[i] * y[i];
|
|
// check = omp_in_parallel();
|
|
//sum += exp(x[i])*log(y[i]);
|
|
}
|
|
// cout << "in parallel: " << check << endl;
|
|
return sum;
|
|
}
|
|
|
|
|
|
double scalar_manual(vector<double> const &x, vector<double> const &y)
|
|
{
|
|
assert(x.size() == y.size());
|
|
size_t const N = x.size();
|
|
double sum = 0.0;
|
|
|
|
#pragma omp parallel default(none) shared(x,y,N,sum)
|
|
{
|
|
int tid = omp_get_thread_num();
|
|
int nt = omp_get_num_threads();
|
|
|
|
// manual splitting of the index area
|
|
size_t start = (N * tid) / nt;
|
|
size_t end = (N * (tid+1)) / nt;
|
|
|
|
double local_sum = 0.0;
|
|
for (size_t i = start; i < end; ++i) {
|
|
local_sum += x[i] * y[i];
|
|
}
|
|
|
|
// local subtotal combined to global sum
|
|
#pragma omp atomic
|
|
sum += local_sum;
|
|
}
|
|
|
|
return sum;
|
|
}
|
|
|
|
|
|
double norm(vector<double> const &x)
|
|
{
|
|
size_t const N = x.size();
|
|
double sum = 0.0;
|
|
#pragma omp parallel for default(none) shared(x,N) reduction(+:sum)
|
|
for (size_t i = 0; i < N; ++i)
|
|
{
|
|
sum += x[i]*x[i];
|
|
}
|
|
return sum;
|
|
}
|
|
|
|
|
|
|
|
vector<int> reduction_vec(int n)
|
|
{
|
|
vector<int> vec(n);
|
|
#pragma omp parallel default(none) shared(cout) reduction(VecAdd:vec)
|
|
{
|
|
#pragma omp barrier
|
|
#pragma omp critical
|
|
cout << omp_get_thread_num() << " : " << vec.size() << endl;
|
|
#pragma omp barrier
|
|
iota( vec.begin(),vec.end(), omp_get_thread_num() );
|
|
#pragma omp barrier
|
|
|
|
}
|
|
return vec;
|
|
}
|
|
|
|
|
|
vector<int> reduction_vec_append(int n)
|
|
{
|
|
vector<int> vec;
|
|
|
|
#pragma omp parallel default(none) shared(n,cout) reduction(VecAppend:vec)
|
|
{
|
|
std::vector<int> local(n); // local vector of each thread
|
|
|
|
// consecutive numbers starting from thread ID
|
|
iota(local.begin(), local.end(), omp_get_thread_num());
|
|
|
|
// output for checking
|
|
#pragma omp critical
|
|
cout << "Thread " << omp_get_thread_num() << " local size = " << local.size() << std::endl;
|
|
|
|
// append local to vec
|
|
vec.insert(vec.end(), local.begin(), local.end());
|
|
}
|
|
|
|
return vec;
|
|
}
|
|
|
|
|
|
|
|
double scalar_trans(vector<double> const &x, vector<double> const &y)
|
|
{
|
|
assert(x.size() == y.size()); // switch off via compile flag: -DNDEBUG
|
|
vector<double> z(x.size());
|
|
//list<double> z(x.size()); // parallel for-loop on iterators not possible (missing 'operator-')
|
|
// c++-20 CLANG_, ONEAPI_:condition of OpenMP for loop must be a relational comparison
|
|
|
|
transform(cbegin(x),cend(x),cbegin(y),begin(z),std::multiplies<>{});
|
|
|
|
double sum = 0.0;
|
|
#pragma omp parallel for default(none) shared(z) reduction(+:sum)
|
|
for (auto pi = cbegin(z); pi!=cend(z); ++pi)
|
|
{
|
|
sum += *pi;
|
|
}
|
|
//for (auto val: z)
|
|
//{
|
|
//sum += val;
|
|
//}
|
|
return sum;
|
|
}
|
|
|
|
|
|
|