Tutorial: Jacobi method in GraphLab

I have implemented today a parallel version of the Jacobi method, a very simple algorithm for solving
a system of linear equations in Graphlab. The whole process including debugging
took me less than an hour and a half.

Just to show that creating new Graphlab applications is easy, I give below the
resulting code, which is about 20 lines of code including debug traces..

In GraphLab you first have to define node data and edge data. In our case each node represents a row in the linear system of equations, and each edge is a non-zero entry in the linear relation
matrix.

struct node_data{
  double x_i; //value from previous iteration, placeholder for next value
  double b_i; //Entry i of the observtion vector b
  double A_ii; //Entry i of the main diagonal of matrix A
}

And the edge data is:

struct edge_data{
   double A_ij; //contains the A_ij value of the matrix A: i != j
}

/**
 * Functionality: The code solves the linear system Ax = b using
 * The Jacobi algorithm. (A is a square matrix). 
 * A assumed to be full column rank.  Algorithm is described
 * http://en.wikipedia.org/wiki/Jacobi_method
 * Written by Danny Bickson 
*/
/***
 * JACOBI UPDATE FUNCTION
 * Update rule is:
 * x_i = (b_i - \sum_j A_ij * x_j)/A_ii
 */
void jacobi_update_function(gl_types::iscope &scope,
                          gl_types::icallback &scheduler) {


  /* GET current vertex data */
  vertex_data& vdata = scope.vertex_data();
  graphlab::edge_list outedgeid = scope.out_edge_ids();

  //store last round x_i for convergence detection
  vdata.prev_mean = vdata.cur_mean;

  double x_i = vdata.b_i;
  double A_ii = vdata.A_ii;
  assert(A_ii != 0); //diagonal of square matrix can not be zero

  // go over each edge (non zero A_ij value)
  for(size_t i = 0; i < outedgeid.size(); ++i) {
      edge_data& out_edge = scope.edge_data(outedgeid[i]);
      const vertex_data & other = scope.neighbor_vertex_data(scope.target(outedgeid[i]));
      x_i -= out_edge.A_ij * other.x_i;
    }

    x_i /= A_ii;
    //store computed result in node data
    vdata.x_i = x_i;
}

And here is an example for running Jacobi:

<129|0>bickson@biggerbro:~/newgraphlab/graphlabapi/debug/demoapps/gabp$ ./gabp 1 mat3x3 --scheduler="round_robin(max_iterations=25)"
Loading mat3x3
Creating 6 edges...
symmetry: 0
.max iterations = 25
step = 1
max_iterations = 25
INFO:     asynchronous_engine.hpp(run:94): Worker 0 started.

INFO:     asynchronous_engine.hpp(run:94): Worker 1 started.

INFO:     asynchronous_engine.hpp(run:102): Worker 0 finished.

INFO:     asynchronous_engine.hpp(run:102): Worker 1 finished.

Jacobi finished in 0.0013
Assuming the linear system is Ax=y, and the correct solution is x*,Jacobi converged to an accuracy norm(x-x*) of 1.41598e-22 
Writing result to file: mat3x3.out
You can read the file in Matlab using the load_c_gl.m matlab script

 === REPORT FOR core() ===
[Numeric]
ncpus:          2
[Other]
affinities:     false
compile_flags:
engine: async
scheduler:      round_robin
schedyield:     true
scope:  edge

 === REPORT FOR engine() ===
[Numeric]
num_edges:              6
num_syncs:              2
num_vertices:           3
updatecount:            75
[Timings]
runtime:                0 s
[Other]
termination_reason:     task depletion (natural)
[Numeric]
updatecount_vector:             75      (count: 2, min: 37, max: 38, avg: 37.5)
updatecount_vector.values:              37,38,

If anyone wants to experiment with the code, you should download GraphLab, the Jacobi example is found on demoapps/gabp/ directory.