/* ========================================================================== */
  /* === UMF_singletons ======================================================= */
  /* ========================================================================== */
  
  /* -------------------------------------------------------------------------- */
  /* UMFPACK Copyright (c) Timothy A. Davis, CISE,                              */
  /* Univ. of Florida.  All Rights Reserved.  See ../Doc/License for License.   */
  /* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
  /* -------------------------------------------------------------------------- */
  
  /* Find and order the row and column singletons of a matrix A.  If there are
   * row and column singletons, the output is a row and column permutation such
   * that the matrix is in the following form:
   *
   *	x x x x x x x x x
   *	. x x x x x x x x
   *	. . x x x x x x x
   *	. . . x . . . . .
   *	. . . x x . . . .
   *	. . . x x s s s s
   *	. . . x x s s s s
   *	. . . x x s s s s
   *	. . . x x s s s s
   *
   * The above example has 3 column singletons (the first three columns and
   * their corresponding pivot rows) and 2 row singletons.  The singletons are
   * ordered first, because they have zero Markowitz cost.  The LU factorization
   * for these first five rows and columns is free - there is no work to do
   * (except to scale the pivot columns for the 2 row singletons), and no
   * fill-in occurs.  The remaining submatrix (4-by-4 in the above example)
   * has no rows or columns with degree one.  It may have empty rows or columns.
   *
   * This algorithm does not perform a full permutation to block triangular
   * form.  If there are one or more singletons, then the matrix can be
   * permuted to block triangular form, but UMFPACK does not perform the full
   * BTF permutation (see also "dmperm" in MATLAB, CSparse cs_dmperm,
   * and SuiteSparse/BTF).
   */
  
  #include "umf_internal.h"
  
  #ifndef NDEBUG
  
  /* ========================================================================== */
  /* === debug routines ======================================================= */
  /* ========================================================================== */
  
  /* Dump the singleton queue */
  
  PRIVATE void dump_singletons
  (
      Int head,		/* head of the queue */
      Int tail,		/* tail of the queue */
      Int Next [ ],	/* Next [i] is the next object after i */
      char *name,		/* "row" or "col" */
      Int Deg [ ],	/* Deg [i] is the degree of object i */
      Int n		/* objects are in the range 0 to n-1 */
  )
  {
      Int i, next, cnt ;
      DEBUG6 (("%s Singleton list: head "ID" tail "ID"
  ", name, head, tail)) ;
      i = head ;
      ASSERT (head >= EMPTY && head < n) ;
      ASSERT (tail >= EMPTY && tail < n) ;
      cnt = 0 ;
      while (i != EMPTY)
      {
  	DEBUG7 ((" "ID": "ID" deg: "ID"
  ", cnt, i, Deg [i])) ;
  	ASSERT (i >= 0 && i < n) ;
  	next = Next [i] ;
  	if (i == tail) ASSERT (next == EMPTY) ;
  	i = next ;
  	cnt++ ;
  	ASSERT (cnt <= n) ;
      }
  }
  
  PRIVATE void dump_mat
  (
      char *xname,
      char *yname,
      Int nx,
      Int ny,
      const Int Xp [ ],
      const Int Xi [ ],
      Int Xdeg [ ],
      Int Ydeg [ ]
  )
  {
      Int x, y, p, p1, p2, xdeg, do_xdeg, ydeg ;
      DEBUG6 (("
   ==== Dump %s mat:
  ", xname)) ;
      for (x = 0 ; x < nx ; x++)
      {
  	p1 = Xp [x] ;
  	p2 = Xp [x+1] ;
  	xdeg = Xdeg [x] ;
  	DEBUG6 (("Dump %s "ID" p1 "ID" p2 "ID" deg "ID"
  ",
  	    xname, x, p1, p2, xdeg)) ;
  	do_xdeg = (xdeg >= 0) ;
  	for (p = p1 ; p < p2 ; p++)
  	{
  	    y = Xi [p] ;
  	    DEBUG7 (("    %s "ID" deg: ", yname, y)) ;
  	    ASSERT (y >= 0 && y < ny) ;
  	    ydeg = Ydeg [y] ;
  	    DEBUG7 ((ID"
  ", ydeg)) ;
  	    if (do_xdeg && ydeg >= 0)
  	    {
  		xdeg-- ;
  	    }
  	}
  	ASSERT (IMPLIES (do_xdeg, xdeg == 0)) ;
      }
  }
  #endif
  
  /* ========================================================================== */
  /* === create_row_form ====================================================== */
  /* ========================================================================== */
  
  /* Create the row-form R of the column-form input matrix A.  This could be done
   * by UMF_transpose, except that Rdeg has already been computed.
   */
  
  PRIVATE void create_row_form
  (
      /* input, not modified: */
      Int n_row,		    /* A is n_row-by-n_col, nz = Ap [n_col] */
      Int n_col,
      const Int Ap [ ],	    /* Ap [0..n_col]: column pointers for A */
      const Int Ai [ ],	    /* Ai [0..nz-1]:  row indices for A */
      Int Rdeg [ ],	    /* Rdeg [0..n_row-1]: row degrees */
  
      /* output, not defined on input: */
      Int Rp [ ],		    /* Rp [0..n_row]: row pointers for R */
      Int Ri [ ],		    /* Ri [0..nz-1]:  column indices for R */
  
      /* workspace, not defined on input or output */
      Int W [ ]		    /* size n_row */
  )
  {
      Int row, col, p, p2 ;
  
      /* create the row pointers */
      Rp [0] = 0 ;
      W [0] = 0 ;
      for (row = 0 ; row < n_row ; row++)
      {
  	Rp [row+1] = Rp [row] + Rdeg [row] ;
  	W [row] = Rp [row] ;
      }
  
      /* create the indices for the row-form */
      for (col = 0 ; col < n_col ; col++)
      {
  	p2 = Ap [col+1] ;
  	for (p = Ap [col] ; p < p2 ; p++)
  	{
  	    Ri [W [Ai [p]]++] = col ;
  	}
      }
  }
  
  /* ========================================================================== */
  /* === order_singletons ===================================================== */
  /* ========================================================================== */
  
  PRIVATE int order_singletons	/* return new number of singletons */
  (
      Int k,	    /* the number of singletons so far */
      Int head,
      Int tail,
      Int Next [ ],
      Int Xdeg [ ], Int Xperm [ ], const Int Xp [ ], const Int Xi [ ],
      Int Ydeg [ ], Int Yperm [ ], const Int Yp [ ], const Int Yi [ ]
  #ifndef NDEBUG
      , char *xname, char *yname, Int nx, Int ny
  #endif
  )
  {
      Int xpivot, x, y, ypivot, p, p2, deg ;
  
  #ifndef NDEBUG
      Int i, k1 = k ;
      dump_singletons (head, tail, Next, xname, Xdeg, nx) ;
      dump_mat (xname, yname, nx, ny, Xp, Xi, Xdeg, Ydeg) ;
      dump_mat (yname, xname, ny, nx, Yp, Yi, Ydeg, Xdeg) ;
  #endif
  
      while (head != EMPTY)
      {
  	/* remove the singleton at the head of the queue */
  	xpivot = head ;
  	DEBUG1 (("------ Order %s singleton: "ID"
  ", xname, xpivot)) ;
  	head = Next [xpivot] ;
  	if (head == EMPTY) tail = EMPTY ;
  
  #ifndef NDEBUG
  	if (k % 100 == 0) dump_singletons (head, tail, Next, xname, Xdeg, nx) ;
  #endif
  
  	ASSERT (Xdeg [xpivot] >= 0) ;
  	if (Xdeg [xpivot] != 1)
  	{
  	    /* This row/column x is empty.  The matrix is singular.
  	     * x will be ordered last in Xperm. */
  	    DEBUG1 (("empty %s, after singletons removed
  ", xname)) ;
  	    continue ;
  	}
  
  	/* find the ypivot to match with this xpivot */
  #ifndef NDEBUG
  	/* there can only be one ypivot, since the degree of x is 1 */
  	deg = 0 ;
  	p2 = Xp [xpivot+1] ;
  	for (p = Xp [xpivot] ; p < p2 ; p++)
  	{
  	    y = Xi [p] ;
  	    DEBUG1 (("%s: "ID"
  ", yname, y)) ;
  	    if (Ydeg [y] >= 0)
  	    {
  		/* this is a live index in this xpivot vector */
  		deg++ ;
  	    }
  	}
  	ASSERT (deg == 1) ;
  #endif
  
  	ypivot = EMPTY ;
  	p2 = Xp [xpivot+1] ;
  	for (p = Xp [xpivot] ; p < p2 ; p++)
  	{
  	    y = Xi [p] ;
  	    DEBUG1 (("%s: "ID"
  ", yname, y)) ;
  	    if (Ydeg [y] >= 0)
  	    {
  		/* this is a live index in this xpivot vector */
  		ypivot = y ;
  		break ;
  	    }
  	}
  
  	DEBUG1 (("Pivot %s: "ID"
  ", yname, ypivot)) ;
  	ASSERT (ypivot != EMPTY) ;
  	DEBUG1 (("deg "ID"
  ", Ydeg [ypivot])) ;
  	ASSERT (Ydeg [ypivot] >= 0) ;
  
  	/* decrement the degrees after removing this singleton */
  	DEBUG1 (("p1 "ID"
  ", Yp [ypivot])) ;
  	DEBUG1 (("p2 "ID"
  ", Yp [ypivot+1])) ;
  	p2 = Yp [ypivot+1] ;
  	for (p = Yp [ypivot] ; p < p2 ; p++)
  	{
  	    x = Yi [p] ;
  	    DEBUG1 (("    %s: "ID" deg: "ID"
  ", xname, x, Xdeg [x])) ;
  	    if (Xdeg [x] < 0) continue ;
  	    ASSERT (Xdeg [x] > 0) ;
  	    if (x == xpivot) continue ;
  	    deg = --(Xdeg [x]) ;
  	    ASSERT (Xdeg [x] >= 0) ;
  	    if (deg == 1)
  	    {
  		/* this is a new singleton, put at the end of the queue */
  		Next [x] = EMPTY ;
  		if (head == EMPTY)
  		{
  		    head = x ;
  		}
  		else
  		{
  		    ASSERT (tail != EMPTY) ;
  		    Next [tail] = x ;
  		}
  		tail = x ;
  		DEBUG1 ((" New %s singleton:  "ID"
  ", xname, x)) ;
  #ifndef NDEBUG
  		if (k % 100 == 0)
  		{
  		    dump_singletons (head, tail, Next, xname, Xdeg, nx) ;
  		}
  #endif
  	    }
  	}
  
  	/* flag the xpivot and ypivot by FLIP'ing the degrees */
  	Xdeg [xpivot] = FLIP (1) ;
  	Ydeg [ypivot] = FLIP (Ydeg [ypivot]) ;
  
  	/* keep track of the pivot row and column */
  	Xperm [k] = xpivot ;
  	Yperm [k] = ypivot ;
  	k++ ;
  
  #ifndef NDEBUG
  	if (k % 1000 == 0)
  	{
  	    dump_mat (xname, yname, nx, ny, Xp, Xi, Xdeg, Ydeg) ;
  	    dump_mat (yname, xname, ny, nx, Yp, Yi, Ydeg, Xdeg) ;
  	}
  #endif
      }
  
  #ifndef NDEBUG
      DEBUGm4 (("%s singletons: k = "ID"
  ", xname, k)) ;
      for (i = k1 ; i < k ; i++)
      {
  	DEBUG1 (("  %s: "ID" %s: "ID"
  ", xname, Xperm [i], yname, Yperm [i])) ;
      }
      ASSERT (k > 0) ;
  #endif
  
      return (k) ;
  }
  
  /* ========================================================================== */
  /* === find_any_singletons ================================================== */
  /* ========================================================================== */
  
  PRIVATE Int find_any_singletons	    /* returns # of singletons found */
  (
      /* input, not modified: */
      Int n_row,
      Int n_col,
      const Int Ap [ ],	    /* size n_col+1 */
      const Int Ai [ ],	    /* size nz = Ap [n_col] */
  
      /* input, modified on output: */
      Int Cdeg [ ],	    /* size n_col */
      Int Rdeg [ ],	    /* size n_row */
  
      /* output, not defined on input: */
      Int Cperm [ ],	    /* size n_col */
      Int Rperm [ ],	    /* size n_row */
      Int *p_n1r,		    /* # of row singletons */
      Int *p_n1c,		    /* # of col singletons */
  
      /* workspace, not defined on input or output */
      Int Rp [ ],		    /* size n_row+1 */
      Int Ri [ ],		    /* size nz */
      Int W [ ],		    /* size n_row */
      Int Next [ ]	    /* size MAX (n_row, n_col) */
  )
  {
      Int n1, col, row, row_form, head, tail, n1r, n1c ;
  
      /* ---------------------------------------------------------------------- */
      /* eliminate column singletons */
      /* ---------------------------------------------------------------------- */
  
      n1 = 0 ;
      n1r = 0 ;
      n1c = 0 ;
      row_form = FALSE ;
  
      head = EMPTY ;
      tail = EMPTY ;
      for (col = n_col-1 ; col >= 0 ; col--)
      {
  	if (Cdeg [col] == 1)
  	{
  	    /* put the column singleton in the queue */
  	    if (head == EMPTY) tail = col ;
  	    Next [col] = head ;
  	    head = col ;
  	    DEBUG1 (("Column singleton: "ID"
  ", col)) ;
  	}
      }
  
      if (head != EMPTY)
      {
  
  	/* ------------------------------------------------------------------ */
  	/* create the row-form of A */
  	/* ------------------------------------------------------------------ */
  
  	create_row_form (n_row, n_col, Ap, Ai, Rdeg, Rp, Ri, W) ;
  	row_form = TRUE ;
  
  	/* ------------------------------------------------------------------ */
  	/* find and order the column singletons */
  	/* ------------------------------------------------------------------ */
  
  	n1 = order_singletons (0, head, tail, Next,
  		Cdeg, Cperm, Ap, Ai,
  		Rdeg, Rperm, Rp, Ri
  #ifndef NDEBUG
  		, "col", "row", n_col, n_row
  #endif
  		) ;
  	n1c = n1 ;
      }
  
      /* ---------------------------------------------------------------------- */
      /* eliminate row singletons */
      /* ---------------------------------------------------------------------- */
  
      head = EMPTY ;
      tail = EMPTY ;
      for (row = n_row-1 ; row >= 0 ; row--)
      {
  	if (Rdeg [row] == 1)
  	{
  	    /* put the row singleton in the queue */
  	    if (head == EMPTY) tail = row ;
  	    Next [row] = head ;
  	    head = row ;
  	    DEBUG1 (("Row singleton: "ID"
  ", row)) ;
  	}
      }
  
      if (head != EMPTY)
      {
  
  	/* ------------------------------------------------------------------ */
  	/* create the row-form of A, if not already created */
  	/* ------------------------------------------------------------------ */
  
  	if (!row_form)
  	{
  	    create_row_form (n_row, n_col, Ap, Ai, Rdeg, Rp, Ri, W) ;
  	}
  
  	/* ------------------------------------------------------------------ */
  	/* find and order the row singletons */
  	/* ------------------------------------------------------------------ */
  
  	n1 = order_singletons (n1, head, tail, Next,
  		Rdeg, Rperm, Rp, Ri,
  		Cdeg, Cperm, Ap, Ai
  #ifndef NDEBUG
  		, "row", "col", n_row, n_col
  #endif
  		) ;
  	n1r = n1 - n1c ;
      }
  
      DEBUG0 (("n1 "ID"
  ", n1)) ;
      *p_n1r = n1r ;
      *p_n1c = n1c ;
      return (n1) ;
  }
  
  /* ========================================================================== */
  /* === find_user_singletons ================================================= */
  /* ========================================================================== */
  
  PRIVATE Int find_user_singletons	/* returns # singletons found */
  (
      /* input, not modified: */
      Int n_row,
      Int n_col,
      const Int Ap [ ],	    /* size n_col+1 */
      const Int Ai [ ],	    /* size nz = Ap [n_col] */
      const Int Quser [ ],    /* size n_col if present */
  
      /* input, modified on output: */
      Int Cdeg [ ],	    /* size n_col */
      Int Rdeg [ ],	    /* size n_row */
  
      /* output, not defined on input */
      Int Cperm [ ],	    /* size n_col */
      Int Rperm [ ],	    /* size n_row */
      Int *p_n1r,		    /* # of row singletons */
      Int *p_n1c,		    /* # of col singletons */
  
      /* workspace, not defined on input or output */
      Int Rp [ ],		    /* size n_row+1 */
      Int Ri [ ],		    /* size nz */
      Int W [ ]		    /* size n_row */
  )
  {
      Int n1, col, row, p, p2, pivcol, pivrow, found, k, n1r, n1c ;
  
      n1 = 0 ;
      n1r = 0 ;
      n1c = 0 ;
      *p_n1r = 0 ;
      *p_n1c = 0 ;
  
      /* find singletons in the user column permutation, Quser */
      pivcol = Quser [0] ;
      found = (Cdeg [pivcol] == 1) ;
      DEBUG0 (("Is first col: "ID" a col singleton?: "ID"
  ", pivcol, found)) ;
      if (!found)
      {
  	/* the first column is not a column singleton, check for a row
  	 * singleton in the first column. */
  	for (p = Ap [pivcol] ; p < Ap [pivcol+1] ; p++)
  	{
  	    if (Rdeg [Ai [p]] == 1)
  	    {
  		DEBUG0 (("Row singleton in first col: "ID" row: "ID"
  ",
  		    pivcol, Ai [p])) ;
  		found = TRUE ;
  		break ;
  	    }
  	}
      }
  
      if (!found)
      {
  	/* no singletons in the leading part of A (:,Quser) */
  	return (0) ;
      }
  
      /* there is at least one row or column singleton.  Look for more. */
      create_row_form (n_row, n_col, Ap, Ai, Rdeg, Rp, Ri, W) ;
  
      n1 = 0 ;
  
      for (k = 0 ; k < n_col ; k++)
      {
  	pivcol = Quser [k] ;
  	pivrow = EMPTY ;
  
  	/* ------------------------------------------------------------------ */
  	/* check if col is a column singleton, or contains a row singleton */
  	/* ------------------------------------------------------------------ */
  
  	found = (Cdeg [pivcol] == 1) ;
  
  	if (found)
  	{
  
  	    /* -------------------------------------------------------------- */
  	    /* pivcol is a column singleton */
  	    /* -------------------------------------------------------------- */
  
  	    DEBUG0 (("Found a col singleton: k "ID" pivcol "ID"
  ", k, pivcol));
  
  	    /* find the pivrow to match with this pivcol */
  #ifndef NDEBUG
  	    /* there can only be one pivrow, since the degree of pivcol is 1 */
  	    {
  		Int deg = 0 ;
  		p2 = Ap [pivcol+1] ;
  		for (p = Ap [pivcol] ; p < p2 ; p++)
  		{
  		    row = Ai [p] ;
  		    DEBUG1 (("row: "ID"
  ", row)) ;
  		    if (Rdeg [row] >= 0)
  		    {
  			/* this is a live index in this column vector */
  			deg++ ;
  		    }
  		}
  		ASSERT (deg == 1) ;
  	    }
  #endif
  
  	    p2 = Ap [pivcol+1] ;
  	    for (p = Ap [pivcol] ; p < p2 ; p++)
  	    {
  		row = Ai [p] ;
  		DEBUG1 (("row: "ID"
  ", row)) ;
  		if (Rdeg [row] >= 0)
  		{
  		    /* this is a live index in this pivcol vector */
  		    pivrow = row ;
  		    break ;
  		}
  	    }
  
  	    DEBUG1 (("Pivot row: "ID"
  ", pivrow)) ;
  	    ASSERT (pivrow != EMPTY) ;
  	    DEBUG1 (("deg "ID"
  ", Rdeg [pivrow])) ;
  	    ASSERT (Rdeg [pivrow] >= 0) ;
  
  	    /* decrement the degrees after removing this col singleton */
  	    DEBUG1 (("p1 "ID"
  ", Rp [pivrow])) ;
  	    DEBUG1 (("p2 "ID"
  ", Rp [pivrow+1])) ;
  	    p2 = Rp [pivrow+1] ;
  	    for (p = Rp [pivrow] ; p < p2 ; p++)
  	    {
  		col = Ri [p] ;
  		DEBUG1 (("    col: "ID" deg: "ID"
  ", col, Cdeg [col])) ;
  		if (Cdeg [col] < 0) continue ;
  		ASSERT (Cdeg [col] > 0) ;
  		Cdeg [col]-- ;
  		ASSERT (Cdeg [col] >= 0) ;
  	    }
  
  	    /* flag the pivcol and pivrow by FLIP'ing the degrees */
  	    Cdeg [pivcol] = FLIP (1) ;
  	    Rdeg [pivrow] = FLIP (Rdeg [pivrow]) ;
  	    n1c++ ;
  
  	}
  	else
  	{
  
  	    /* -------------------------------------------------------------- */
  	    /* pivcol may contain a row singleton */
  	    /* -------------------------------------------------------------- */
  
  	    p2 = Ap [pivcol+1] ;
  	    for (p = Ap [pivcol] ; p < p2 ; p++)
  	    {
  		pivrow = Ai [p] ;
  		if (Rdeg [pivrow] == 1)
  		{
  		    DEBUG0 (("Row singleton in pivcol: "ID" row: "ID"
  ",
  			pivcol, pivrow)) ;
  		    found = TRUE ;
  		    break ;
  		}
  	    }
  
  	    if (!found)
  	    {
  		DEBUG0 (("End of user singletons
  ")) ;
  		break ;
  	    }
  
  #ifndef NDEBUG
  	    /* there can only be one pivrow, since the degree of pivcol is 1 */
  	    {
  		Int deg = 0 ;
  		p2 = Rp [pivrow+1] ;
  		for (p = Rp [pivrow] ; p < p2 ; p++)
  		{
  		    col = Ri [p] ;
  		    DEBUG1 (("col: "ID" cdeg::: "ID"
  ", col, Cdeg [col])) ;
  		    if (Cdeg [col] >= 0)
  		    {
  			/* this is a live index in this column vector */
  			ASSERT (col == pivcol) ;
  			deg++ ;
  		    }
  		}
  		ASSERT (deg == 1) ;
  	    }
  #endif
  
  	    DEBUG1 (("Pivot row: "ID"
  ", pivrow)) ;
  	    DEBUG1 (("pivcol deg "ID"
  ", Cdeg [pivcol])) ;
  	    ASSERT (Cdeg [pivcol] > 1) ;
  
  	    /* decrement the degrees after removing this row singleton */
  	    DEBUG1 (("p1 "ID"
  ", Ap [pivcol])) ;
  	    DEBUG1 (("p2 "ID"
  ", Ap [pivcol+1])) ;
  	    p2 = Ap [pivcol+1] ;
  	    for (p = Ap [pivcol] ; p < p2 ; p++)
  	    {
  		row = Ai [p] ;
  		DEBUG1 (("    row: "ID" deg: "ID"
  ", row, Rdeg [row])) ;
  		if (Rdeg [row] < 0) continue ;
  		ASSERT (Rdeg [row] > 0) ;
  		Rdeg [row]-- ;
  		ASSERT (Rdeg [row] >= 0) ;
  	    }
  
  	    /* flag the pivcol and pivrow by FLIP'ing the degrees */
  	    Cdeg [pivcol] = FLIP (Cdeg [pivcol]) ;
  	    Rdeg [pivrow] = FLIP (1) ;
  	    n1r++ ;
  	}
  
  	/* keep track of the pivot row and column */
  	Cperm [k] = pivcol ;
  	Rperm [k] = pivrow ;
  	n1++ ;
  
  #ifndef NDEBUG
  	dump_mat ("col", "row", n_col, n_row, Ap, Ai, Cdeg, Rdeg) ;
  	dump_mat ("row", "col", n_row, n_col, Rp, Ri, Rdeg, Cdeg) ;
  #endif
  
      }
  
      DEBUGm4 (("User singletons found: "ID"
  ", n1)) ;
      ASSERT (n1 > 0) ;
  
      *p_n1r = n1r ;
      *p_n1c = n1c ;
      return (n1) ;
  }
  
  /* ========================================================================== */
  /* === finish_permutation =================================================== */
  /* ========================================================================== */
  
  /* Complete the permutation for the pruned submatrix.  The singletons are
   * already ordered, but remove their flags.  Place rows/columns that are empty
   * in the pruned submatrix at the end of the output permutation.  This can only
   * occur if the matrix is singular.
   */
  
  PRIVATE Int finish_permutation
  (
      Int n1,
      Int nx,
      Int Xdeg [ ],
      const Int Xuser [ ],
      Int Xperm [ ],
      Int *p_max_deg
  )
  {
      Int nempty, x, deg, s, max_deg, k ;
      nempty = 0 ;
      s = n1 ;
      max_deg = 0 ;
      DEBUG0 (("n1 "ID" nempty "ID"
  ", n1, nempty)) ;
      for (k = 0 ; k < nx ; k++)
      {
  	x = (Xuser != (Int *) NULL) ? Xuser [k] : k ;
  	DEBUG0 (("finish perm k "ID" x "ID" nx "ID"
  ", k, x, nx)) ;
  	deg = Xdeg [x] ;
  	if (deg == 0)
  	{
  	    /* this row/col is empty in the pruned submatrix */
  	    ASSERT (s < nx - nempty) ;
  	    DEBUG0 (("empty k "ID"
  ", k)) ;
  	    nempty++ ;
  	    Xperm [nx - nempty] = x ;
  	}
  	else if (deg > 0)
  	{
  	    /* this row/col is nonempty in the pruned submatrix */
  	    ASSERT (s < nx - nempty) ;
  	    Xperm [s++] = x ;
  	    max_deg = MAX (max_deg, deg) ;
  	}
  	else
  	{
  	    /* This is a singleton row/column - it is already ordered.
  	     * Just clear the flag. */
  	    Xdeg [x] = FLIP (deg) ;
  	}
      }
      ASSERT (s == nx - nempty) ;
      *p_max_deg = max_deg ;
      return (nempty) ;
  }
  
  /* ========================================================================== */
  /* === UMF_singletons ======================================================= */
  /* ========================================================================== */
  
  GLOBAL Int UMF_singletons
  (
  
      /* input, not modified: */
      Int n_row,
      Int n_col,
      const Int Ap [ ],	    /* size n_col+1 */
      const Int Ai [ ],	    /* size nz = Ap [n_col] */
      const Int Quser [ ],    /* size n_col if present */
      Int strategy,	    /* strategy requested by user */
  
      /* output, not defined on input: */
      Int Cdeg [ ],	/* size n_col */
      Int Cperm [ ],	/* size n_col */
      Int Rdeg [ ],	/* size n_row */
      Int Rperm [ ],	/* size n_row */
      Int InvRperm [ ],	/* size n_row, the inverse of Rperm */
      Int *p_n1,		/* # of col and row singletons */
      Int *p_n1c,		/* # of col singletons */
      Int *p_n1r,		/* # of row singletons */
      Int *p_nempty_col,	/* # of empty columns in pruned submatrix */
      Int *p_nempty_row,	/* # of empty columns in pruned submatrix */
      Int *p_is_sym,	/* TRUE if pruned submatrix is square and has been
  			 * symmetrically permuted by Cperm and Rperm */
      Int *p_max_rdeg,	/* maximum Rdeg in pruned submatrix */
  
      /* workspace, not defined on input or output */
      Int Rp [ ],		/* size n_row+1 */
      Int Ri [ ],		/* size nz */
      Int W [ ],		/* size n_row */
      Int Next [ ]	/* size MAX (n_row, n_col) */
  )
  {
      Int n1, s, col, row, p, p1, p2, cdeg, last_row, is_sym, k,
  	nempty_row, nempty_col, max_cdeg, max_rdeg, n1c, n1r ;
  
      /* ---------------------------------------------------------------------- */
      /* initializations */
      /* ---------------------------------------------------------------------- */
  
  #ifndef NDEBUG
      UMF_dump_start ( ) ;
      DEBUGm4 (("Starting umf_singletons
  ")) ;
  #endif
  
      /* ---------------------------------------------------------------------- */
      /* scan the columns, check for errors and count row degrees */
      /* ---------------------------------------------------------------------- */
  
      if (Ap [0] != 0 || Ap [n_col] < 0)
      {
  	return (UMFPACK_ERROR_invalid_matrix) ;
      }
      for (row = 0 ; row < n_row ; row++)
      {
  	Rdeg [row] = 0 ;
      }
      for (col = 0 ; col < n_col ; col++)
      {
  	p1 = Ap [col] ;
  	p2 = Ap [col+1] ;
  	cdeg = p2 - p1 ;
  	if (cdeg < 0)
  	{
  	    return (UMFPACK_ERROR_invalid_matrix) ;
  	}
  	last_row = EMPTY ;
  	for (p = p1 ; p < p2 ; p++)
  	{
  	    row = Ai [p] ;
  	    if (row <= last_row || row >= n_row)
  	    {
  		return (UMFPACK_ERROR_invalid_matrix) ;
  	    }
  	    Rdeg [row]++ ;
  	    last_row = row ;
  	}
  	Cdeg [col] = cdeg ;
      }
  
      /* ---------------------------------------------------------------------- */
      /* find singletons */
      /* ---------------------------------------------------------------------- */
  
      if (Quser != (Int *) NULL)
      {
  	/* user has provided an input column ordering */
  	if (strategy == UMFPACK_STRATEGY_UNSYMMETRIC)
  	{
  	    /* look for singletons, but respect the user's input permutation */
  	    n1 = find_user_singletons (n_row, n_col, Ap, Ai, Quser,
  		    Cdeg, Rdeg, Cperm, Rperm, &n1r, &n1c, Rp, Ri, W) ;
  	}
  	else
  	{
  	    /* do not look for singletons if Quser given and strategy is
  	     * not unsymmetric */
  	    n1 = 0 ;
  	    n1r = 0 ;
  	    n1c = 0 ;
  	}
      }
      else
      {
  	/* look for singletons anywhere */
  	n1 = find_any_singletons (n_row, n_col, Ap, Ai,
  		Cdeg, Rdeg, Cperm, Rperm, &n1r, &n1c, Rp, Ri, W, Next) ;
      }
  
      /* ---------------------------------------------------------------------- */
      /* eliminate empty columns and complete the column permutation */
      /* ---------------------------------------------------------------------- */
  
      nempty_col = finish_permutation (n1, n_col, Cdeg, Quser, Cperm, &max_cdeg) ;
  
      /* ---------------------------------------------------------------------- */
      /* eliminate empty rows and complete the row permutation */
      /* ---------------------------------------------------------------------- */
  
      if (Quser != (Int *) NULL && strategy == UMFPACK_STRATEGY_SYMMETRIC)
      {
  	/* rows should be symmetrically permuted according to Quser */
  	ASSERT (n_row == n_col) ;
  	nempty_row = finish_permutation (n1, n_row, Rdeg, Quser, Rperm,
  	    &max_rdeg) ;
      }
      else
      {
  	/* rows should not be symmetrically permuted according to Quser */
  	nempty_row = finish_permutation (n1, n_row, Rdeg, (Int *) NULL, Rperm,
  	    &max_rdeg) ;
      }
  
      /* ---------------------------------------------------------------------- */
      /* compute the inverse of Rperm */
      /* ---------------------------------------------------------------------- */
  
      for (k = 0 ; k < n_row ; k++)
      {
  	ASSERT (Rperm [k] >= 0 && Rperm [k] < n_row) ;
  	InvRperm [Rperm [k]] = k ;
      }
  
      /* ---------------------------------------------------------------------- */
      /* see if pruned submatrix is square and has been symmetrically permuted */
      /* ---------------------------------------------------------------------- */
  
      if (n_row == n_col && nempty_row == nempty_col)
      {
  	/* is_sym is true if the submatrix is square, and
  	 * Rperm [n1..n_row-nempty_row-1] = Cperm [n1..n_col-nempty_col-1] */
  	is_sym = TRUE ;
  	for (s = n1 ; s < n_col - nempty_col ; s++)
  	{
  	    if (Cperm [s] != Rperm [s])
  	    {
  		is_sym = FALSE ;
  		break ;
  	    }
  	}
      }
      else
      {
  	is_sym = FALSE ;
      }
      DEBUGm4 (("Submatrix square and symmetrically permuted? "ID"
  ", is_sym)) ;
      DEBUGm4 (("singletons "ID" row "ID" col "ID"
  ", n1, n1r, n1c)) ;
      DEBUGm4 (("Empty cols "ID" rows "ID"
  ", nempty_col, nempty_row)) ;
      *p_n1 = n1 ;
      *p_n1r = n1r ;
      *p_n1c = n1c ;
      *p_is_sym = is_sym ;
      *p_nempty_col = nempty_col ;
      *p_nempty_row = nempty_row ;
      *p_max_rdeg = max_rdeg ;
      return (UMFPACK_OK) ;
  }