• Aprender formas de repartir información entre procesos
• Aprender a resolver problemas resultados de la repartición de datos
• Comprobar las ventajas de una distribución de trabajo paralelo

Vamos a implementar un programa que resuelva caminos cortos mediante el algoritmo de Floyd.

Problema de los caminos más cortos

Inicialmente se presenta un algoritmo secuencial para resolver el problema de encontrar los caminos más cortos en un grafo.
Sea un grafo etiquetado G=(V, E, long), donde:

• V={vi} es un conjunto de N vértices (|V| = N).
• E está totalmente incluido en el conjunto VxV, que es un conjunto de arcos que conectan vértices de V.
• long: E -> Z es una función que asigna una etiqueta entera a cada arista de E.

Podemos representar un grafo mediante una matriz de adyacencia A en la que:

Un camino desde el vértice vi al vértice vj es una secuencia de vértices (vi, vk), (vk, vl), ..., (vm, vj), donde ningún vértice aparece más de una vez. El camino más corto entre dos vértices vi y vj es el camino cuya suma de las etiquetas en sus aristas es menor.
El problema del camino más corto sencillo requiere encontrar el camino más corto desde un único vértice a todos los demás vértices del grafo. El problema de todos los caminos más cortos requiere encontrar los caminos más cortos entre todos los pares de vértices del grafo. El algoritmo para resolver este último problema toma como entrada la matriz de incidencia A y calcula la matriz S de tamaño N x N, donde S(i,j) es la longitud del camino más corto desde vi, a vj, o un valor infinito si no hay camino entre dichos vértices.

Floyd secuencial

El algoritmo de Floyd deriva la matriz S en N pasos, construyendo en cada paso K una matriz intermedia I(k) con el camino más corto conocido entre cada par de nodos. Inicialmente:
Ii,j(0) = Ai,j
El k-ésimo paso del algoritmo considera cada Ii,j y determina si el camino más corto conocido desde vi a vj es mayor que las longitudes combinadas de los caminos desde vi a vk y desde vk a vj, en cuyo caso se actualizará la entrada Ii,j. La operación de comparación se realiza un total de N^3 veces (N al cubo), por lo que aproximamos el coste secuencial del algoritmo como tcN^3, siendo tc el coste de una operación de comparación.

procedure floyd sequential
begin
  Ii,j = 0 si i=j
  Ii,j = long(vi, vk) si (vi, vj) pertenece a E
  Ii,j = infinito en otro caso
  for k=0 to N-1
    for i=0 to N-1
      for j=0 to N-1
        Ii,j(k+1) = min(Ii,j(k), Ii,k(k)+ Ik,j(k))
      endfor
    endfor
  endfor
  S = I(N)
end

Floyd descomposición Unidimensional

Se asume que el número de vértices N es múltiplo del número de procesos P.
La primera versión paralela del algoritmo de Floyd está basada en una descomposición unidimensional por bloques de filas de la matriz intermedia I y de la matriz de salida S. Podremos utilizar N procesadores como máximo. Cada tarea es responsable de una o más filas adyacentes de I y ejecutará el siguiente código:

for k=0 to N-1
  for i=local_i_start to local_i_end
    for j=0 to N-1
      Ii,j(k+1) = min(Ii,j(k), Ii,k(k)+ Ik,j(k))
    endfor
  endfor
endfor

En la iteración k, cada tarea, además de sus datos locales, necesita los valores Ik,0, Ik,1, ..., Ik,N-1. Es decir, la fila k de I. Por ello la tarea que tenga asignada la fila K deberá difundirla (MPI_Broadcast) a todas las demás.

Descomposición ID de la matriz I entre 4 procesadores para la primera versión paralela del algoritmo de Floyd.

Para distribuir la matriz A entre los procesadores basta con utilizar la operación colectiva MPI_Scatter.

Floyd descomposición Bidimensional

Se asume que el número de vértices N es múltiplo de la raíz del número de procesos P. Es decir, si hay 9 procesos, el número de vértices (N) debe ser múltiplo de 3 (Raíz de 9 = 3).
Esta versión del algoritmo de Floyd utiliza una descomposición bidimensional, pudiendo utilizar hasta N^2 (N elevado a 2) procesadores:

for k=0 to N-1
  for i=local_i_start to local_i_end
    for j= local_j_start to local_j_end
      Ii,j(k+1) = min(Ii,j(k), Ii,k(k)+ Ik,j(k))
    endfor
  endfor
endfor

En cada paso, además de los datos locales, cada tarea necesita N/Raiz(P) valores de dos tareas localizadas en la misma fila y en la misma columna respectivamente (Ver imagen de descomposición bidimensional). Por tanto, los requerimientos de comunicación en la etapa K son dos operaciones de broadcast (MPI_Bcast):

• Desde la tarea en cada fila que contiene parte de la columna k a todas las demás tareas en esta fila y.
• Desde la tarea en cada columna que contiene parte de la fila k a todas las demás tareas en esta columna.

En cada uno de los N pasos, N/Raiz(P) valores deben ser difundidos a las Raiz(P) tareas en cada fila y en cada columna. Nótese que cada tarea debe servir como el origen para al menos un broadcast a cada tarea en la misma fila y a cada columna del array 2D.

En el caso general, cada uno de los P procesadores almacena una submatriz de A de tamaño ( N/Raiz(P) x N/Raiz(P) ) (Por simplicidad, supondremos que Raiz(P) divide a N).
Inicialmente el procesador P0 contiene la matriz completa y a cada procesador le correspondería N/Raiz(P) elementos de N/Raiz(P) filas de dicha matriz.
Para permitir que la matriz A pueda ser repartida con una operación colectiva entre los procesadores, podemos definir un tipo de datos para especificar submatrices. Para ello, es necesario considerar que los elementos de una matriz bidimensional se almacenan en posiciones consecutivas de memoria por orden de fila, en primer lugar, y por orden de columna en segundo lugar. Así cada submatriz será un conjunto de N/Raiz(P) bloques de N/Raiz(P) elementos cada uno, con un desplazamiento de N elementos entre cada bloque.
La función MPI_Type_vector permite asociar una submatriz cuadrada como un tipo de datos. De esta manera, el procesador P0 podrá enviar bloques no contiguos de datos a los demás procesadores en un solo mensaje. También es necesario calcular, para cada procesador, la posición de comienzo de la submatriz que le corresponde.
Para poder alojar de forma contigua y ordenada los elementos del nuevo tipo creado (las submatrices cuadradas), con objeto de poder repartirlos con MPI_Scatter entre los procesadores, podemos utilizar la función MPI_Pack. Utilizando esta función, se empaquetan las submatrices de forma consecutiva, de tal forma que al repartirlas (usando el tipo MPI_PACKED), se le asigna una submatriz cuadrada a cada procesador. A continuación, se muestra una porción de código C con la secuencia de operaciones necesarias para empaquetar todos los bloques de una matriz NxN de forma ordenada y repartirlos con un MPI_Scatter entre los procesadores:

MPI_Datatype MPI_BLOQUE;
......
......
 
raiz_P=sqrt(P);
tam=N/raiz_P;
 
/*Creo buffer de envío para almacenar los datos empaquetados*/
buf_envio=reservar_vector(N*N); 
 
if (rank==0)
  {
  /* Obtiene matriz local a repartir*/
  Inicializa_matriz(N,N,matriz_I); 
   /*Defino el tipo bloque cuadrado */
   MPI_Type_vector (tam, tam, N, MPI_INT, &MPI_BLOQUE);
   /* Creo el nuevo tipo */
   MPI_Type_commit (&MPI_BLOQUE);
 
   /* Empaqueta bloque a bloque en el buffer de envío*/
   for (i=0, posicion=0; i<size; i++)
     {
        /* Calculo la posicion de comienzo de cada submatriz */	
        fila_P=i/raiz_P;  
        columna_P=i%raiz_P;
        comienzo=(columna_P*tam)+(fila_P*tam*tam*raiz_P);
        MPI_Pack (matriz_I(comienzo), 1, MPI_BLOQUE, 
                          buf_envio,sizeof(int)*N*N, &posicion, MPI_COMM_WORLD);
     }
 
     /*Destruye la matriz local*/	
     free(matriz_I);
     /* Libero el tipo bloque*/	
     MPI_Type_free (&MPI_BLOQUE);
   }
 
/*Creo un buffer de recepcion*/ 
buf_recep=reservar_vector(tam*tam);
/* Distribuimos la matriz entre los procesos */
MPI_Scatter (buf_envio, sizeof(int)*tam*tam, MPI_PACKED, 
                      buf_recep, tam*tam, MPI_INT, 0, MPI_COMM_WORLD);

Para obtener la matriz resultado en el procesador P0, se utiliza un MPI_Gather seguido de un MPI_Unpack.

El programa recibe por parámetros el fichero donde se encuentran almacenados los caminos a recorrer.
La forma de ejecutar el fichero secuencial es:

Para generar más ficheros de prueba, se dispone de un código denominado "creaejemplo.cpp", este se debe compilar

Una vez creado el ejecutable, se pueden crear ficheros de entrada pasando por parámetros en número de vértices que se desean

Hay que tener que el programa es solo una ayuda, al ser secuencial puede tardar bastante tiempo en generar ejemplos, además de que consume memoria exponencialmente (cuantos más nodos se desean, muchísima más memoria consume).

Los ficheros guardan primero el número de vértices (o nodos) y a continuación se guardan en ternas el vértice origen, el destino y el coste de ese recorrido (en ese orden).

Si analizamos el grafo resultante, podríamos comprobar que la solución es llegar a 2 desde 0, llegar a 3 desde 2, llegar a 1 desde 3 y llegar a 0 desde 1.

floyd.cc

#include <iostream>
#include <fstream>
#include <string.h>
#include <cstdlib>
#include "Graph.h"
#include "mpi.h"
 
using namespace std;
 
//**************************************************************************
 
int main (int argc, char *argv[]){
	MPI::Init(argc,argv);
 
	if (argc != 2) 
		{
		 cerr << "Sintaxis: " << argv[0] << " <archivo de grafo>" << endl;
		 exit(-1);
		}
 
	Graph G;
	G.lee(argv[1]);		// Read the Graph
	cout << "EL Grafo de entrada es:"<<endl;
	G.imprime();
 
	int nverts=G.vertices;
 
	Graph *I=&G;
 
	// BUCLE PPAL DEL ALGORITMO
	int i,j,k,vikj;
	for(k=0;k<nverts;k++){
		for(i=0;i<nverts;i++)
			for(j=0;j<nverts;j++)
				if (i!=j && i!=k && j!=k){
					vikj=G.arista(i,k)+G.arista(k,j);
					vikj=min(vikj,G.arista(i,j));
					I->inserta_arista(i,j,vikj);   
				}
		G=*I;
	}
 
  MPI::Finalize();
 
 
  cout << endl<<"EL Grafo con las distancias de los caminos más cortos es:"<<endl<<endl;
  I->imprime();
}

Graph.h

//**************************************************************************
#ifndef GRAPH_H
#define GRAPH_H
 
#include <cstdlib>
#include <string.h>
//**************************************************************************
const int INF= 100000;
 
//**************************************************************************
class Graph //Adjacency List clas
{
	private:
	   int *A;
	public:
           Graph();
           int vertices;
	   void fija_nverts(const int verts);
	   void inserta_arista(const int ptA,const int ptB, const int edge);
	   int arista(const int ptA,const int ptB);
           void imprime();
           void lee(char *filename);
};
 
//**************************************************************************
#endif
//**************************************************************************

Graph.cc

//***********************************************************************
#include "Graph.h"
#include <iostream>
#include <fstream>
using namespace std;
 
//***********************************************************************
Graph::Graph ()		// Constructor
{
}
//***********************************************************************
void Graph::fija_nverts (const int nverts){
	A=new int[nverts*nverts];
	vertices=nverts;
}
//***********************************************************************
void Graph::inserta_arista(const int vertA, const int vertB, const int edge) // inserta
A->B
{
	A[vertA*vertices+vertB]=edge;
}
//***********************************************************************
int Graph::arista(const int ptA,const int ptB)
{
	return A[ptA*vertices+ptB];
}
//***********************************************************************
void Graph::imprime()
{
	int i,j,vij;
	for(i=0;i<vertices;i++){
		cout << "A["<<i << ",*]= ";
 
		for(j=0;j<vertices;j++){
			if (A[i*vertices+j]==INF) 
				cout << "INF";
			else  
				cout << A[i*vertices+j];
			if (j<vertices-1) 
				cout << ",";
			else
				cout << endl;
		}
	}
}
//***********************************************************************
void Graph::lee(char *filename){
	const int BUF_SIZE 100
	std::ifstream infile(filename);
 
	if (!infile){
		cerr << "Nombre de archivo inválido \"" << filename << "\" !!" << endl;
		cerr << "Saliendo........." << endl;
		exit(-1);
	}
	//Obten el numero de vertices
	char buf[BUF_SIZE];
	infile.getline(buf,BUF_SIZE,'\n');
	vertices=atoi(buf);
	A=new int[vertices*vertices];
 
	int i,j;
	for(i=0;i<vertices;i++)
		for(j=0;j<vertices;j++)
			if (i==j) A[i*vertices+j]=0;
			else A[i*vertices+j]=INF;
 
	while (infile.getline(buf,BUF_SIZE) && infile.good() && !infile.eof()){
		char *vertname2 = strpbrk(buf, " \t");
		*vertname2++ = '\0';
		char *buf2 = strpbrk(vertname2, " \t");
		*buf2++ = '\0';
		int weight = atoi(buf2);
		i=atoi(buf);
		j=atoi(vertname2);
		 A[i*vertices+j]=weight;
	}
}

• Implementar los algoritmos paralelos de cálculo de todos los caminos más cortos que han sido descritos previamente. Se aconseja realizar cambios sobre la clase Graph, que encapsula la gestión del grafo etiquetado, para implementar todas las operaciones de comunicación dentro de dicha clase. Por tanto, la clase encapsularía un grafo distribuido entre los procesadores en base a una determinada distribución (por bloques de filas para la descomposición unidimensional, o por bloques cuadrados para la descomposición bidimensional). Se debe crear un directorio diferente para cada versión paralela (Floyd_1, Floyd_2). Cada directorio mantendrá los mismos archivos que se usan en la versión secuencial pero con diferente código. De esa forma el Makefile se puede reutilizar y se percibe claramente la diferencia entre versiones.
• Realizar medidas de tiempo de ejecución sobre los algoritmos implementados. Para medir tiempos de ejecución, podemos utilizar la función MPI_Wtime. Para asegurarnos de que todos los procesadores comienzan su computación al mismo tiempo, podemos utilizar MPI_Barrier. Las medidas deberán excluir las fases de entrada/salida. Deberán realizarse las siguientes medidas:

(A) Medidas para el algoritmo secuencial (P=1).
(B) Medidas para el algoritmo paralelo (P=4). Las medidas deberán excluir las fases de entrada/salida, así como la fase de distribución inicial de la matriz A desde P0 y la fase de reunión del resultado en P0.
Las medidas deberán realizarse para diferentes tamaños de problema, para así poder comprobar el efecto de la granularidad sobre el rendimiento de los algoritmos. Se presentará una tabla con el siguiente formato: