#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <unistd.h>
#ifdef _OPENMP
#include <omp.h>
#else
void omp_set_num_threads(int num) {
return;
}
double omp_get_wtime() {
struct timeval T11;
gettimeofday(&T11, NULL);
return (T11.tv_sec) * 1000 + (T11.tv_usec) / 1000;
}
#endif
#define A 819
#define DBL_MAX __DBL_MAX__
void sort (double *M2, int from, int to) {
for (int i = from; i < to - 1; i++) {
int indexMin = i;
for (int j = i + 1; j < to; j++) {
if (M2[indexMin] > M2[j])
indexMin = j;
}
if (indexMin != i) {
double swap = M2[i];
M2[i] = M2[indexMin];
M2[indexMin] = swap;
}
}
}
void ex(int argc, char *argv[], int *stage) {
const int N = atoi(argv[1]);
const int N2 = N / 2;
//const int N2m1 = N2 - 1;
const int randIndM1 = RAND_MAX / A;
const int randIndM2 = RAND_MAX / (9 * A);
const int CountTest = 50;
double timeStart, timeEnd;
double delta_ms;
timeStart = omp_get_wtime();
omp_set_num_threads(4);
double X = 0.0;
for (unsigned int numTest = 1; numTest <= CountTest; ++numTest) {
double min;
double M1[N];
double M2[N2];
double M2copy[N2];
double randM1[N];
double randM2[N2];
unsigned int myseed;
myseed = numTest;
for (int i = 0; i < N; i++) {
randM1[i] = rand_r(&myseed);
}
for (int i = 0; i < N2; i++) {
randM2[i] = rand_r(&myseed);
}
#pragma omp parallel for shared(M1)
for (int i = 0; i < N; i++) {
M1[i] = randM1[i] / randIndM1;
}
#pragma omp parallel for shared(M2)
for (int i = 0; i < N2; i++) {
M2[i] = A + (randM2[i] / randIndM2);
}
*stage = 20;
// Stage 2#7 begin
#pragma omp parallel for shared(M1)
for (int i = 0; i < N; i++) {
M1[i] = exp(sqrt(M1[i]));
}
for (int i = 0; i < N2; i++) {
double cur = M2[i];
double prev = cur;
if (i == 0) {
prev = 0.0;
} else {
prev = M2[i - 1];
}
M2[i] = log(fabs(tan(prev + cur)));
}
*stage = 40;
// Stage 3 (merging) begin
#pragma omp parallel for default(none) shared(M1, M2)
for (int i = 0; i < N2; i++) {
M2[i] *= M1[i];
}
*stage = 60;
//Stage 4 (sorting) begin
#pragma omp parallel sections
{
#pragma omp section
{
sort(M2, 0, N2 / 2);
}
#pragma omp section
{
sort(M2, N2 / 2, N2);
}
}
for (int i = 0, j = N2 / 2, p = 0; p < N2; p++) {
if ((j == N2) || (i < N2 / 2 && M2[i] < M2[j])) {
M2copy[p] = M2[i++];
} else {
M2copy[p] = M2[j++];
}
}
*stage = 80;
//Stage 5 (getting X)
min = DBL_MAX;
#pragma omp parallel for default(none) shared(M2) reduction(min:min)
for (int i = 0; i < N2; i++) {
double current = M2[i];
if (current < min && current > 0) {
min = current;
}
}
#pragma omp parallel for default(none) shared(M2, min) reduction(+:X)
for (int i = 0; i < N2; i++) {
if ((int) (M2[i] / min) % 2 == 0) {
X += sin(M2[i]);
}
}
*stage = 99;
}
timeEnd = omp_get_wtime();
delta_ms = 1000 * (timeEnd - timeStart);
printf("\nN=%d. Milliseconds passed: %f\n", N, delta_ms);
printf("\nX=%f. \n", X);
*stage = 100;
}
int main(int argc, char *argv[]) {
int stage = 0;
omp_set_max_active_levels(2);
#pragma omp parallel sections
{
#pragma omp section
{
while (1) {
printf("\r%d%%", stage);
fflush(stdout);
sleep(1);
if (stage == 100) {
break;
}
}
}
#pragma omp section
{
ex(argc, argv, &stage);
}
}
return 0;
}