rcm.c

OpenFVM-v1.1 open source cfd code

	{
	  mask[i] = 0;
	}
      xflags &= (~RCM_USE_MASK);
    }
  DEGREE (n, xflags, xadj, adj, mask, deg);

  for (i = 0; i < n; ++i)
    {
      if (mask[i])
	{
	  continue;
	}
      root = i + base;

      /*
       * Find a pseudo-peripheral node and calculate RCM.
       * Note that we reuse perm to also store the level
       * structure (perm+num is parameter ls of FNROOT()).
       */
      FNROOT (&root, xflags, xadj, adj, deg, &ccsize, mask, perm + num);
      RCM (root, xflags, xadj, adj, deg, mask, perm + num, &ccsize);

      num += ccsize;
      CHECK (num <= n);
      if (num >= n)
	{
	  break;
	}

    }

}



/**
 * @brief Generates rooted level structure.
 *
 * Generates rooted level structure.
 * Only those nodes for
 * which @a mask is nonzero will be considered.
 *
 * @param[in] root
 *            the node at which the level structure is to
 *            be rooted
 * @param[in] flags
 *            tailors behavior. Use @c 0 to not set any flag.
 *            @c ROOTLS() supports the @c #RCM_FORTRAN_INDICES flag.
 * @param[in] xadj
 *            pointers to the adjacency lists in @a adj
 * @param[in] adj
 *            adjacency lists
 * @param[in,out] mask
 *            nodes @a p with @a mask[@a p] not zero are
 *            considered masked out.
 *            @a mask is modified, but restored to the original
 *            values on output.
 * @param[out] ccsize
 *            size of the component rooted by @a root
 * @param[out] nlvl
 *            number of levels found
 * @param[out] last_lvl
 *            pointer to the beginning of the last level in @a ls
 * @param[out] ls
 *            nodes in the level sets, ordered by level sets.
 *            The last level is stored in @a ls[@a k], where
 *            @a last_lvl <= @a k < @a ccsize.
 *
 */

inline static void
ROOTLS (const INT root, const int flags,
	const INT * restrict xadj, const INT * restrict adj,
	signed char *restrict mask,
	INT * restrict ccsize, INT * restrict nlvl,
	INT * restrict last_lvl, INT * restrict ls)
{
  INT i, j;
  INT nbr, node;

  INT lvl_begin, lvl_end;	/* C style indices */

  if (flags & RCM_FORTRAN_INDICES)
    {
      --mask;
      --xadj;
      --adj;
    }

  ASSERT (!mask[root]);
  mask[root] = 1;		/* mask root */

  /* init first level */
  *nlvl = 0;
  ls[0] = root;
  lvl_begin = 0;
  lvl_end = 1;
  *ccsize = 1;

  /*
   * lvl_begin is the pointer to the beginning of the current
   * level, and lvl_end-1 points to the end of this level.
   */
  do
    {
      /*
       * Generate the next level by finding all the non-masked
       * neighbors of all the nodes in the current level.
       */
      for (i = lvl_begin; i < lvl_end; ++i)
	{
	  node = ls[i];
	  for (j = xadj[node]; j < xadj[node + 1]; ++j)
	    {
	      nbr = adj[j];
	      if (!mask[nbr])
		{
		  mask[nbr] = 1;
		  ls[*ccsize] = nbr;
		  ++(*ccsize);
		}
	    }
	}

      /*
       * Reset lvl_begin and lvl_end to iterate over the now
       * generated next level.
       */
      *last_lvl = lvl_begin;
      lvl_begin = lvl_end;
      lvl_end = *ccsize;
      ++(*nlvl);
    }
  while (lvl_end - lvl_begin > 0);

  /*
   * Reset mask to zero for the nodes in the level structure,
   * i.e. the nodes are again unmasked.
   */
  for (i = 0; i < *ccsize; ++i)
    {
      mask[ls[i]] = 0;
    }

}

/**
 * @brief Generic @c fnrootX() code (see @c fnrooti()).
 *
 * Generic @c fnrootX() code (see @c fnrooti()).
 * @n
 * <i>The following documentation is copied from @c fnrooti():</i>
 *
 * @copydoc fnrooti()
 */
void
FNROOT (INT * restrict root, const int flags,
	const INT * restrict xadj, const INT * restrict adj,
	const INT * restrict deg,
	INT * restrict ccsize, signed char *restrict mask, INT * restrict ls)
{
  INT j, new_nlvl, last_lvl, node, ndeg, mindeg;
  INT nlvl = 1;

  if (flags & RCM_FORTRAN_INDICES)
    {
      --deg;
    }

  nlvl = 1;

  while (1)
    {
      ROOTLS (*root, flags, xadj, adj,
	      mask, ccsize, &new_nlvl, &last_lvl, ls);
      /*
       * Return if the number of levels stays the same
       * or if we have reached the maximum level count.
       * The number of levels can not decrease.
       */
      CHECK (new_nlvl > 0);
      CHECK (new_nlvl >= nlvl);
      CHECK (new_nlvl <= *ccsize);
      if (new_nlvl == nlvl || new_nlvl == *ccsize)
	{
	  break;
	}
      nlvl = new_nlvl;
      /*
       * Pick a node with minimum degree from the last level.
       * last_lvl points to the start of the last level.
       */
      mindeg = *ccsize;
      for (j = last_lvl; j < *ccsize; ++j)
	{
	  node = ls[j];
	  ndeg = deg[node];
	  if (ndeg < mindeg)
	    {
	      *root = node;
	      mindeg = ndeg;
	    }
	}
    }
}


/**
 * @brief Sifts down a heap element, helper for @c HEAPSORT().
 *
 * Helper function for @c HEAPSORT().
 * Sifts down the value @a t in the heap between @a i and <i>n</i>-1.
 *
 */
static void
SIFT (const INT i, const INT n, const INT t, INT * restrict perm,
      const INT * restrict deg)
{
  INT j, w;
  j = i;
  w = i * 2 + 1;

  while (w < n)
    {
      if ((w + 1 < n) && (deg[perm[w]] < deg[perm[w + 1]]))
	{
	  ++w;
	}
      if (deg[t] < deg[perm[w]])
	{
	  /* we have to exchange/rotate and to go further down */
	  perm[j] = perm[w];
	  j = w;
	  w = j * 2 + 1;
	}
      else
	{
	  /* v has heap property */
	  break;
	}
    }
  perm[j] = t;
}


/**
 * @brief Heapsort for RCM; sorts after the degree of nodes.
 *
 * Heapsort function to sorts an array.
 * This version is tailored for the RCM algorithm and sort
 * the entries in @a perm after their degree (stored in @a deg).
 * The @c #RCM_INSERTION_SORT flag controls if @c HEAPSORT()
 * or @c INSERTION_SORT() is used for sorting.
 *
 * @param[in] sz        size of the array
 * @param[in,out] perm  data to be sorted
 * @param[in] deg       degrees of the nodes
 *
 */
inline static void
HEAPSORT (const INT sz, INT * restrict perm, const INT * restrict deg)
{
  INT n = sz;
  INT i, t;

  for (i = n / 2; i > 0;)
    {
      --i;
      SIFT (i, n, perm[i], perm, deg);
    }
  for (; n > 1;)
    {
      --n;
      t = perm[n];
      perm[n] = perm[0];
      SIFT (i, n, t, perm, deg);
    }
}

/**
 * @brief Linear insertion sort for RCM; sorts after the degree of nodes.
 *
 * Linear insertion sort function for RCM.
 * It sorts the entries in @a perm after their degree (stored in @a deg).
 * The @c #RCM_INSERTION_SORT flag controls if @c HEAPSORT()
 * or @c INSERTION_SORT() is used for sorting.
 *
 * @param[in] sz        size of the array
 * @param[in,out] perm  data to be sorted
 * @param[in] deg       degrees of the nodes
 *
 */
inline static void
INSERTION_SORT (const INT sz, INT * restrict perm, const INT * restrict deg)
{
  INT k, l;
  INT nbr;

  for (k = sz - 1; k > 0; --k)
    {
      nbr = perm[k - 1];
      for (l = k; l < sz && deg[perm[l]] < deg[nbr]; ++l)
	{
	  perm[l - 1] = perm[l];
	}
      perm[l - 1] = nbr;
    }
}



/**
 * @brief Generic @c rcmX() code (see @c rcmi()).
 *
 * Generic @c rcmX() code (see @c rcmi()).
 * @n
 * <i>The following documentation is copied from @c rcmi():</i>
 *
 * @copydoc rcmi()
 */
void
RCM (const INT root, const int flags,
     const INT * restrict xadj, const INT * restrict adj,
     const INT * restrict deg,
     signed char *restrict mask, INT * restrict perm, INT * restrict ccsize)
{
  INT i, j;
  INT nbr, node;

  INT lvl_begin, lvl_end;
  INT fnbr, lnbr;

  int base = (flags & RCM_FORTRAN_INDICES) ? 1 : 0;
  int sort = !(flags & RCM_NO_SORT);

  /*
   * Fill the first level set with the root node,
   * and initialize lvl_begin, lvl_end and lnbr
   * accordingly.
   * Adjust the arrays if Fortran style indices are used
   * for xadj, adj and perm.
   */

  if (flags & RCM_FORTRAN_INDICES)
    {
      --xadj;
      --adj;
      --mask;
      --perm;
      --deg;
    }

  perm[base] = root;
  ASSERT (!mask[root]);
  mask[root] = -1;

  lvl_begin = 0 + base;
  lvl_end = 1 + base;
  lnbr = lvl_end;

  /*
   * Iterate over all level sets. lvl_begin and lvl_end point to the
   * beginning and behind the end of the current level respectively.
   */
  do
    {
      /* iterate over all nodes in the current level */
      for (i = lvl_begin; i < lvl_end; ++i)
	{
	  node = perm[i];
	  CHECK (mask[node] == -1);
	  /*
	   * Find the non-masked neighbors of node.
	   * fnbr and lnbr keep track of where these nodes are
	   * copied to in perm.
	   */
	  fnbr = lnbr;
	  for (j = xadj[node]; j < xadj[node + 1]; ++j)
	    {
	      nbr = adj[j];
	      if (!mask[nbr])
		{
		  mask[nbr] = -1;
		  perm[lnbr] = nbr;
		  ++lnbr;
		}
	    }
	  if (lnbr - fnbr > 1 && sort)
	    {
	      if ((flags & RCM_INSERTION_SORT))
		{
		  INSERTION_SORT (lnbr - fnbr, perm + fnbr, deg);
		}
	      else
		{
		  HEAPSORT (lnbr - fnbr, perm + fnbr, deg);
		}
	    }
	}
      lvl_begin = lvl_end;
      lvl_end = lnbr;
    }
  while (lvl_end - lvl_begin > 0);

  ASSERT (lnbr > 0);		/* check that no overflow occured */

  /* we now have a Cuthill-McKee ordering in perm */

  if (flags & RCM_FORTRAN_INDICES)
    {
      ++perm;
      --lnbr;
    }
  *ccsize = lnbr;

  /* reverse perm if reversal is not explicitly unwanted */
  if (!(flags & RCM_NO_REVERSE))
    {
      j = lnbr - 1;
      for (i = 0; i < lnbr / 2; ++i, --j)
	{
	  node = perm[j];
	  perm[j] = perm[i];
	  perm[i] = node;
	}
    }
}


#endif /* RCM__DRIVER */