Un poco de algebra para transformar un algoritmo
Con esto del doctorado que estoy realizando me ha tocado en una asignatura recordar cosas de álgebra que ya tenía bastante olvidades, pero que me ha permitido valorarlas de una forma muy distinta a como lo hice en su día cuando tuve que estudiarlas. En aquel momento la asignatura de álgebra y en menor medida también la de física, las vi como un mal menor que tenía que sufrir para pasar de primer curso y poder avanzar en la carrera estudiando otras cosas que seguramente me iban a gustar más, como así fue de hecho. Desde aquel entonces no había vuelto a entrar en contacto, al menos de forma consciente, con nada de álgebra.
Esta semana en el curso de Computación Distribuida de Altas Prestaciones hemos estado aprendiendo cómo funcionan las librerías MPI y realizando algún ejemplo de utilización. De forma muy elemental, estas librerías permiten realizar programas basados en paso de mensajes entre procesos, que permiten resolver problemas utilizando varios procesos, que pueden residir en diferentes máquinas, trabajando de forma conjunta. Una vez compilado el programa es suficiente utilizar un programa especial indicándole las máquinas y el número de procesos que queremos generar, para que este programa cree los procesos necesarios en cada una de las máquinas indicadas. La parte complicada por supuesto, es reformular el problema de forma que se adapte lo mejor posible a la resolución en paralelo por varios procesos, lo cual será más o menos sencillo en función del problema.
El problema que más me ha llamado la atención y que ha provocado este post, es el del maximum segment sum
(mss). En este problema se dispone de una entrada que consta de una lista de números reales y se trata de obtener la mayor suma de cualquiera de sus segmentos. Por ejemplo, para los datos de entrada 1,-2,3,-4,6,2
tendríamos los segmentos (1,-2), (-2,3), (3,-4), (-4,6), (6,2), (1,-2,3), (-2,3,-4), (3,-4,6), (-4,6,2), (1,-2,3,-4), (-2,3,-4,6), (3,-4,6,2), (1,-2,3,-4,6), (-2,3,-4,6,2), (1,-2,3,-4,6,2)
y el mss sería 8, correspondiente al segmento (6,2)
.
Utilizando notación de programación funcional la forma de resolver el problema que hemos empleado, sería la siguiente:
Siendo segs:
En esta formulación segs
es la pieza que nos permite obtener todo el conjunto posible de segmentos. Sobre esta lista de elementos la primera de las fórmulas hace una reducción utilizando el operador suma, de forma que obtendríamos para cada posible segmento su suma, con lo que únicamente nos queda obtener la mayor de ellas.
Un algoritmo descrito de esta forma tendría una coste computacional de O(n3), siendo n el tamaño del problema de entrada. Pues bien, lo interesante de todo esto, es que utilizando un poco de algebra es posible transformar este algoritmo en otro que nos permite resolver el problema con un coste de O(n+logp p), utilizando el concepto de Homomorfismo sobre listas.
Una vez dicho esto, a partir de la formulación inicial del problema podemos realizar transformaciones que nos llevan finalmente a la siguiente solución, preparada ya para trabajar en varios procesos.
par
(hacer pareja) a cada elemento, por ejemplo, a partir de la entrada (2,5), obtendríamos ((2,2),(2,2)),((5,5),(5,5)). Tras esto cada proceso aplica el operador de reducción a sus cuadruplas, de forma que como resultado cada procesador obtiene una nueva cuadrupla. El operador de reducción actúa de la siguiente forma sobre dos cuadruplas, siendo ((r,s),(t,u)) una cuadrupla.
La implementación de este algoritmo en C y utilizando las librerías MPI es sencilla.
#include <stdio.h> #include <mpi.h> #define WORLD MPI_COMM_WORLD #define BLOCK_LOW(id,p,n)((id)*(n)/(p)) #define BLOCK_SIZE(id,p,n)(BLOCK_LOW((id)+1,p,n)-BLOCK_LOW(id,p,n)) //#define MAX(A,B) ( (A) > (B) ? (A):(B)) double MAX(double a, double b){ return (a > b) ? a : b; } typedef struct{ double r; double s; double t; double u; } Cuadrupla; Cuadrupla hmss_f(); Cuadrupla hmss_operation(); void hmss_reduce(); // proceso 0 carga número de datos -> n // proceso 0 broadcast n // proceso 0 carga datos // El resto de los procesos espera con Barrier // Cada proceso calcula el número de elementos a procesar // Cada proceso envía el número de elementos a 0 (Gather) // ScaterV proceso 0 para enviar bloques al resto de procesos // Cada proceso trabaja con sus datos // Se hace una reducción con un operador propio sobre Cuadruplas int main(int argc, char* argv[]){ MPI_Status status; FILE* fh; int id, p; int n; double* data; int i; int local_n; // para calcular el número de elementos a tratar por el proceso int* size_p; // para almacenar el tamaño del problema para cada proceso int* displs; // para los desplazamientos en MPI_Scatterv double* dist_p; // para almacenar los datos del problema local del proceso Cuadrupla mss_local; Cuadrupla mss_global; MPI_Datatype mpi_type; MPI_Op mpi_op; MPI_Init(&argc, &argv); MPI_Comm_size(WORLD, &p); MPI_Comm_rank(WORLD, &id); MPI_Type_contiguous(4, MPI_DOUBLE, &mpi_type); MPI_Type_commit(&mpi_type); MPI_Op_create(hmss_reduce, 0, &mpi_op); // NO ES CONMUTATIVA, ASÍ QUE UN 0 EN EL 2º PARÁMETRO if(id==0){ // El proceso 0 es el encargado de cargar los datos // y enviar el número de elementos al resto de procesos // abro el archivo if((fh=fopen(argv[1],"e;r"e;))==NULL){ printf("e;No se puede abrir archivo %s"e;, argv[1]); return(0); } // leo el número de elementos fscanf(fh, "e;%d\n"e;, &n); printf("e;Número de elementos: %d\n"e;, n); // creo la matriz de datos data = (double*)malloc(sizeof(double)* n); // leo los datos del archivo, no hago ningún tipo de comprobación for(i=0;i<n;i++){ fscanf(fh, "e;%lf\n"e;, &data[i]); } printf("e;Datos de la entrada\n====================\n"e;); for(i=0;i<n;i++) printf("e;%lf\n"e;,data[i]); } // se hace un broadcast de n desde root MPI_Bcast(&n, 1, MPI_INT, 0, WORLD); // este es el tamaño del problema que tiene que resolver este proceso local_n = BLOCK_SIZE(id,p,n); // creo la matriz para almacenar los datos a evaluar, que me los enviará 0 dist_p = (double*)malloc(sizeof(double)*local_n); printf("e;Proceso[%d] valor de local_n: %d\n"e;, id, local_n);fflush(stdout); if(id == 0){ size_p = (int*)malloc(sizeof(int)*(p)); displs = (int*)malloc(sizeof(int)*(p)); } MPI_Gather(&local_n, 1, MPI_INT, size_p, 1, MPI_INT, 0, WORLD); if(id==0){ displs[0] = 0; for(i=1;i<p;i++){ displs[i] = displs[i-1]+size_p[i-1]; printf("e;\tsize_p[%d]: %d\tdispls[%d]: %d\n"e;, i, size_p[i], i, displs[i]);fflush(stdout); } } // distribuyo a cada proceso la porción de problema que le corresponde MPI_Scatterv(data, size_p, displs, MPI_DOUBLE, dist_p, local_n, MPI_DOUBLE, 0, WORLD); for(i=0;i<local_n;i++){ printf("e;Proceso[%d] - dist_p[%d]: %g\n"e;, id, i, dist_p[i]);fflush(stdout); } // { cada proceso ya tiene su trabajo listo en dist_p } // ahora hay que aplicar la reducción a la lista de Cuadruplas, // obteniendo como resultado otra Cuadrupla mss_local = hmss_f(dist_p, local_n, id); printf("e;Proceso[%d] mss_local: ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf))\n"e;, id, mss_local.r, mss_local.s, mss_local.t, mss_local.u);fflush(stdout); // ahora hago la reducción MPI_Reduce(&mss_local, &mss_global, 1, mpi_type, mpi_op, 0, WORLD); MPI_Finalize(); if(id == 0){ // ahora hago (pi1·pi1) printf("e;======================================\n\n"e;); printf("e;mss_global: ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf))\n"e;, mss_global.r, mss_global.s, mss_global.t, mss_global.u); printf("e;Resultado: %lf\n\n"e;, mss_global.r); } } // Esta función crea las Cuadruplas a partir de los elementos // de la lista y aplica el operador hmss_operation sobre sus // elementos obteniendo como resultado una Cuadrupla Cuadrupla hmss_f(double* list, int n, int id){ int i; Cuadrupla* rstu; // Cuadruplas locales generadas a partir de dist_p Cuadrupla res = {0.0,0.0,0.0,0.0}; Cuadrupla r1, r2, r3, r4, r5, r6; // Creo las Cuadruplas: map(par·par) => ((r,s),(t,u)) // y voy operando sobre ellas rstu = (Cuadrupla*)malloc(sizeof(Cuadrupla)*n); for(i=0;i<n;i++){ rstu[i].r = list[i]; rstu[i].s = list[i]; rstu[i].t = list[i]; rstu[i].u = list[i]; printf("e;Proceso[%d] cuadrupla %d: ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf))\n"e;, id, i, rstu[i].r, rstu[i].s, rstu[i].t, rstu[i].u); if(i==0) res = rstu[i]; else res = hmss_operation(res, rstu[i]); } printf("e;Proceso[%d] hmss_f res: ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf))\n"e;, id, res.r, res.s, res.t, res.u); return res; } // Esta función efectúa la operación hmss sobre dos Cuadruplas y // devuelve la Cuadrupla resultante Cuadrupla hmss_operation(Cuadrupla a, Cuadrupla b){ int i; Cuadrupla res = {0.0,0.0,0.0,0.0}; res.r = MAX(a.r, MAX(b.r, (a.t + b.s))); res.s = MAX(a.s, (a.u + b.s)); res.t = MAX(b.t, (a.t + b.u)); res.u = a.u + b.u; /* printf("e;hmss_operation( ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf)), ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf)) ) => ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf))\n"e;, a.r, a.s, a.t, a.u, b.r, b.s, b.t, b.u, res.r, res.s, res.t, res.u); */ return res; } void hmss_reduce(Cuadrupla* invec, Cuadrupla* inoutvec, int* n, MPI_Datatype *dptr){ Cuadrupla c = hmss_operation(invec[0], inoutvec[0]); printf("e;hmss_reduce( ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf)), ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf)) ) => ((%2.2lf, %2.2lf), (%2.2lf, %2.2lf))\n"e;, invec[0].r, invec[0].s, invec[0].t, invec[0].u, inoutvec[0].r, inoutvec[0].s, inoutvec[0].t, inoutvec[0].u, c.r, c.s, c.t, c.u); fflush(stdout); *inoutvec = c; } |
Hace mucho que no programaba en C, así que el código no es ninguna maravilla, y además he aprovechado para probar alguna cosa más del paso de mensajes, por lo que se podría evitar por ejemplo el primer paso de comunicación en el que cada proceso calcula el tamaño del problema que le va a tocar calcular.
En cualquier caso, a continuación muestro unas trazas de ejecución para los datos mostrados inicialmente,1,-2,3,-4,6,2
, primero con un proceso, y luego con dos y con tres:
1 proceso
Número de elementos: 6 Datos de la entrada ==================== 1.000000 -2.000000 3.000000 -4.000000 6.000000 2.000000 Proceso[0] valor de local_n: 6 size_p[0]: 6 displs[0]: 0 Proceso[0] - dist_p[0]: 1 Proceso[0] - dist_p[1]: -2 Proceso[0] - dist_p[2]: 3 Proceso[0] - dist_p[3]: -4 Proceso[0] - dist_p[4]: 6 Proceso[0] - dist_p[5]: 2 Proceso[0] cuadrupla 0: ((1.00, 1.00), (1.00, 1.00)) Proceso[0] cuadrupla 1: ((-2.00, -2.00), (-2.00, -2.00)) Proceso[0] cuadrupla 2: ((3.00, 3.00), (3.00, 3.00)) Proceso[0] cuadrupla 3: ((-4.00, -4.00), (-4.00, -4.00)) Proceso[0] cuadrupla 4: ((6.00, 6.00), (6.00, 6.00)) Proceso[0] cuadrupla 5: ((2.00, 2.00), (2.00, 2.00)) Proceso[0] hmss_f res: ((8.00, 6.00), (8.00, 6.00)) Proceso[0] mss_local: ((8.00, 6.00), (8.00, 6.00)) ====================================== mss_global: ((8.00, 6.00), (8.00, 6.00)) Resultado: 8.000000 |
2 proceso
Número de elementos: 6 Datos de la entrada ==================== 1.000000 -2.000000 3.000000 -4.000000 6.000000 2.000000 Proceso[0] valor de local_n: 3 Proceso[1] valor de local_n: 3 size_p[0]: 3 displs[0]: 0 size_p[1]: 3 displs[1]: 3 Proceso[0] - dist_p[0]: 1 Proceso[0] - dist_p[1]: -2 Proceso[0] - dist_p[2]: 3 Proceso[0] cuadrupla 0: ((1.00, 1.00), (1.00, 1.00)) Proceso[0] cuadrupla 1: ((-2.00, -2.00), (-2.00, -2.00)) Proceso[0] cuadrupla 2: ((3.00, 3.00), (3.00, 3.00)) Proceso[0] hmss_f res: ((3.00, 2.00), (3.00, 2.00)) Proceso[0] mss_local: ((3.00, 2.00), (3.00, 2.00)) hmss_reduce( ((3.00, 2.00), (3.00, 2.00)), ((8.00, 4.00), (8.00, 4.00)) ) => ((8.00, 6.00), (8.00, 6.00)) Proceso[1] - dist_p[0]: -4 Proceso[1] - dist_p[1]: 6 Proceso[1] - dist_p[2]: 2 Proceso[1] cuadrupla 0: ((-4.00, -4.00), (-4.00, -4.00)) Proceso[1] cuadrupla 1: ((6.00, 6.00), (6.00, 6.00)) Proceso[1] cuadrupla 2: ((2.00, 2.00), (2.00, 2.00)) Proceso[1] hmss_f res: ((8.00, 4.00), (8.00, 4.00)) Proceso[1] mss_local: ((8.00, 4.00), (8.00, 4.00)) ====================================== mss_global: ((8.00, 6.00), (8.00, 6.00)) Resultado: 8.000000 |
3 proceso
Número de elementos: 6 Datos de la entrada ==================== 1.000000 -2.000000 3.000000 -4.000000 6.000000 2.000000 Proceso[0] valor de local_n: 2 Proceso[1] valor de local_n: 2 Proceso[2] valor de local_n: 2 size_p[0]: 2 displs[0]: 0 size_p[1]: 2 displs[1]: 2 size_p[2]: 2 displs[2]: 4 Proceso[0] - dist_p[0]: 1 Proceso[0] - dist_p[1]: -2 Proceso[0] cuadrupla 0: ((1.00, 1.00), (1.00, 1.00)) Proceso[0] cuadrupla 1: ((-2.00, -2.00), (-2.00, -2.00)) Proceso[0] hmss_f res: ((1.00, 1.00), (-1.00, -1.00)) Proceso[0] mss_local: ((1.00, 1.00), (-1.00, -1.00)) hmss_reduce( ((1.00, 1.00), (-1.00, -1.00)), ((3.00, 3.00), (-1.00, -1.00)) ) => ((3.00, 2.00), (-1.00, -2.00)) hmss_reduce( ((3.00, 2.00), (-1.00, -2.00)), ((8.00, 8.00), (8.00, 8.00)) ) => ((8.00, 6.00), (8.00, 6.00)) Proceso[1] - dist_p[0]: 3 Proceso[1] - dist_p[1]: -4 Proceso[1] cuadrupla 0: ((3.00, 3.00), (3.00, 3.00)) Proceso[1] cuadrupla 1: ((-4.00, -4.00), (-4.00, -4.00)) Proceso[1] hmss_f res: ((3.00, 3.00), (-1.00, -1.00)) Proceso[1] mss_local: ((3.00, 3.00), (-1.00, -1.00)) Proceso[2] - dist_p[0]: 6 Proceso[2] - dist_p[1]: 2 Proceso[2] cuadrupla 0: ((6.00, 6.00), (6.00, 6.00)) Proceso[2] cuadrupla 1: ((2.00, 2.00), (2.00, 2.00)) Proceso[2] hmss_f res: ((8.00, 8.00), (8.00, 8.00)) Proceso[2] mss_local: ((8.00, 8.00), (8.00, 8.00)) ====================================== mss_global: ((8.00, 6.00), (8.00, 6.00)) Resultado: 8.000000 |
Y hasta aquí hemos llegado, hemos conseguido, partiendo de una solución intuitiva del problema, obtener mediante transformaciones algebraicas otra solución que a primera vista no es nada intuitiva, y que nos permite resolver el problema de una manera mucho más eficiente.
Los comentarios están cerrados.