75 lines
2.1 KiB
C
75 lines
2.1 KiB
C
// metoda za adaptivno numericno integracijo
|
|
// prevedi s: gcc -O2 quad.c -fopenmp -lm -o quad
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <omp.h>
|
|
|
|
#define TOL 1e-8
|
|
|
|
double func(double x)
|
|
{
|
|
return sin(x*x);
|
|
}
|
|
|
|
double quad(double (*f)(double), double lower, double upper, double tol)
|
|
// funkcija za integracijo, spodnja meja, zgornja meja, dovoljena napaka
|
|
{
|
|
double quad_res; // rezultat
|
|
double h; // dolzina intervala
|
|
double middle; // sredina intervala
|
|
double quad_coarse; // groba aproksimacija
|
|
double quad_fine; // fina aproksimacija (two trapezoids)
|
|
double quad_lower; // rezultat na spodnjem intervalu
|
|
double quad_upper; // rezultat na zgornjem intervalu
|
|
double eps; // razlika
|
|
|
|
h = upper - lower;
|
|
middle = (lower + upper) / 2;
|
|
|
|
// izracunaj integral z obema aproksimacijama trapezoidnega pravila
|
|
quad_coarse = h * (f(lower) + f(upper)) / 2.0;
|
|
quad_fine = h/2 * (f(lower) + f(middle)) / 2.0 + h/2 * (f(middle) + f(upper)) / 2.0;
|
|
eps = fabs(quad_coarse - quad_fine);
|
|
|
|
//ce se nismo dosegli zelene natancnosti, razdelimo interval na pol in ponovimo
|
|
if (eps > tol)
|
|
{
|
|
#pragma omp parallel
|
|
{
|
|
#pragma omp sections
|
|
{
|
|
#pragma omp section
|
|
{
|
|
# pragma omp task firstprivate(lower, middle, tol)
|
|
quad_lower = quad(f, lower, middle, tol/2);
|
|
}
|
|
#pragma omp section
|
|
{
|
|
# pragma omp task firstprivate(middle, upper, tol)
|
|
quad_upper = quad(f, middle, upper, tol/2);
|
|
}
|
|
}
|
|
quad_res = quad_lower + quad_upper;
|
|
}
|
|
}
|
|
else
|
|
quad_res = quad_fine;
|
|
|
|
return quad_res;
|
|
}
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
double quadrature;
|
|
double dt = omp_get_wtime();
|
|
|
|
quadrature = quad(func, 0.0, 50.0, TOL);
|
|
|
|
dt = omp_get_wtime() - dt;
|
|
|
|
printf("Integral: %lf\nCas: %lf\n", quadrature, dt);
|
|
|
|
return 0;
|
|
}
|