Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_c_kernel_local -p 8888:8888 -d palmoreck/jupyterlab_c_kernel:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_c_kernel_localDocumentación de la imagen de docker palmoreck/jupyterlab_c_kernel:1.1.0 en liga.
Nota generada a partir de liga
A diferencia con los paquetes de parallel de R o de multiprocessing o dask de Python en los que se sugiere el lanzamiento de subprocesos a partir de un proceso principal para el cómputo en paralelo en tareas principalmente CPU bound (cada subproceso son instancias de Python y R), en el lenguaje C lanzamos threads desde un proceso principal conocido con el nombre de threading en los que los threads realizan fork's y join's vistos en 2.1.Un_poco_de_historia_y_generalidades.
Una herramienta que permite el fork y join de threads a partir de un proceso es OpenMP. Otra más es POSIX threads Pthreads. Revisaremos en esta nota algunos conceptos y usos básicos de OpenMP.
Documentación: openMP
Es una API para shared memory parallel programming. Las siglas MP se refieren a multiprocessing un sinónimo de shared memory parallel computing. Al programar con OpenMP consideramos que en nuestro sistema cada thread y core pueden potencialmente tener acceso a toda la memoria disponible.
Algunas características de OpenMP son:
Paralelización de ciclos for secuenciales en los que las iteraciones son independientes una de la otra de forma simple.
Paralelización de tareas y sincronización explícita de threads.
OpenMP es una API que provee directivas para shared_memory programming. Una directiva es una instrucción especial para indicar al preprocesador (vía la compilación) que ejecutaremos una instrucción que no se encuentra en la especificación básica del lenguaje C (ver C preprocessor). En OpenMP utilizamos el pragma: #pragma omp.
Comentario: los pragma's son añadidos al sistema para comportamientos que no son parte de la especificación básica de C. Cabe señalar que las versiones más recientes del compilador gcc sí soportan a los pragmas y todas las preprocessor directive son por default de longitud una línea. Si no cabe en una línea usamos un escaping (precedemos con un backslash "\" la "nueva línea").
In [13]:
//%cflags:-fopenmp
//%cflags:-Wall
#include<stdio.h>
#include<stdlib.h>
#include<omp.h> //header file con prototipos y macro definitions para
//la librería de funciones y macros de openMP.
void Hello(void); //prototipo de función a ejecutar por los threads.
int main(){
//Siempre iniciamos con un #pragma omp ... :
#pragma omp parallel //directiva parallel
//structured block que sólo consiste del llamado a la función Hello:
Hello();
return 0;
}
//función que será ejecutada por los threads
void Hello(void){
int mi_rango = omp_get_thread_num(); //obtenemos el rank dado por
//el run-time system a cada thread
int conteo_threads = omp_get_num_threads(); //obtenemos el número de threads
//que realizaron un fork
//del master thread
printf("Hola del thread: %d de %d\n", mi_rango, conteo_threads);
}
Obs:
-fopenmp para soporte de OpenMP. De modo que la compilación se realiza:gcc -Wall -fopenmp hello_world_omp.c -o hello_world_omp.out
suponiendo que el programa anterior se guarda en un archivo con nombre hello_world_omp.c.
Dependiendo del número de cores de nuestro sistema tendremos diferentes número de printf's.
Lo que continúa a la línea de #pragma omp parallel es un structured block, esto es, un statement o conjunto de statements que tienen un punto de entrada y un punto de salida, no se permiten statements como el siguiente:
#pragma omp parallel
if(...) break;ni tampoco:
#pragma omp parallel
{
if(variable == valor) return 1;
return -1;
}A continuación del nombre parallel podemos usar diferentes tipos de clauses, una clause en OpenMP es un texto que modifica una directiva. Por ejemplo, podemos usar la clause num_threads para especificar el número de threads que ejecutarán el structured block:
In [16]:
//%cflags:-fopenmp
//%cflags:-Wall
#include<stdio.h>
#include<stdlib.h>
#include<omp.h> //header file con prototipos y macro definitions para
//la librería de funciones y macros de openMP.
void Hello(void); //prototipo de función a ejecutar por los threads.
int main(){
int numero_threads=5;
//Siempre iniciamos con un #pragma omp ... :
#pragma omp parallel num_threads(numero_threads) //directive parallel con clause num_threads
//structured block que sólo consiste del llamado a la función Hello:
Hello();
return 0;
}
//función que será ejecutada por los threads
void Hello(void){
int mi_rango = omp_get_thread_num(); //obtenemos el rank dado por
//el run-time system a cada thread
int conteo_threads = omp_get_num_threads(); //obtenemos el número de threads
//que realizaron un fork
//del master thread
printf("Hola del thread: %d de %d\n", mi_rango, conteo_threads);
}
Obs: obsérvese el no determinismo en el printf.
Forma secuencial
In [1]:
//%cflags:-lm
#include<stdio.h>
#include<stdlib.h>
#include<math.h> //header para funciones de mate
#include<time.h>
#include <sys/time.h>
void Rcf(double a, double b, int n,\
double *suma_global_p);
double f(double nodo);
int main(int argc, char *argv[]){
double suma_global = 0.0;
double a=0.0, b=1.0;
int n=1e6; //número de subintervalos
double objetivo=0.7468241328124271;
struct timeval start;
struct timeval end;
long seconds;
long long mili;
gettimeofday(&start, NULL);
Rcf(a,b,n,&suma_global);
gettimeofday(&end, NULL);
seconds = (end.tv_sec - start.tv_sec);
mili = 1000*(seconds) + (end.tv_usec - start.tv_usec)/1000;
printf("Integral de %f a %f = %1.15e\n", a,b,suma_global);
printf("Error relativo de la solución: %1.15e\n", fabs(suma_global-objetivo)/fabs(objetivo));
printf("Tiempo de ejecución: %lld miliseconds", mili);
return 0;
}
void Rcf(double a, double b, int n, double *sum){
double h_hat=(b-a)/n;
double x=0.0;
int i=0;
*sum = 0.0;
for(i=0;i<=n-1;i++){
x = a+(i+1/2.0)*h_hat;
*sum+=f(x);
}
*sum =h_hat*(*sum);
}
double f(double nodo){
double valor_f;
valor_f = exp(-pow(nodo,2));
return valor_f;
}
Forma en paralelo:
Nota: los siguientes resultados se obtuvieron con una máquina con 8 cores, así que pueden no coincidir con los resultados previos de esta sección.
In [2]:
//%cflags:-fopenmp
//%cflags:-lm
//%cflags:-Wall
#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
#include<math.h> //header para funciones de mate
#include<time.h>
#include <sys/time.h>
//prototipo de funciones
double Rcf_parallel(double a, double h_hat, int ns_p);
double f(double nodo);
//programa principal
int main(){
double suma_global = 0.0; //variable que es shared, en esta se sumarán los resultados
//de las aproximaciones de cada uno de los threads.
double a=0.0, b=1.0;
int n=1e6; //número de subintervalos
double h_hat=(b-a)/n;
int ns_p;
int numero_threads[5]={0};//conteo_threads debe dividir de forma exacta a n
int long_numero_threads = 0;
double objetivo=0.7468241328124271;
struct timeval start;
struct timeval end;
long seconds;
long long mili;
int i;
numero_threads[0]=1;
numero_threads[1]=2;
numero_threads[2]=4;
numero_threads[3]=5;
numero_threads[4]=8;
long_numero_threads = sizeof(numero_threads)/sizeof(numero_threads[0]);
for(i=0;i<long_numero_threads;i++){
ns_p=n/numero_threads[i];
gettimeofday(&start, NULL);
#pragma omp parallel num_threads(numero_threads[i]) reduction(+: suma_global)
suma_global+=Rcf_parallel(a,h_hat,ns_p);
suma_global = h_hat*suma_global;
gettimeofday(&end, NULL);
seconds = (end.tv_sec - start.tv_sec);
mili = 1000*(seconds) + (end.tv_usec - start.tv_usec)/1000;
printf("Integral de %f a %f = %1.15e\n", a,b,suma_global);
printf("Error relativo de la solución: %1.15e\n", fabs(suma_global-objetivo)/fabs(objetivo));
printf("Tiempo de ejecución con %d threads: %lld miliseconds\n", numero_threads[i],mili);
printf("----------------------\n");
suma_global=0.0;
}
return 0;
}
//definición de funciones:
double Rcf_parallel(double a, double h_hat, int ns_p){
int begin, end;
int mi_rango = omp_get_thread_num();
double local_int=0;
int i;
double x;
begin = mi_rango*ns_p;
end = begin + ns_p;
for(i=begin;i<=end-1;i++){
x = a+(i+1/2.0)*h_hat;
local_int+=f(x);
}
return local_int;
}
double f(double nodo){
double valor_f;
valor_f = exp(-pow(nodo,2));
return valor_f;
}
Ejercicio: Calcular tiempo de ejecución para la regla de Simpson. No olviden medir errores relativos. Tal regla está en 1.5.Integracion_numerica,