#ifndef __CSGRAPH_H__
#define __CSGRAPH_H__

#include <vector>

/*
 * Determine connected compoments of a compressed sparse graph.
 * Note:
 *   Output array flag must be preallocated
 */
template <class I>
I cs_graph_components(const I n_nod,
		      const I Ap[],
		      const I Aj[],
		            I flag[])
{
  // pos is a work array: list of nodes (rows) to process.
  std::vector<I> pos(n_nod,01);
  I n_comp = 0;
  I n_tot, n_pos, n_pos_new, n_pos0, n_new, n_stop;
  I icomp, ii, ir, ic;

  n_stop = n_nod;
  for (ir = 0; ir < n_nod; ir++) {
    flag[ir] = -1;
    if (Ap[ir+1] == Ap[ir]) {
      n_stop--;
      flag[ir] = -2;
    }
  }

  n_tot = 0;
  for (icomp = 0; icomp < n_nod; icomp++) {
    // Find seed.
    ii = 0;
    while ((flag[ii] >= 0) || (flag[ii] == -2)) {
      ii++;
      if (ii >= n_nod) {
	/* Sanity check, if this happens, the graph is corrupted. */
	return -1;
      }
    }

    flag[ii] = icomp;
    pos[0] = ii;
    n_pos0 = 0;
    n_pos_new = n_pos = 1;

    for (ii = 0; ii < n_nod; ii++) {
      n_new = 0;
      for (ir = n_pos0; ir < n_pos; ir++) {
	for (ic = Ap[pos[ir]]; ic < Ap[pos[ir]+1]; ic++) {
	  if (flag[Aj[ic]] == -1) {
	    flag[Aj[ic]] = icomp;
	    pos[n_pos_new] = Aj[ic];
	    n_pos_new++;
	    n_new++;
	  }
	}
      }
      n_pos0 = n_pos;
      n_pos = n_pos_new;
      if (n_new == 0) break;
    }
    n_tot += n_pos;

    if (n_tot == n_stop) {
      n_comp = icomp + 1;
      break;
    }
  }

  return n_comp;
}

#endif