A master-slave implementation of kmeans algorithm using PVM3.

kmeans-master.c

#include <pthread.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include <stdio.h>
#include <sys/types.h>
#include <pvm3.h>

#define CANT_SLAVES 20

#define N 100000
#define K 100
#define TOLERANCE 0.00001
#define SEED 11

typedef struct {
  double x;
  double y;
} Point;

void generateOctaveScript(Point* data, Point* means, long* groupForPoint, int step) {
  long i, j, first;
  
  printf("# STEP %d\n", step);
  printf("function step%d()\n", step);
  for (i = 0; i < K; i++) {
    first = 0;
    printf("means(%d, :) = [%f %f];\n", i+1, means[i].x, means[i].y);
    for (j = 0; j < N; j++) {
      if (groupForPoint[j] == i) {
        if (first == 0){
          first++;
          printf("data%d = [", i+1);
        }
        printf("[%f %f] ; ", data[j].x, data[j].y);
      }
    }
    printf("];\n");
  }

  printf("plot(");
  for (i = 1; i <= K; i++) {
    printf("means(%d, 1), means(%d, 2), '.', data%d(:,1), data%d(:,2), '+'", i, i, i, i);
    if (i < K)
      printf(",");
  }
  printf(");\n");
  printf("endfunction\n\n");
}

int kmeans(Point* data, Point* means, long* groupForPoint) {
  int i;

  long cantPoints = N;
  long cantMeans = K;

  long index;

  double prevError, error = DBL_MAX;
  double error_aux;

  int res;
  int myTID = pvm_mytid();
  int slaves[CANT_SLAVES];
  
  if ((res = pvm_spawn("kmeans_slave",NULL,PvmTaskDefault,"",CANT_SLAVES,slaves)) < 1){
    printf("Master: pvm_spawn error\n");
    pvm_exit();
    exit(1);
  }

  long cantN = cantPoints/CANT_SLAVES;
  long cantK = cantMeans/CANT_SLAVES;
  for (i = 0; i< CANT_SLAVES; i++) {
    pvm_initsend(PvmDataDefault);
    pvm_pklong(&cantN, 1, 1);
    pvm_pklong(&cantK, 1, 1);
    pvm_pklong(&cantPoints, 1, 1);
    pvm_pklong(&cantMeans, 1, 1);
    pvm_pkbyte((char*)data, sizeof(Point)*cantPoints, 1);
    pvm_send(slaves[i], 1);
  }

  int iter = 0;
  do {
    prevError = error;
    error = 0;
    
    for (i=0; i < CANT_SLAVES; i++) {
      index = i*cantN;
      pvm_initsend(PvmDataDefault);
      pvm_pklong(&index, 1, 1);
      pvm_pkbyte((char*)means, sizeof(Point)*cantMeans, 1);
      pvm_send(slaves[i], 1);
    }

    for (i=0; i < CANT_SLAVES; i++) {
      pvm_recv(-1, -1);
      pvm_upkdouble(&error_aux, 1, 1);
      error += error_aux;

      pvm_upklong(&index, 1, 1);
      pvm_upklong(groupForPoint + index, cantN, 1);
    }

    for (i=0; i < CANT_SLAVES; i++) {
      index = i*cantK;
      pvm_initsend(PvmDataDefault);
      pvm_pklong(&index, 1, 1);
      pvm_pklong(groupForPoint, cantPoints, 1);
      pvm_send(slaves[i], 1);
    }

    for (i=0; i < CANT_SLAVES; i++) {
      pvm_recv(-1, -1);
      pvm_upklong(&index, 1, 1);
      pvm_upkbyte((char*)means + (sizeof(Point)*index), sizeof(Point)*cantK, 1);
    }

    iter++;

  } while (fabs(error - prevError) > TOLERANCE);

  for (i=0; i < CANT_SLAVES; i++) {
    pvm_kill(slaves[i]);
  }

  return iter;
}

int main(void) {

  long i;

  Point* data;
  Point* means;
  long* groupForPoint;

  data =  (Point*)malloc(N*sizeof(Point));
  means =  (Point*)malloc(K*sizeof(Point));
  groupForPoint = (long*)malloc(N*sizeof(long));

  srand(SEED);

  /* Cargamos puntos */
  for (i = 0; i < N; i++) {
    data[i].x = (double)rand()/RAND_MAX;
    data[i].y = (double)rand()/RAND_MAX;
  }

  /* Inicializamos las medias */
  for (i = 0; i < K; i++) {
    means[i].x = (double)rand()/RAND_MAX;
    means[i].y = (double)rand()/RAND_MAX;
  }

  int iterations = kmeans(data, means, groupForPoint);

  generateOctaveScript(data, means, groupForPoint, iterations);

  free(groupForPoint);
  free(data);
  free(means);

}

kmeans-slave.c

#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include <stdio.h>
#include <sys/types.h>
#include <pvm3.h>

typedef struct {
  double x;
  double y;
} Point;

double distance(Point p1, Point p2) {
  Point diff;
  diff.x = p1.x - p2.x;
  diff.y = p1.y - p2.y;
  return (double) sqrt(diff.x*diff.x + diff.y*diff.y);
}

int main(void) {
  long maxIndex, p, p_aux, k, n, closerGroup;
  double minDist, dist, error;

  int myTID = pvm_mytid();

  Point* data;
  Point* means;
  long* localGroupForPoint;
  long* groupForEveryPoint;

  Point pAdd;

  long index, cantNSlave, cantKSlave, cantNTotal, cantKTotal;


  pvm_recv(-1, -1);
  pvm_upklong(&cantNSlave, 1, 1);
  pvm_upklong(&cantKSlave, 1, 1);
  pvm_upklong(&cantNTotal, 1, 1);
  data =  (Point*)malloc(cantNTotal*sizeof(Point));
  memset(data, 0, cantNTotal);
  
  pvm_upklong(&cantKTotal, 1, 1);
  pvm_upkbyte((char*)data, sizeof(Point)*cantNTotal, 1);

  means =  (Point*)malloc(cantKTotal*sizeof(Point));
  groupForEveryPoint = (long*)malloc(cantNTotal*sizeof(long));
  localGroupForPoint = (long*)malloc(cantNSlave*sizeof(long));

  memset(means, 0, cantKTotal);
  memset(groupForEveryPoint, 0, cantNTotal);
  memset(localGroupForPoint, 0, cantNSlave);

  /*
  pvm_initsend(PvmDataDefault);
  pvm_pkbyte((char*)data + (sizeof(Point)*1), sizeof(Point)*(cantNTotal-1), 1);
  pvm_send(pvm_parent(), 1);
  */

  while(1) {
    pvm_recv(-1, -1);
    pvm_upklong(&index, 1, 1);
    pvm_upkbyte((char*)means, sizeof(Point)*cantKTotal, 1);

    error = 0;

    maxIndex = ((index + cantNSlave) < cantNTotal) ? index + cantNSlave : cantNTotal;
    for (p = index, p_aux = 0; p < maxIndex; p++, p_aux++) {
      minDist = DBL_MAX;
      closerGroup = 0;

      for (k = 0; k < cantKTotal; k++) {
        dist = distance(data[p], means[k]);
        if (dist < minDist) {
          minDist = dist;
          closerGroup = k;
        }
      }
      error += minDist;
      localGroupForPoint[p_aux] = closerGroup;
    }

    pvm_initsend(PvmDataDefault);
    pvm_pkdouble(&error, 1, 1);
    pvm_pklong(&index, 1, 1);
    pvm_pklong(localGroupForPoint, cantNSlave, 1);
    pvm_send(pvm_parent(), 1);

    pvm_recv(-1, -1);
    pvm_upklong(&index, 1, 1);
    pvm_upklong(groupForEveryPoint, cantNTotal, 1);
    
    maxIndex = ((index + cantKSlave) < cantKTotal) ? index + cantKSlave : cantKTotal;
    for (k = index; k < maxIndex; k++) {
      pAdd.x = 0; pAdd.y = 0;
      n = 0;
      for (p = 0; p < cantNTotal; p++) {
        if (groupForEveryPoint[p] == k) {
          pAdd.x += data[p].x;
          pAdd.y += data[p].y;
          n++; 
        }
      }

      if (n > 0) {
        means[k].x = pAdd.x / n;
        means[k].y = pAdd.y / n;
      }
    }

    pvm_initsend(PvmDataDefault);
    pvm_pklong(&index, 1, 1);
    pvm_pkbyte((char*)means + sizeof(Point)*index, sizeof(Point)*cantKSlave, 1);
    pvm_send(pvm_parent(), 1);

  }

}