1. Indice de pagina
  2. 1. Objetivos
  3. 2. Descripción
  4. 2.1. Solución mediante Branch and Bound
  5. 2.2. Algoritmo paralelo
  6. 2.3. El programa principal
  7. 2.4. Equilibrado de carga
  8. 2.5. Detección de fin en anillo
  9. 2.6. Difusión de la cota superior
  10. 2.7. Aspectos de implementación con MPI
  11. 3. Modo de uso
  12. 4. Código inicial
  13. 4.1. libbb.h
  14. 4.2. libbb.cc
  15. 4.3. bbseq.cc
  16. 5. Ejercicios propuestos

Objetivos

Descripción

Vamos a implementar un programa que resuelva el problema del viajante de comercio de forma paralela y distribuida.

Solución mediante Branch and Bound

El problema del viajante de comercio (Traveling Salesman Problem, TSP) puede ser representado como un grafo dirigido compuesto de un conjunto de vértices (ciudades) y arcos etiquetados (distancias entre ciudades). Una solución optimizada al TSP es un camino de coste mínimo en el cual todos los vértices son visitados exactamente una vez.
El problema se resuelve mediante un algoritmo de acotación y ramificación (Branch-and-Bound, BB) que, dinámicamente, construye un árbol de búsqueda cuya raíz es el problema inicial y sus nodos hoja son caminos entre todas las ciudades (no necesariamente de coste óptimo).
En los algoritmos BB se definen las siguientes funciones:

Cálculo de la cota inferior

El cálculo de la cota inferior de un nodo del árbol de búsqueda se lleva a cabo obteniendo la matriz de coste reducido. Para ello introducimos las siguientes definiciones:

  1. Una fila o columna de una matriz se dice reducida si tiene al menos una entrada cero y las demás entradas son no negativas.
  2. Una matriz se dice reducida si todas sus filas y columnas están reducidas.

Si se resta un valor v a todas las entradas de una fila o columna de la matriz del coste del (sub)problema representado por la matriz se reduce exactamente en v. Si P es un subproblema entonces su cota inferior ci(P) es la suma de las cantidades restadas a sus filas y columnas para obtener la matriz de coste reducido. Un recorrido completo de costo mínimo continúa siéndolo tras la operación de sustración. De esta forma, toda la información relativa a un subproblema está contenida en su matriz de coste reducido y en su valor de cota inferior.
Ejemplo. Supongamos que partimos de un problema A a partir del cual se obtiene una matriz reducida B:

Imagen explicacion matriz reducida

La matriz B se obtiene sustrayendo 10, 2, 2, 3, 4 de las filas 1, 2, 3, 4 y 5, así como 1 y 3 de las columnas 1 y 3. Al final el total sustraido es de 25. Esto significa que todos los recorridos completos del TSP para el grafo que representan las matrices A y B, tienen al menos 25 unidades (ci(A) = 25).

Construcción del árbol de búsqueda

Para un nodo del árbol de búsqueda A, se generan dos subproblemas independientes mediante la elección de un arco (i, j):

Nos interesa obtener árboles derechos que puedan ser acotados pronto, para así reducir el espacio de búsqueda. Para ello utilizamos una heurística que elija un arco (i, j) que produzca un hijo derecho con un valor de cota inferior tan grande como sea posible. Si elegimos un arco con coste positivo, el valor de cota inferior del hijo derecho no se incrementará ya que su matriz de coste seguirá siendo reducida al eliminar la entrada (i,j). Si elegimos un arco con coste nulo, al eliminar dicha entrada quizás tengamos que restar alguna cantidad a las entradas de la fila i y de la columna j para obtener la matriz de coste reducido. El valor de cota inferior para el hijo derecho se incrementará en

Imagen formula minimo

ya que ésto es lo que se necesita sustraer de la fila i y de la columna j para introducir un 0 en ambas. Esta heurística permite maximizar el valor de cota inferior para el subárbol derecho.

Algoritmo paralelo

Tras un estudio de los aspectos de descomposición y asignación, llegamos a un diseño de programa paralelo con las siguientes características:

  1. Para pedir nuevos nodos (cuando su parte del árbol de búsqueda haya sido ya explorada.
  2. Para ofrecer parte de sus nodos (cuando reciba peticiones de otros procesos).

Nuestra implementación del algoritmo de BB paralelo consistirá en un único proceso por procesador. La utilización de varios comunicadores y de variantes asíncronas de las primitivas de paso de mensajes nos van a permitir implementar los diferentes requerimientos de comunicación de forma modular.

El programa principal

Cada proceso ejecuta el siguiente algoritmo:

main () {
  U = INFINITO; // inicializa la cota superior
  if(id == 0){
    Leer_Problema_Inicial(&nodo);
  }else{
    Equilibrar_Carga(&pila, &fin);
    if(!fin) Pop(&pila, &nodo);
  }
 
  while(!fin){   // ciclo del Branch&Bound
    Ramifica(&nodo, &nodo_izq, &nodo_dch);
    if(Solucion(&nodo_dch)){
      if(ci(nodo_dch) < U) U = ci(nodo_dch); //Actualiza cota sup
    }else{  // Si no es un nodo hoja
      if(ci(nodo_dch) < U) Push(&pila, &nodo_dch);
    }
 
    if(Solucion(&nodo_izq)){
      if(ci(nodo_izq) < U) U = ci(nodo_izq); // Actualiza cota sup
    }else{ // No es un nodo hoja
      if(ci(nodo_izq) < U) Push(&pila, &nodo_izq);
    }
 
    hay_nueva_cota_superior = Difusion_Cota_Superior(&U);
    if(hay_nueva_cota_superior) Acotar(&pila, U);
 
    Equilibrado_Carga(&pila, &fin);
    if(!fin) Pop(&pila, &nodo);
  }
}

Inicialmente sólo el proceso 0 expande el problema inicial. El resto de procesos llaman al procedimiento de equilibrado de carga, lo que provocará que se realicen peticiones de trabajo al proceso 0. Según se vayan generando nuevos nodos, las peticiones de trabajo se irán atendiendo y los procesos obtendrán nodos para expandir.
En el ciclo principal del algoritmo, la expansión de un nodo genera dos nuevos nodos, nodo_izq y nodo_dch. Si un nodo generado es una solución al problema se comprueba si su coste mejora la cota superior actual, U, en cuyo caso se actualiza ésta. Si el nodo generado no es una solución y su cota inferior es menor que la cota superior actual, se almacena para su posterior expansión.
Cada proceso mantiene localmente una pila de nodos para expandir. La elección de una estructura de datos tipo pila permite que el árbol de búsqueda sea explorado por cada proceso de acuerdo con un esquema primero en profundidad, para así encontrar soluciones lo antes posible. Obsérvese que en el proceso, el almacenamiento del nodo derecho se lleva a cabo antes que el del nodo izquierdo, con el fin de expandir primero éste último y llegar antes a una solución.
Cuando un proceso produce un nuevo valor de cota superior, dicho valor deberá ser difundido para que sea conocido por los demás procesos. Por ello, cada proceso deberá comprobar si algún otro proceso generó un nuevo valor de cota superior. El procedimiento Difundir_Cota_Superior implementa estas cuestiones.
Si se generó un nuevo valor de cota superior (local o remotamente) podremos eliminar de la pila aquellos nodos con valor de la cota inferior mayor que la nueva cota superior (procedimiento Acotar).
Antes de extraer un nuevo nodo de la pila (procedimiento Pop), la llamada a Equilibrar_Carga permitirá pedir nodos a otros procesos si la pila está vacía, o bien ceder nodos a otros procesos si se reciben peticiones (y hay nodos para dar) o detectar la situación de fin.

Equilibrado de carga

Cuando la pila de un proceso está vacía, se envía un mensaje de petición de trabajo, indicando el identificador del proceso, id, como el origen de la petición. Las peticiones de trabajo circulan por los procesos siguiendo una estructura de anillo, por lo que el mensaje es enviado al siguiente proceso (id+1)%P. A continuación el proceso entra en un bucle del que saldrá cuando haya conseguido trabajo para seguir ejecutándose, o bien cuando, al fallar sus peticiones, se detecte la situación de fin. El proceso espera recibir un mensaje de respuesta. Si recibe un mensaje de trabajo, la petición tuvo éxito, almacenándose los nodos recibidos en la pila local. Si recibe un mensaje de petición (desde el proceso anterior (id-1+P)%P), esta petición no podrá ser atendida al estar vacía la pila del proceso, por lo que dicho mensaje es reenviado al siguiente proceso (id+1)%P. Además si el origen de la petición es el mismo proceso, la petición habrá viajado por todo el anillo sin poder ser atendida, por lo que podríamos encontrarnos en una situación de fin (ver sección 4.5).
Si no se llegó a la situación de fin, el proceso sondea mensajes de petición de trabajo de otros procesos. Si no hay mensajes pendientes el procedimiento de equilibrado de carga acaba y el proceso puede continuar ejecutando el algoritmo de BB. Mientras haya peticiones de trabajo, el proceso intenta dividir su pila local para ceder nodos al proceso solicitante. Si la pila no tiene nodos suficientes para ceder, reenvía las peticiones al proceso siguiente (podemos decidir no ceder nodos si el tamaño de la pila está por debajo de un umbral, por ejemplo: la función que divide la pila falla si no hay al menos dos nodos en la pila).

Equilibrado_Carga(tPila *pila, bool *fin){
  if(Vacia(pila)){   // el proceso no tiene trabajo, pide a otros procesos
    Enviar peticion de trabajo al proceso (id+1)%P;
    while(Vacia(pila) && !fin){
      Esperar mensaje de otro proceso;
      switch(tipo de mensaje){
        case PETICION:    // Peticion de trabajo
          Recibir mensaje de peticion de trabajo;
          Reenviar peticion de trabajo al proceso (id+1)%P;
          if(solicitante == id){  // Peticion devuelta
            Iniciar deteccion de posible situacion de fin;
          }
          break;
        case NODOS:   // Resultado de una peticion de trabajo
          Recibir nodos del proceso donante;
          Almacenar nodos recibidos en la pila;
          break;
      }
    }
  }
  if(!fin){ // el proceso tiene nodos para trabajar
    Sondear si hay mensajes pendientes de otros procesos;
    while(hay mensajes){ // atiende peticiones mientras haya mensajes
      Recibir mensaje de peticion de trabajo;
      if(hay suficientes nodos en la pila para ceder)
        Enviar nodos al proceso solicitante;
      else
        Reenviar peticion de trabajo al proceso (id+1)%P;
 
      Sondear si hay mensajes pendientes de otros procesos;
    }
  }
}

Una forma gráfica de ver el equilibrado de carga es la siguiente.

Imagen equilibrado de carga (Pulsar sobre la imagen para ver la animación)

Detección de fin en anillo

El algoritmo de BB acabará cuando el espacio de búsqueda haya sido explorado, es decir, cuando no queden más nodos para ramificar. El programa debe acabar cuando todos los procesos agoten sus pilas locales y no puedan obtener nodos de otros procesos para seguir trabajando.

Algoritmo de terminación de Dijkstra

Supongamos un conjunto de procesos organizados en anillo. Un proceso puede estar en dos estados: activo y pasivo. Cuando el proceso P0 pasa al estado pasivo envía un testigo al anterior proceso del anillo: Pn-1. Cuando un proceso recibe el testigo, si su estado es pasivo, lo enviará al proceso anterior, o, si su estado es activo, lo mantendrá hasta que se pase al estado pasivo. Cuando el proceso P0 reciba el testigo de nuevo, sabrá que todos los demás procesos han terminado.
Este esquema no funciona si un proceso en estado pasivo puede ser reactivado por un mensaje de trabajo enviado por otro proceso. Supongamos que el proceso Pi envía un mensaje de trabajo al proceso Pj (en estado pasivo), para después pasar al estado pasivo. Cuando Pi reciba el testigo, lo pasará al proceso Pi-1. Si i>j, el mensaje de trabajo llegará al proceso Pj antes que el testigo. Dicho mensaje de trabajo reactivará al proceso Pj, por lo que al recibir el testigo lo retendrá. Sin embargo si i<j y el testigo visitó al proceso Pj antes de recibir el mensaje de trabajo, Pj podrá ser reactivado y Pi recibirá el testigo, creyendo que todos los procesos Pk con k>i están pasivos (entre ellos Pj). En este caso se podría producir una situación de fin errónea.

Imagen mala deteccion de fin (Pulsar sobre la imagen para ver la animación)

Para evitar este escenario, se propone la siguiente modificación: cada proceso, además de un estado, tiene un color: blanco o negro (inicialmente blanco). El testigo, también tiene un color (inicialmente blanco).
Al comienzo el proceso P0 posee el testigo de detección de fin. Cuando el proceso P0 pasa al estado pasivo, envía un testigo blanco al proceso Pn-1. Cuando un proceso Pi envía un mensaje de trabajo a otro proceso Pj, siendo i<j, cambia su color a negro. Un proceso pasivo con color blanco reenviará el testigo sin cambiar el color de éste. Un proceso pasivo con color negro reenviará el testigo con color negro; tras ello, el proceso cambiará su color a blanco. El algoritmo terminará cuando el proceso P0, siendo su color blanco, reciba un testigo blanco. Si el proceso P0 recibe un testigo negro, volverá a enviarlo con color blanco. Si el mensaje de trabajo que provocó que el testigo fuera cambiando a negro reactivó a algún proceso, el testigo será retenido por el proceso activo hasta que pase al estado pasivo.

P(id)
  ...
  Esperar evento;
  switch(tipo_evento){
    case MENSAJE_TRABAJO:
      estado = ACTIVO;
    case MENSAJE_PETICION:
      if(hay trabajo para ceder){
        j=origen(PETICION);
        Enviar TRABAJO a P(j);
        if(id<j) mi_color = NEGRO;
      }
    case MENSAJE_TOKEN:
      token_presente = TRUE;
      if(estado == PASIVO){
        if(id==0 && mi_color == BLANCO && color(TOKEN) == BLANCO)
          TERMINACION DETECTADA;
        else{
          if(id==0) color(TOKEN) = BLANCO;
          else if (mi_color == NEGRO) color(TOKEN) = NEGRO;
          Enviar TOKEN a P(id-1);
          mi_color = BLANCO;
          token_presente = FALSE;
        }
      }
    case TRABAJO_AGOTADO:
      estado = PASIVO;
      if(token_presente){
        if(id==0) color(TOKEN) = BLANCO;
        else if(mi_color == NEGRO) color(TOKEN) = NEGRO;
        Enviar TOKEN a P(id-1);
        mi_color = BLANCO;
        token_presente = FALSE;
      }
  }

La representación de la solución dada está recogida en la siguiente animación.

Imagen deteccion de fin (Pulsar sobre la imagen para ver la animación)

Integración del algoritmo de detección de fin en el procedimiento de equilibrado de carga

Un proceso se considerará pasivo al agotar su pila local y emitir una petición de trabajo, ésta sea devuelta al no poder ser atendida por ningún otro proceso.
Cuando un proceso pase al estado pasivo volverá a enviar el mensaje de petición de trabajo al siguiente proceso, y además, si tiene el testigo, reiniciará la detección de fin enviándolo al proceso anterior en el anillo.
Cuando un proceso activo reciba el testigo de detección de fin, lo guardará para reenviarlo cuando pase al estado pasivo.
Cuando el proceso 0 detecte la situación de fin, recogerá los mensajes de petición de trabajo que haya circulando, y emitirá un segundo testigo de confirmación para que los demás procesos conozcan dicha situación y puedan terminar correctamente. Los demás procesos, según vayan recibiendo el testigo de confirmación, lo retransmitirán al siguiente proceso y terminarán. El proceso 0 terminará cuando reciba este testigo.

Difusión de la cota superior

La difusión y mantenimiento de valores globales de cota superior se lleva a cabo haciendo circular dichos valores por los procesos en una estructura de comunicación en anillo. Es conveniente que la difusión de la cota superior se realice de forma asíncrona con respecto a la computación local de cada proceso. No es necesario que todos los procesos tengan en cada momento el mismo valor de cota superior: en este algoritmo la consistencia no afecta a la corrección. Solamente se requiere que los valores estén tan actualizados como sea posible para así optimizar la búsqueda.
Cuando un proceso tiene un nuevo valor de cota superior, deberá difundir dicho valor enviándolo al siguiente proceso del anillo. A continuación, el proceso comprueba la existencia de mensajes pendientes desde el proceso anterior del anillo. Si hay algún mensaje de cota superior, lo recibe, actualiza su valor local al menor de ambos, y reenvía el mensaje al proceso siguiente del anillo.
Para evitar un número excesivo de mensajes y llegar así a una situación de interbloqueo, necesitamos alguna forma de asegurar que cada valor de cota superior dé una sola vuelta al anillo, y que en cada momento, haya como máximo un solo mensaje por proceso circulando. Para ello, cada valor de cota superior enviado al anillo irá acompañado del identificador del proceso que lo generó. Cuando un proceso reciba un valor de cota superior enviado por él mismo, no lo reenviará hacia el siguiente proceso del anillo. Además, cada proceso deberá registrar si está pendiente de recibir desde el proceso anterior del anillo aǵun mensaje enviado por él mismo:
Un proceso no difundirá un nuevo mensaje de cota superior (encontrado por él mismo) mientras no reciba el mensaje anteriormente enviado por él mismo. Cuando un proceso recibe un valor desde el proceso anterior del anillo, actualiza su valor local y lo reenvía si su origen era otro proceso. Si el origen del mensaje era el mismo proceso que lo ha recibido, se considera terminada la difusión del valor de cota superior, habilitándose la difusión de nuevos valores.

Difusion_Cota_Superior(){
  if(difundir_cs_local && !pendiente_retorno_cs){
    Enviar valor local de cs al proceso (id+1)%P;
    pendiente_retorno_cs = TRUE;
    difundir_cs_local = FALSE;
  }
  Sondear si hay mensajes de cota superior pendientes;
  while(hay mensajes){
    Recibir mensaje con valor de cota superior desde el proceso (id-1+P)%P;
    Actualizar valor local de cota superior;
 
    if(origen_mensaje == id && difundir_cs_local){
      Enviar valor local de cs al proceso (id+1)%P;
      pendiente_retorno_cs = TRUE;
      difundir_cs_local = FALSE;
    }else if(origen_mensaje == id && !difundir_cs_local)
      pendiente_retorno_cs = FALSE;
    else              // origen mensaje == otro proceso
      Reenviar mensaje al proceso (id+1)%P;
 
    Sondear si hay mensajes de cota superior pendientes;
  }
}

Aspectos de implementación con MPI

Para implementar nuestro programa de forma modular podemos seguir algunas recomendaciones.
Además del programa principal, implementaremos dos funciones separadas, una para el equilibrado de carga/detección de fin, y otra para la difusión de la cota superior. Para evitar interferencias entre los mensajes de estas funciones, utilizaremos comunicadores separados para ellas. Esto nos permite, por ejemplo, esperar mensajes relacionados con el equilibrado de carga (utilizando MPI_ANY_SOURCE y/o MPI_ANY_TAG) independientemente de los mensajes de difusión de la cota superior.
Dentro de cada función utilizaremos etiquetas diferentes para cada tipo de mensaje. Así en la función de equilibrado de carga podríamos tener los siguientes tipos de mensaje: petición, trabajo, testigo de detección de fin y fin detectado. En la función de difusión de cota superior podríamos utilizar la etiqueta para identificar al origen del mensaje de la cota superior.
Es importante mantener adecuadamente la información de estado de cada proceso. El estado de un proceso podría incluir:

Prevención de interbloqueos

El algoritmo de BB paralelo muestra un esquema de comunicación asíncrono que debe ser implementado utilizando paso de mensajes. La situación de interbloqueo se produce cuando se forma un ciclo de procesos en el que cada proceso está bloqueado intentando enviar un mensaje al siguiente.
Para evitar los interbloqueos necesitamos limitar el número máximo de mensajes que los procesos pueden enviar. En el caso del equilibrado de carga, como máximo hay un mensaje por proceso: una petición de trabajo que, en caso de ser atendida, es cambiada por un mensaje de trabajo. Simultáneamente con los mensajes de trabajo, podría haber circulando un testigo de detección de fin. En el caso de la difusión de la cota superior, como máximo hay tantos mensajes como procesos, ya que nunca un proceso hace circular simultáneamente más de un valor de cota superior.
Por otro lado, es necesario suministrar el espacio de búfer suficiente para albergar el número máximo de mensajes permitidos. En nuestra aplicación, aún limitando el número de mensajes que cada proceso puede enviar, necesitamos cierta cantidad de buferización para evitar interbloqueos.

Modo de uso

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

./bbseq X tsp_problems/tspX.1

Código inicial

libbb.h

/* ************************************************************************ */
/*    Libreria para Branch-Bound, manejo de pila y Gestion de memoria       */
/* ************************************************************************ */
 
#include <cstdio>
const unsigned int MAXPILA =150;
 
#ifndef NULO_T
#define NULO_T
const int NULO = -1;
#endif
 
#ifndef INFINITO_T
#define INFINITO_T
const long int INFINITO = 100000;
#endif
 
extern unsigned int NCIUDADES;
 
/* ************************************************************************ */
/* ******** Representacion del problema T.S.P. **************************** */
/* ************************************************************************ */
 
/* ************************************************************************ */
/*
 # Una instancia del problema TSP es una matriz bidimensional de
   Nciudades X Nciudades entradas, donde:
     * las entradas (i,j) corresponden a la distancia desde la ciudad
       i a la ciudad j (siendo 0<=i<Nciudades, 0<=j<Nciudades, i<>j)
     * las entradas (i,i) tienen valor infinito (0<=i<N.ciudades)
 
 # VECTOR DE DISTANCIAS (REPRESENTACION COMPLETA DEL TSP)
   ------------------------------------------------------
   Los procedimientos de la libreria libbb representan la matriz del TSP
   como un vector de distancias:
         int tsp[NCIUDADES][NCIUDADES];
   donde la entrada (i,j) corresponde al elemento tsp[i][j]
 
 
 # DESCRIPCION ABREVIADA DE TRABAJO.
   ---------------------------------
   Representa un (sub)problema del TSP, y contiene informacion de los
   arcos includos e incluidos.
   Una descripcion completa del (sub)problema (vector de distancias) se
   reconstruye a partir de la descripcion abreviada y del vector de
   distancias del problema inicial.
 
   Una descripcion abreviada de trabajo se almacena en la siguiente
   estructura:
                struct sNodo {
                  int   ci;
                  int   incl[NCIUDADES];
                  int   orig_excl;
                  int   dest_excl[NCIUDADES-2];
                };
                typedef struct sNodo tNodo;
                tNodo nodo;
   donde:
    * nodo.ci  es la cota inferior del subproblema
    * nodo.incl[0], nodo.incl[1], ... , nodo.incl[NCIUDADES-1]
      tienen el siguiente significado:
      - si nodo.incl[i] = j   (0<=j<NCIUDADES)
        se ha incluido el arco <i,j>
      - si nodo.incl[i] = NULO
        no hay ningun arco <i,*> incluido en el camino
      - cualquier otro valor sera considerado erroneo
   *  Los arcos excluidos explicitamente en la busqueda corresponden siempre
      a la misma fila (la busqueda se realiza de ciudad en ciudad).
      Cuando en una fila i tengamos (explicita o implicitamente) excluidos
      NCIUDADES-2 arcos (sin contar el arco <i,i>) esto supondra la inclusion
      forzosa del arco que queda.
      Por tanto solo tenemos que anotar como maximo NCIUDADES-2 arcos
      explicitamente excluidos.
      - si nodo.orig_excl = k  (0<=k<N.ciudades)
        los arcos excluidos explicitamente son de la forma <k,*>
      - si  nodo.orig_excl = NULO
        no quedan mas arcos por excluir explicitamente
   * nodo.dest_excl[0], nodo.dest_excl[1], ... ,nodo.dest_excl[NCIUDADES-2]
     tienen el siguiente significado:
     - si nodo.dest_excl[l]=NULO
       no corresponde a ningun arco excluido
     - si nodo.dest_excl[l]=m  (0<=m<NCIUDADES)
       se ha excluido (de forma explicita) el arco <nodo.orig_excl , m>
     - cualquier otro valor se considera erroneo
*/
/* *********************************************************************** */
 
#ifndef TARCO_T
#define TARCO_T
struct tArco {
	int   v;
	int   w;
};
#endif
 
 
#ifndef TNODO_T
#define TNODO_T
class tNodo {
	public:
		int ci;
		int* incl;
		int orig_excl;
		int* dest_excl;
 
		tNodo() {incl = new int[NCIUDADES]; dest_excl = new int[NCIUDADES-2];};
		~tNodo() {delete [] incl; delete [] dest_excl;};
};
#endif
 
 
 
/* ******************************************************************* */
/*  Pila de nodos para el problema del viajante (TSP) con estrategia   */
/*  primero en profundidad.                                            */
/* ******************************************************************* */
/*                                                                     */
/*  Una pila es almacenada, declarada e inicializada de la forma       */
/*  siguiente:                                                         */
/*                                                                     */
/*        struct sPila {                                               */
/*          int tope;                                                  */
/*          struct sNodo nodos[MAXPILA];                               */
/*        };                                                           */
/*        typedef struct sPila tPila;                                  */
/*        sPila pila_nodos;                                            */
/*        PilaInic (*pila_nodos);                                      */
/*                                                                     */
/*  Cada elemento de pila_nodos almacena un nodo del arbol             */
/*  de busqueda (una descripcion abreviadada de trabajo).              */
/*                                                                     */
/***********************************************************************/
 
#ifndef TPILA_T
#define TPILA_T
struct tPila {
	int tope;
	tNodo nodos[MAXPILA];
};
#endif
 
 
/* ********************************************************************* */
/* *** Cabeceras de funciones para el algoritmo de Branch-and-Bound *** */
/* ********************************************************************* */
 
void LeerMatriz (char archivo[], int** tsp) ;
 
bool Inconsistente  (int** tsp);
  /* tsp   -  matriz de inicidencia del problema o subproblema            */
  /* Si un subproblema tiene en alguna fila o columna todas sus entradas  */
  /* a infinito entonces es inconsistente: no conduce a ninguna solucion  */
  /* del TSP                                                              */
 
void Reduce (int** tsp, int *ci);
  /* tsp   -  matriz de incidencia del problema o subproblema        */
  /* ci    -  cota Inferior del problema                             */
  /* Reduce la matriz tsp e incrementa ci con las cantidades         */
  /* restadas a las filas/columnas de la matriz tsp.                 */
  /* La matriz reducida del TSP se obtiene restando a cada           */
  /* entrada de cada fila/columna de la matriz de incidencia,        */
  /* la entrada minima de dicha fila/columna. La cota inferior del   */
  /* (sub)problema es la suma de las cantidades restadas en todas    */
  /* las filas/columnas.                                             */
 
bool EligeArco (tNodo *nodo, int** tsp, tArco *arco);
  /*  nodo  -  descripcion de trabajo de un nodo del arbol (subproblema)   */
  /*  tsp   -  matriz de incidencia del subproblema (se supone reducida)   */
  /*  Busca una arco con valor 0 en una fila de la que no se haya incluido */
  /*  todavia ningun arco.                                                 */
 
void IncluyeArco(tNodo *nodo, tArco arco);
  /* nodo    - descripcion abreviada de trabajo                  */
  /* arco    - arco a incluir en la descripcion de trabajo       */
  /* Incluye el arco 'arco' en la descripcion de trabajo 'nodo'  */
 
bool ExcluyeArco(tNodo *nodo, tArco arco);
  /* nodo  - descripcion abreviada de trabajo                           */
  /* arco  - arco a excluir (explicitam) en la descripcion de trabajo   */
  /* Excluye el arco en la descripcion de trabajo 'nodo'                */
 
void PonArco(int** tsp, tArco arco);
  /*  tsp   -   matriz de incidencia                             */
  /*  Pone las entradas <v,?> y <?,w> a infinito, excepto <v,w>  */
 
void QuitaArco(int** tsp, tArco arco);
  /* tsp   -    matriz de incidencia                          */
  /* Pone la entrada <v,w> a infinito (excluye este arco)     */
 
void EliminaCiclos(tNodo *nodo, int** tsp);
  /* Elimina en 'tsp' los posibles ciclos que puedan formar los arcos */
  /* incluidos en 'nodo'                                              */
 
void ApuntaArcos(tNodo *nodo, int** tsp);
/* Dada una descripcion de trabajo 'nodo' y una matriz de inicidencia 'tsp' */
/* llama a Pon.Arco() para los arcos incluidos en 'nodo' y a Quita.Arco()   */
/* para los excluidos. Despues llama a Elimina.Ciclos para eliminar ciclos  */
/* en los caminos                                                           */
 
void InfiereArcos(tNodo *nodo, int** tsp);
  /* Infiere nuevos arcos a incluir en 'nodo' para aquellas filas que  */
  /* tengan N.ciudades-2 arcos a infinito en 'tsp'                     */
 
void Reconstruye (tNodo *nodo, int** tsp0, int** tsp);
  /* A partir de la descripcion del problema inicial 'tsp0' y de la  */
  /* descripcion abreviada de trabajo 'nodo', construye la matriz de */
  /* incidencia reducida 'tsp' y la cota Inferior 'ci'.              */
 
void HijoIzq (tNodo *nodo, tNodo *lnodo, int** tsp, tArco arco);
  /* Dada la descripcion de trabajo 'nodo', la matriz de incidencia        */
  /* reducida 'tsp', construye la descripcion de trabajo 'l.nodo'          */
  /* a partir de la inclusion del arco                                     */
 
void HijoDch (tNodo *nodo, tNodo *rnodo, int** tsp, tArco arco);
  /* Dada la descripcion de trabajo 'nodo' y la matriz de incidencia */
  /* reducida 'tsp', construye la descripcion de trabajo 'r.nodo'    */
  /* a partir de la exclusion del arco <v,w>.                        */
 
void Ramifica (tNodo *nodo, tNodo *rnodo, tNodo *lnodo, int** tsp0);
  /* Expande nodo, obteniedo el hijo izquierdo lnodo    */
  /* y el hijo derecho rnodo                            */
 
bool Solucion(tNodo *nodo);
  /* Devuelve TRUE si la descripcion de trabajo 'nodo' corresponde a una */
  /* solucion del problema (si tiene N.ciudades arcos incluidos).        */
 
int Tamanio (tNodo *nodo);
  /* Devuelve el tamanio del subproblema de la descripcion de trabajo nodo */
  /* (numero de arcos que quedan por incluir)                              */
 
void InicNodo (tNodo *nodo);
  /* Inicializa la descripcion de trabajo 'nodo' */
 
void CopiaNodo (tNodo *origen, tNodo *destino);
  /* Copia una descripcion de trabajo origen en otra destino */
 
void EscribeNodo (tNodo *nodo);
  /* Escribe en pantalla el contenido de la descripcion de trabajo nodo */
 
 
/* ********************************************************************* */
/* ***   Cabeceras de Funciones para manejo de la pila de nodos      *** */
/* ********************************************************************* */
 
void PilaInic (tPila *pila);
  /* Inicializa la pila */
bool PilaLlena (tPila *pila);
  /* Devuelve TRUE si la pila esta llena, FALSE en caso contrario */
bool PilaVacia (tPila *pila);
  /* Devuelve TRUE si la pila esta vacia, FALSE en caso contrario */
int PilaTamanio (tPila *pila);
  /* Devuelve el numero de elementos en la pila */
 
bool PilaPush (tPila *pila, tNodo *nodo);
  /* Inserta nuevo elemento en la pila.       */
  /* devuelve FALSE si pila llena.            */
 
bool PilaPop (tPila *pila, tNodo *nodo);
  /* Saca un elemento de la pila.  */
  /* devuelve FALSE si pila vacia.  */
 
bool PilaDivide (tPila *pila1, tPila *pila2);
  /* Divide pila1 en dos (pila1 y pila2) tomando elementos desde el      */
  /* fondo al tope.                                                      */
  /* Si la pila esta vacia o contiene un solo elemento devuelve FALSE.   */
 
void PilaAcotar (tPila *pila, int U);
  /* Elimina los elementos que tengan valor de cota inferior >= U */
 
 
/* ******************************************************************** */
// Funciones para reserva dinamica de memoria
int ** reservarMatrizCuadrada(unsigned int orden);
void liberarMatriz(int** m);
/* ******************************************************************** */

libbb.cc

/* ************************************************************************ */
/*  Libreria de funciones para el Branch-Bound y manejo de la pila          */
/* ************************************************************************ */
 
#include <cstdio>
#include <cstdlib>
#include "libbb.h"
 
extern unsigned int NCIUDADES;
 
/* ********************************************************************* */
/* ****************** Funciones para el Branch-Bound  ********************* */
/* ********************************************************************* */
 
void LeerMatriz (char archivo[], int** tsp) {
  FILE *fp;
  int i, j;
 
  if (!(fp = fopen(archivo, "r" ))) {
    printf ("ERROR abriendo archivo %s en modo lectura.\n", archivo);
    exit(1);
  }
  printf ("-------------------------------------------------------------\n");
  for (i=0; i<NCIUDADES; i++) {
    for (j=0; j<NCIUDADES; j++) {
      fscanf( fp, "%d", &tsp[i][j] );
      printf ("%3d", tsp[i][j]);
    }
    fscanf (fp, "\n");
    printf ("\n");
  }
  printf ("-------------------------------------------------------------\n");
}
 
 
bool Inconsistente (int** tsp) {
  int  fila, columna;
  for (fila=0; fila<NCIUDADES; fila++) {   /* examina cada fila */
    int i, n_infinitos;
    for (i=0, n_infinitos=0; i<NCIUDADES; i++)
      if (tsp[fila][i]==INFINITO && i!=fila)
        n_infinitos ++;
    if (n_infinitos == NCIUDADES-1)
      return true;
  }
  for (columna=0; columna<NCIUDADES; columna++) { /* examina columnas */
    int i, n_infinitos;
    for (i=0, n_infinitos=0; i<NCIUDADES; i++)
      if (tsp[columna][i]==INFINITO && i!=columna)
        n_infinitos++;               /* increm el num de infinitos */
    if (n_infinitos == NCIUDADES-1)
      return true;
  }
  return false;
}
 
void Reduce (int** tsp, int *ci) {
  int min, v, w;
  for (v=0; v<NCIUDADES; v++) {
    for (w=0, min=INFINITO; w<NCIUDADES; w++)
      if (tsp[v][w] < min && v!=w)
        min = tsp[v][w];
    if (min!=0) {
      for (w=0; w<NCIUDADES; w++)
        if (tsp[v][w] != INFINITO && v!=w)
          tsp[v][w] -= min;
      *ci += min;       /* acumula el total restado para calc c.i. */
    }
  }
  for (w=0; w<NCIUDADES; w++) {
    for (v=0, min=INFINITO; v<NCIUDADES; v++)
      if (tsp[v][w] < min && v!=w)
        min = tsp[v][w];
    if (min !=0) {
      for (v=0; v<NCIUDADES; v++)
        if (tsp[v][w] != INFINITO && v!=w)
          tsp[v][w] -= min;
      *ci += min;     /* acumula cantidad restada en ci */
    }
  }
}
 
bool EligeArco (tNodo *nodo, int** tsp, tArco *arco) {
  int i, j;
  for (i=0; i<NCIUDADES; i++)
    if (nodo->incl[i] == NULO)
      for (j=0; j<NCIUDADES; j++)
        if (tsp[i][j] == 0 && i!=j) {
          arco->v = i;
          arco->w = j;
          return true;
        }
  return false;
}
 
void IncluyeArco(tNodo *nodo, tArco arco) {
  nodo->incl[arco.v] = arco.w;
  if (nodo->orig_excl == arco.v) {
    int i;
    nodo->orig_excl++;
    for (i=0; i<NCIUDADES-2; i++)
      nodo->dest_excl[i] = NULO;
  }
}
 
 
bool ExcluyeArco(tNodo *nodo, tArco arco) {
  int i;
  if (nodo->orig_excl != arco.v)
    return false;
  for (i=0; i<NCIUDADES-2; i++)
    if (nodo->dest_excl[i]==NULO) {
      nodo->dest_excl[i] = arco.w;
      return true;
    }
  return false;
}
 
void PonArco(int** tsp, tArco arco) {
  int j;
  for (j=0; j<NCIUDADES; j++) {
    if (j!=arco.w)
      tsp[arco.v][j] = INFINITO;
    if (j!=arco.v)
      tsp[j][arco.w] = INFINITO;
  }
}
 
void QuitaArco(int** tsp, tArco arco) {
  tsp[arco.v][arco.w] = INFINITO;
}
 
void EliminaCiclos(tNodo *nodo, int** tsp) {
  int cnt, i, j;
  for (i=0; i<NCIUDADES; i++)
    for (cnt=2, j=nodo->incl[i]; j!=NULO && cnt<NCIUDADES;
         cnt++, j=nodo->incl[j])
      tsp[j][i] = INFINITO; /* pone <nodo[j],i> infinito */
}
 
void ApuntaArcos(tNodo *nodo, int** tsp) {
  int i;
  tArco arco;
 
  for (arco.v=0; arco.v<NCIUDADES; arco.v++)
    if ((arco.w=nodo->incl[arco.v]) != NULO)
      PonArco (tsp, arco);
  for (arco.v=nodo->orig_excl, i=0; i<NCIUDADES-2; i++)
    if ((arco.w=nodo->dest_excl[i]) != NULO)
      QuitaArco (tsp, arco);
  EliminaCiclos (nodo, tsp);
}
 
void InfiereArcos(tNodo *nodo, int** tsp) {
  bool cambio;
  int cont, i, j;
  tArco arco;
 
  do {
    cambio = false;
    for  (i=0; i<NCIUDADES; i++)     /* para cada fila i */
      if (nodo->incl[i] == NULO) {   /* si no hay incluido un arco <i,?> */
        for (cont=0, j=0; cont<=1 && j<NCIUDADES; j++)
          if (tsp[i][j] != INFINITO && i!=j) {
            cont++;  /* contabiliza entradas <i,?> no-INFINITO */
            arco.v = i;
            arco.w = j;
          }
        if (cont==1) {  /* hay una sola entrada <i,?> no-INFINITO */
          IncluyeArco(nodo, arco);
          PonArco(tsp, arco);
          EliminaCiclos (nodo, tsp);
          cambio = true;
        }
      }
  } while (cambio);
}
 
void Reconstruye (tNodo *nodo, int** tsp0, int** tsp) {
  int i, j;
  for (i=0; i<NCIUDADES; i++)
    for (j=0; j<NCIUDADES; j++)
      tsp[i][j] = tsp0[i][j];
  ApuntaArcos (nodo, tsp);
  EliminaCiclos (nodo, tsp);
  nodo->ci = 0;
  Reduce (tsp,&nodo->ci);
}
 
void HijoIzq (tNodo *nodo, tNodo *lnodo, int** tsp, tArco arco) {
  int** tsp2 = reservarMatrizCuadrada(NCIUDADES);;
  int i, j;
 
  CopiaNodo (nodo, lnodo);
  for (i=0; i<NCIUDADES; i++)
    for (j=0; j<NCIUDADES; j++)
      tsp2[i][j] = tsp[i][j];
  IncluyeArco(lnodo, arco);
  ApuntaArcos(lnodo, tsp2);
  InfiereArcos(lnodo, tsp2);
  Reduce(tsp2, &lnodo->ci);
  liberarMatriz(tsp2);
}
 
void HijoDch (tNodo *nodo, tNodo *rnodo, int** tsp, tArco arco) {
  int** tsp2 = reservarMatrizCuadrada(NCIUDADES);
  int i, j;
 
  CopiaNodo (nodo, rnodo);
  for (i=0; i<NCIUDADES; i++)
    for (j=0; j<NCIUDADES; j++)
      tsp2[i][j] = tsp[i][j];
  ExcluyeArco(rnodo, arco);
  ApuntaArcos(rnodo, tsp2);
  InfiereArcos(rnodo, tsp2);
  Reduce(tsp2, &rnodo->ci);
 
	liberarMatriz(tsp2);
 
}
 
void Ramifica (tNodo *nodo, tNodo *lnodo, tNodo *rnodo, int** tsp0) {
  int** tsp = reservarMatrizCuadrada(NCIUDADES);
  tArco arco;
  Reconstruye (nodo, tsp0, tsp);
  EligeArco (nodo, tsp, &arco);
  HijoIzq (nodo, lnodo, tsp, arco);
  HijoDch (nodo, rnodo, tsp, arco);
 
	liberarMatriz(tsp);
 
}
 
bool Solucion(tNodo *nodo) {
  int i;
  for (i=0; i<NCIUDADES; i++)
    if (nodo->incl[i] == NULO)
      return false;
  return true;
}
 
int Tamanio (tNodo *nodo) {
  int i, cont;
  for (i=0, cont=0; i<NCIUDADES; i++)
    if (nodo->incl[i] == NULO)
      cont++;
  return cont;
}
 
void InicNodo (tNodo *nodo) {
  int i;
  nodo->ci = 0;
  for (i=0; i<NCIUDADES; i++)
    nodo->incl[i] = NULO;
  nodo->orig_excl = 0;
  for (i=0; i<NCIUDADES-2; i++)
    nodo->dest_excl[i] = NULO;
}
 
void CopiaNodo (tNodo *origen, tNodo *destino) {
  int i;
  destino->ci = origen->ci;
  for (i=0; i<NCIUDADES; i++)
    destino->incl[i] = origen->incl[i];
  destino->orig_excl = origen->orig_excl;
  for (i=0; i<NCIUDADES-2; i++)
    destino->dest_excl[i] = origen->dest_excl[i];
}
 
void EscribeNodo (tNodo *nodo) {
  int i;
  printf ("ci=%d : ",nodo->ci);
  for (i=0; i<NCIUDADES; i++)
    if (nodo->incl[i] != NULO)
      printf ("<%d,%d> ",i,nodo->incl[i]);
  if (nodo->orig_excl < NCIUDADES)
    for (i=0; i<NCIUDADES-2; i++)
      if (nodo->dest_excl[i] != NULO)
        printf ("!<%d,%d> ",nodo->orig_excl,nodo->dest_excl[i]);
  printf ("\n");
}
 
 
 
/* ********************************************************************* */
/* ********** Funciones para manejo de la pila  de nodos *************** */
/* ********************************************************************* */
void PilaInic (tPila *pila) {
  pila->tope = 0;
}
 
bool PilaLlena (tPila *pila) {
  return (pila->tope == MAXPILA);
}
 
bool PilaVacia (tPila *pila) {
  return (pila->tope == 0);
}
 
int PilaTamanio (tPila *pila) {
  return pila->tope;
}
 
bool PilaPush (tPila *pila, tNodo *nodo) {
  if (PilaLlena (pila))
    return false;
  CopiaNodo (nodo, &pila->nodos[pila->tope]);
  pila->tope++;
  return true;
}
 
bool PilaPop (tPila *pila, tNodo *nodo) {
  if (PilaVacia(pila))
    return false;
  pila->tope--;
  CopiaNodo (&pila->nodos[pila->tope], nodo);
  return true;
}
 
 
bool PilaDivide (tPila *pila1, tPila *pila2) {
  int i;
 
  if (PilaVacia(pila1) || PilaTamanio(pila1)==1)
    return false;
  for (i=0; i<pila1->tope/2; i++) {
    CopiaNodo (&pila1->nodos[i*2], &pila1->nodos[i]);
    CopiaNodo (&pila1->nodos[i*2+1], &pila2->nodos[i]);
  }
  if (pila1->tope%2 == 0) {
    pila1->tope = pila1->tope/2;
    pila2->tope = pila1->tope;
  }
  else {
    CopiaNodo (&pila1->nodos[pila1->tope-1], &pila1->nodos[pila1->tope/2]);
    pila1->tope = pila1->tope/2 + 1;
    pila2->tope = pila1->tope - 1;
  }
  return true;
}
 
void PilaAcotar (tPila *pila, int U) {
  int tope2;
  int i;
 
  for (i=0, tope2=0; i<pila->tope; i++)
    if (pila->nodos[i].ci < U) {
      CopiaNodo (&pila->nodos[i], &pila->nodos[tope2]); // correccion 19/2/99
      tope2++;
    }
  pila->tope = tope2;
}
 
 
/* ******************************************************************** */
//         Funciones de reserva dinamica de memoria
/* ******************************************************************** */
 
// Reserva en el HEAP una matriz cuadrada de dimension "orden".
int ** reservarMatrizCuadrada(unsigned int orden) {
	int** m = new int*[orden];
	m[0] = new int[orden*orden];
	for (unsigned int i = 1; i < orden; i++) {
		m[i] = m[i-1] + orden;
	}
 
	return m;
}
 
// Libera la memoria dinamica usada por matriz "m"
void liberarMatriz(int** m) {
	delete [] m[0];
	delete [] m;
}

bbseq.cc

/* ******************************************************************** */
/*               Algoritmo Branch-And-Bound Secuencial                  */
/* ******************************************************************** */
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <mpi.h>
#include "libbb.h"
 
using namespace std;
 
unsigned int NCIUDADES;
 
main (int argc, char **argv) {
        MPI::Init(argc,argv);
	switch (argc) {
		case 3:		NCIUDADES = atoi(argv[1]);
					break;
		default:	cerr << "La sintaxis es: bbseq <tamaño> <archivo>" << endl;
					exit(1);
					break;
	}
 
	int** tsp0 = reservarMatrizCuadrada(NCIUDADES);
	tNodo	nodo,         // nodo a explorar
			lnodo,        // hijo izquierdo
			rnodo,        // hijo derecho
			solucion;     // mejor solucion
	bool activo,        // condicion de fin
		nueva_U;       // hay nuevo valor de c.s.
	int  U;             // valor de c.s.
	tPila pila;         // pila de nodos a explorar
 
	U = INFINITO;                  // inicializa cota superior
	InicNodo (&nodo);              // inicializa estructura nodo
	PilaInic (&pila);              // inicializa pila
	LeerMatriz (argv[2], tsp0);    // lee matriz de fichero
	activo = !Inconsistente(tsp0);
 
	while (activo) {       // ciclo del Branch&Bound
		Ramifica (&nodo, &lnodo, &rnodo, tsp0);
		nueva_U = false;
		if (Solucion(&rnodo)) {
			if (rnodo.ci < U) {    // se ha encontrado una solucion mejor
				U = rnodo.ci;
				nueva_U = true;
				CopiaNodo (&rnodo, &solucion);
			}
		}
		else {                    //  no es un nodo solucion
			if (rnodo.ci < U) {     //  cota inferior menor que cota superior
				if (!PilaPush (&pila, &rnodo)) {
					printf ("Error: pila agotada\n");
					liberarMatriz(tsp0);
					exit (1);
				}
			}
		}
		if (Solucion(&lnodo)) {
			if (lnodo.ci < U) {    // se ha encontrado una solucion mejor
				U = lnodo.ci;
				nueva_U = true;
				CopiaNodo (&lnodo,&solucion);
			}
		}
		else {                     // no es nodo solucion
			if (lnodo.ci < U) {      // cota inferior menor que cota superior
				if (!PilaPush (&pila, &lnodo)) {
					printf ("Error: pila agotada\n");
					liberarMatriz(tsp0);
					exit (1);
				}
			}
		}
		if (nueva_U) PilaAcotar (&pila, U);
		activo = PilaPop (&pila, &nodo);
	}
 
        MPI::Finalize();
	EscribeNodo(&solucion);
	liberarMatriz(tsp0);
}

Ejercicios propuestos

El alumno deberá implementar el algoritmo paralelo de Branch-and-Bound siguiendo las indicaciones de esta página. Deberá hacerse un estudio experimental con las siguientes medidas y cálculos relacionados:

(A) Incluyendo la difusión de la cota superior.
(B) Excluyendo la difusión de la cota superior.

Los archivos proporcionados son los siguientes:

libbb.h archivo de cabecera para las funciones del algoritmo de Branch-and-Bound, de manejo de la pila de nodos y gestión de memoria dinámica.
libbb.cc Implementación de las funciones declaradas en libbb.h.
bbseq.cc Implementación del algoritmo de Branch-and-Bound secuencial que puede servir como plantilla para implementar el algoritmo paralelo (bbpar.cc).
tspXX.X Archivos de ejemplo con problemas TSP. Se incluyen archivos de hasta 40 ciudades.
Creado por: Daniel Guerrero Martínez y Sergio Rodríguez Lumley 2010

Valid HTML 4.01 Transitional