ps_homeworks/dn4/quad.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;
}