Un poco de algebra para transformar un algoritmo - 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(n³), 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+log_p 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.

En esta formulación cada proceso tiene una parte de la lista de datos de entrada, y lo primero que se hace es obtener una cuadrupla aplicando dos veces el operador 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.

Una vez que cada proceso ha realizado su reducción, un proceso central puede coger todas estas cuadruplas y realizar una nueva reducción sobre ellas. Tras esto, el primer elemento de la cuadrupla contiene el resultado de nuestro problema, y para obtenerlo no tenemos más que aplicar el operador pi, que nos devuelve el elemento indicado de una lista.

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;%dn"e;, &n);
    
       printf("e;Número de elementos: %dn"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;%lfn"e;, &data[i]);
       }
    
       printf("e;Datos de la entradan====================n"e;);
       for(i=0;i<n;i++) printf("e;%lfn"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: %dn"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]: %dtdispls[%d]: %dn"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]: %gn"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;======================================nn"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: %lfnn"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.