-
-
Save vicenteguerra/e963f4bd63b4882051a79ed24f07fba6 to your computer and use it in GitHub Desktop.
/* | |
* Polynomial Interpolation * | |
* This program demonstrates a function that performs polynomial | |
* interpolation. The function is taken from "Numerical Recipes | |
* in C", Second Edition, by William H. Press, Saul A. Teukolsky, | |
* William T. Vetterling, and Brian P. Flannery. | |
* | |
*/ | |
#include <math.h> | |
#include <stdlib.h> | |
#include <stdio.h> | |
#include <omp.h> | |
#define N 20 | |
#define X 14.5 | |
#define NUM_THREADS 4 | |
/* Number of function sample points */ | |
/* Interpolate at this value of x */ | |
/* Function 'vector' is used to allocate vectors with subscript | |
range v[nl..nh] */ | |
double *vector (long nl, long nh) | |
{ | |
double *v; | |
v = (double *) malloc(((nh-nl+2)*sizeof(double))); | |
return v-nl+1; | |
} | |
/* Function 'free_vector' is used to free up memory allocated | |
with function 'vector' */ | |
void free_vector(double *v, long nl, long nh) | |
{ | |
free ((char *) (v+nl-1)); | |
} | |
/* Function 'polint' performs a polynomial interpolation */ | |
void polint (double xa[], double ya[], int n, double x, double *y, double *dy) { | |
int i, m, ns=1; | |
double den,dif,dift,ho,hp,w; | |
double *c, *d; | |
dif = fabs(x-xa[1]); | |
c = vector(1,n); | |
d = vector(1,n); | |
for (i=1; i<= n; i++) { | |
dift = fabs (x - xa[i]); | |
if (dift<dif) { | |
ns = i; | |
dif = dift; | |
} | |
c[i] = ya[i]; | |
d[i] = ya[i]; | |
} | |
*y = ya[ns--]; | |
for (m = 1; m < n; m++) { | |
for (i = 1; i<= n-m; i++) { | |
ho = xa[i] - x; | |
hp = xa[i+m] - x; | |
w = c[i+1] - d[i]; den = ho - hp; den = w / den; d[i] = hp * den; c[i] = ho * den; } | |
*y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--])); } | |
free_vector (d, 1, n); | |
free_vector (c, 1, n); | |
} | |
/* Functions 'sign' and 'init' are used to initialize the | |
x and y vectors holding known values of the function. | |
*/ | |
int sign (int j) | |
{ | |
if (j % 2 == 0) return 1; | |
else return -1; | |
} | |
void init (int i, double *x, double *y) | |
{ | |
// int j; | |
*x = (double) i; | |
*y = sin(i); | |
} | |
/* Function 'main' demonstrates the polynomial interpolation function | |
by generating some test points and then calling 'polint' with a | |
value of x between two of the test points. */ | |
int main (int argc, char *argv[]) | |
{ | |
double x, y, dy; | |
double *xa, *ya; | |
int i; | |
double start_time = omp_get_wtime(); | |
omp_set_num_threads(NUM_THREADS); //define el numero de hilos | |
xa = vector (1, N); | |
ya = vector (1, N); | |
/* Initialize xa's and ya's */ | |
#pragma omp parallel | |
{ | |
#pragma omp for | |
for (i = 1; i<= N; i++) { | |
init (i, &xa[i], &ya[i]); | |
printf ("f(%4.2f) = %13.11f\n", xa[i], ya[i]); | |
} | |
} | |
/* Interpolate polynomial at X */ | |
polint (xa, ya, N, X, &y, &dy); | |
printf ("\nf(%6.3f) = %13.11f with error bound %13.11f\n", X, y, fabs(dy)); free_vector (xa, 1, N); | |
free_vector (ya, 1, N); | |
printf("Time: \t %f \n", omp_get_wtime()-start_time); | |
return 0; | |
} |
/*
- Polynomial Interpolation *
- This program demonstrates a function that performs polynomial
- interpolation. The function is taken from "Numerical Recipes
- in C", Second Edition, by William H. Press, Saul A. Teukolsky,
- William T. Vetterling, and Brian P. Flannery.
*/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <omp.h>
#define N 20
#define X 14.5
#define NUM_THREADS 7
/* Number of function sample points /
/ Interpolate at this value of x /
/ Function 'vector' is used to allocate vectors with subscript
range v[nl..nh] */
double *vector (long nl, long nh)
{
double *v;
v = (double *) malloc(((nh-nl+2)sizeof(double)));
return v-nl+1;
}
/ Function 'free_vector' is used to free up memory allocated
with function 'vector' */
void free_vector(double *v, long nl, long nh)
{
free ((char ) (v+nl-1));
}
/ Function 'polint' performs a polynomial interpolation */
void polint (double xa[], double ya[], int n, double x, double *y, double *dy) {
int i, m, ns=1;
double den,dif,dift,ho,hp,w;
double *c, *d;
dif = fabs(x-xa[1]);
c = vector(1,n);
d = vector(1,n);
#pragma omp parallel
#pragma omp parallel for private (i)
for (i=1; i<= n; i++) {
#pragma omp critical
dift = fabs (x - xa[i]);
if (dift<dif) {
ns = i;
dif = dift;
}
c[i] = ya[i];
d[i] = ya[i];
}
*y = ya[ns--];
#pragma omp parallel
#pragma omp parallel for private (i)
for (m = 1; m < n; m++) {
for (i = 1; i<= n-m; i++) {
#pragma omp critical
ho = xa[i] - x;
hp = xa[i+m] - x;
w = c[i+1] - d[i]; den = ho - hp; den = w / den; d[i] = hp * den; c[i] = ho * den; }
*y += (dy=(2ns < (n-m) ? c[ns+1] : d[ns--])); }
free_vector (d, 1, n);
free_vector (c, 1, n);
}
/* Functions 'sign' and 'init' are used to initialize the
x and y vectors holding known values of the function.
*/
int sign (int j)
{
if (j % 2 == 0) return 1;
else return -1;
}
void init (int i, double *x, double *y)
{
// int j;
*x = (double) i;
y = sin(i);
}
/ Function 'main' demonstrates the polynomial interpolation function
by generating some test points and then calling 'polint' with a
value of x between two of the test points. */
int main (int argc, char *argv[])
{
double x, y, dy;
double *xa, *ya;
int i;
double start_time = omp_get_wtime();
omp_set_num_threads(NUM_THREADS); //define el numero de hilos
xa = vector (1, N);
ya = vector (1, N);
/* Initialize xa's and ya's */
#pragma omp parallel
{
#pragma omp for
for (i = 1; i<= N; i++) {
#pragma omp critical
init (i, &xa[i], &ya[i]);
printf ("f(%4.2f) = %13.11f\n", xa[i], ya[i]);
}
}
/* Interpolate polynomial at X */
polint (xa, ya, N, X, &y, &dy);
printf ("\nf(%6.3f) = %13.11f with error bound %13.11f\n", X, y, fabs(dy)); free_vector (xa, 1, N);
free_vector (ya, 1, N);
printf("Time: \t %f \n", omp_get_wtime()-start_time);
return 0;
}
0.000243
0.000783
0.000771
0.019338
0.001289
0.000950
0.001554