/* ========================================================================= */
/* === AMD_2 =============================================================== */
/* ========================================================================= */

/* ------------------------------------------------------------------------- */
/* AMD, Copyright (c) Timothy A. Davis,					     */
/* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
/* email: DrTimothyAldenDavis@gmail.com                                      */
/* ------------------------------------------------------------------------- */

/* AMD_2:  performs the AMD ordering on a symmetric sparse matrix A, followed
 * by a postordering (via depth-first search) of the assembly tree using the
 * AMD_postorder routine.
 */

module amd_2;

nothrow @nogc extern(C):

import glob_opts;
import SuiteSparse_config;
import amd_internal;
import amd;
import amd_postorder;
import core.stdc.math;

/* ========================================================================= */
/* === clear_flag ========================================================== */
/* ========================================================================= */

static Int clear_flag (Int wflg, Int wbig, Int *W, Int n)
{
    Int x ;
    if (wflg < 2 || wflg >= wbig)
    {
	for (x = 0 ; x < n ; x++)
	{
	    if (W [x] != 0) W [x] = 1 ;
	}
	wflg = 2 ;
    }
    /*  at this point, W [0..n-1] < wflg holds */
    return (wflg) ;
}


/* ========================================================================= */
/* === AMD_2 =============================================================== */
/* ========================================================================= */


/* amd_2 is the primary AMD ordering routine.  It is not meant to be
 * user-callable because of its restrictive inputs and because it destroys
 * the user's input matrix.  It does not check its inputs for errors, either.
 * However, if you can work with these restrictions it can be faster than
 * amd_order and use less memory (assuming that you can create your own copy
 * of the matrix for AMD to destroy).  Refer to AMD/Source/amd_2.c for a
 * description of each parameter. */


/*
 * Given a representation of the nonzero pattern of a symmetric matrix, A,
 * (excluding the diagonal) perform an approximate minimum (UMFPACK/MA38-style)
 * degree ordering to compute a pivot order such that the introduction of
 * nonzeros (fill-in) in the Cholesky factors A = LL' is kept low.  At each
 * step, the pivot selected is the one with the minimum UMFAPACK/MA38-style
 * upper-bound on the external degree.  This routine can optionally perform
 * aggresive absorption (as done by MC47B in the Harwell Subroutine
 * Library).
 *
 * The approximate degree algorithm implemented here is the symmetric analog of
 * the degree update algorithm in MA38 and UMFPACK (the Unsymmetric-pattern
 * MultiFrontal PACKage, both by Davis and Duff).  The routine is based on the
 * MA27 minimum degree ordering algorithm by Iain Duff and John Reid.
 *
 * This routine is a translation of the original AMDBAR and MC47B routines,
 * in Fortran, with the following modifications:
 *
 * (1) dense rows/columns are removed prior to ordering the matrix, and placed
 *	last in the output order.  The presence of a dense row/column can
 *	increase the ordering time by up to O(n^2), unless they are removed
 *	prior to ordering.
 *
 * (2) the minimum degree ordering is followed by a postordering (depth-first
 *	search) of the assembly tree.  Note that mass elimination (discussed
 *	below) combined with the approximate degree update can lead to the mass
 *	elimination of nodes with lower exact degree than the current pivot
 *	element.  No additional fill-in is caused in the representation of the
 *	Schur complement.  The mass-eliminated nodes merge with the current
 *	pivot element.  They are ordered prior to the current pivot element.
 *	Because they can have lower exact degree than the current element, the
 *	merger of two or more of these nodes in the current pivot element can
 *	lead to a single element that is not a "fundamental supernode".  The
 *	diagonal block can have zeros in it.  Thus, the assembly tree used here
 *	is not guaranteed to be the precise supernodal elemination tree (with
 *	"funadmental" supernodes), and the postordering performed by this
 *	routine is not guaranteed to be a precise postordering of the
 *	elimination tree.
 *
 * (3) input parameters are added, to control aggressive absorption and the
 *	detection of "dense" rows/columns of A.
 *
 * (4) additional statistical information is returned, such as the number of
 *	nonzeros in L, and the flop counts for subsequent LDL' and LU
 *	factorizations.  These are slight upper bounds, because of the mass
 *	elimination issue discussed above.
 *
 * (5) additional routines are added to interface this routine to MATLAB
 *	to provide a simple C-callable user-interface, to check inputs for
 *	errors, compute the symmetry of the pattern of A and the number of
 *	nonzeros in each row/column of A+A', to compute the pattern of A+A',
 *	to perform the assembly tree postordering, and to provide debugging
 *	ouput.  Many of these functions are also provided by the Fortran
 *	Harwell Subroutine Library routine MC47A.
 *
 * (6) both int and SuiteSparse_long versions are provided.  In the
 *      descriptions below and integer is and int or SuiteSparse_long depending
 *      on which version is being used.

 **********************************************************************
 ***** CAUTION:  ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT.  ******
 **********************************************************************
 ** If you want error checking, a more versatile input format, and a **
 ** simpler user interface, use amd_order or amd_l_order instead.    **
 ** This routine is not meant to be user-callable.                   **
 **********************************************************************

 * ----------------------------------------------------------------------------
 * References:
 * ----------------------------------------------------------------------------
 *
 *  [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern multifrontal
 *	method for sparse LU factorization", SIAM J. Matrix Analysis and
 *	Applications, vol. 18, no. 1, pp. 140-158.  Discusses UMFPACK / MA38,
 *	which first introduced the approximate minimum degree used by this
 *	routine.
 *
 *  [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An approximate
 *	minimum degree ordering algorithm," SIAM J. Matrix Analysis and
 *	Applications, vol. 17, no. 4, pp. 886-905, 1996.  Discusses AMDBAR and
 *	MC47B, which are the Fortran versions of this routine.
 *
 *  [3] Alan George and Joseph Liu, "The evolution of the minimum degree
 *	ordering algorithm," SIAM Review, vol. 31, no. 1, pp. 1-19, 1989.
 *	We list below the features mentioned in that paper that this code
 *	includes:
 *
 *	mass elimination:
 *	    Yes.  MA27 relied on supervariable detection for mass elimination.
 *
 *	indistinguishable nodes:
 *	    Yes (we call these "supervariables").  This was also in the MA27
 *	    code - although we modified the method of detecting them (the
 *	    previous hash was the true degree, which we no longer keep track
 *	    of).  A supervariable is a set of rows with identical nonzero
 *	    pattern.  All variables in a supervariable are eliminated together.
 *	    Each supervariable has as its numerical name that of one of its
 *	    variables (its principal variable).
 *
 *	quotient graph representation:
 *	    Yes.  We use the term "element" for the cliques formed during
 *	    elimination.  This was also in the MA27 code.  The algorithm can
 *	    operate in place, but it will work more efficiently if given some
 *	    "elbow room."
 *
 *	element absorption:
 *	    Yes.  This was also in the MA27 code.
 *
 *	external degree:
 *	    Yes.  The MA27 code was based on the true degree.
 *
 *	incomplete degree update and multiple elimination:
 *	    No.  This was not in MA27, either.  Our method of degree update
 *	    within MC47B is element-based, not variable-based.  It is thus
 *	    not well-suited for use with incomplete degree update or multiple
 *	    elimination.
 *
 * Authors, and Copyright (C) 2004 by:
 * Timothy A. Davis, Patrick Amestoy, Iain S. Duff, John K. Reid.
 *
 * Acknowledgements: This work (and the UMFPACK package) was supported by the
 * National Science Foundation (ASC-9111263, DMS-9223088, and CCR-0203270).
 * The UMFPACK/MA38 approximate degree update algorithm, the unsymmetric analog
 * which forms the basis of AMD, was developed while Tim Davis was supported by
 * CERFACS (Toulouse, France) in a post-doctoral position.  This C version, and
 * the etree postorder, were written while Tim Davis was on sabbatical at
 * Stanford University and Lawrence Berkeley National Laboratory.

 * ----------------------------------------------------------------------------
 * INPUT ARGUMENTS (unaltered):
 * ----------------------------------------------------------------------------

 * n:  The matrix order.  Restriction:  n >= 1.
 *
 * iwlen:  The size of the Iw array.  On input, the matrix is stored in
 *	Iw [0..pfree-1].  However, Iw [0..iwlen-1] should be slightly larger
 *	than what is required to hold the matrix, at least iwlen >= pfree + n.
 *	Otherwise, excessive compressions will take place.  The recommended
 *	value of iwlen is 1.2 * pfree + n, which is the value used in the
 *	user-callable interface to this routine (amd_order.c).  The algorithm
 *	will not run at all if iwlen < pfree.  Restriction: iwlen >= pfree + n.
 *	Note that this is slightly more restrictive than the actual minimum
 *	(iwlen >= pfree), but AMD_2 will be very slow with no elbow room.
 *	Thus, this routine enforces a bare minimum elbow room of size n.
 *
 * pfree: On input the tail end of the array, Iw [pfree..iwlen-1], is empty,
 *	and the matrix is stored in Iw [0..pfree-1].  During execution,
 *	additional data is placed in Iw, and pfree is modified so that
 *	Iw [pfree..iwlen-1] is always the unused part of Iw.
 *
 * Control:  A c_float array of size AMD_CONTROL containing input parameters
 *	that affect how the ordering is computed.  If NULL, then default
 *	settings are used.
 *
 *	Control [AMD_DENSE] is used to determine whether or not a given input
 *	row is "dense".  A row is "dense" if the number of entries in the row
 *	exceeds Control [AMD_DENSE] times sqrt (n), except that rows with 16 or
 *	fewer entries are never considered "dense".  To turn off the detection
 *	of dense rows, set Control [AMD_DENSE] to a negative number, or to a
 *	number larger than sqrt (n).  The default value of Control [AMD_DENSE]
 *	is AMD_DEFAULT_DENSE, which is defined in amd.h as 10.
 *
 *	Control [AMD_AGGRESSIVE] is used to determine whether or not aggressive
 *	absorption is to be performed.  If nonzero, then aggressive absorption
 *	is performed (this is the default).

 * ----------------------------------------------------------------------------
 * INPUT/OUPUT ARGUMENTS:
 * ----------------------------------------------------------------------------
 *
 * Pe:  An integer array of size n.  On input, Pe [i] is the index in Iw of
 *	the start of row i.  Pe [i] is ignored if row i has no off-diagonal
 *	entries.  Thus Pe [i] must be in the range 0 to pfree-1 for non-empty
 *	rows.
 *
 *	During execution, it is used for both supervariables and elements:
 *
 *	Principal supervariable i:  index into Iw of the description of
 *	    supervariable i.  A supervariable represents one or more rows of
 *	    the matrix with identical nonzero pattern.  In this case,
 *	    Pe [i] >= 0.
 *
 *	Non-principal supervariable i:  if i has been absorbed into another
 *	    supervariable j, then Pe [i] = FLIP (j), where FLIP (j) is defined
 *	    as (-(j)-2).  Row j has the same pattern as row i.  Note that j
 *	    might later be absorbed into another supervariable j2, in which
 *	    case Pe [i] is still FLIP (j), and Pe [j] = FLIP (j2) which is
 *	    < EMPTY, where EMPTY is defined as (-1) in amd_internal.h.
 *
 *	Unabsorbed element e:  the index into Iw of the description of element
 *	    e, if e has not yet been absorbed by a subsequent element.  Element
 *	    e is created when the supervariable of the same name is selected as
 *	    the pivot.  In this case, Pe [i] >= 0.
 *
 *	Absorbed element e:  if element e is absorbed into element e2, then
 *	    Pe [e] = FLIP (e2).  This occurs when the pattern of e (which we
 *	    refer to as Le) is found to be a subset of the pattern of e2 (that
 *	    is, Le2).  In this case, Pe [i] < EMPTY.  If element e is "null"
 *	    (it has no nonzeros outside its pivot block), then Pe [e] = EMPTY,
 *	    and e is the root of an assembly subtree (or the whole tree if
 *	    there is just one such root).
 *
 *	Dense variable i:  if i is "dense", then Pe [i] = EMPTY.
 *
 *	On output, Pe holds the assembly tree/forest, which implicitly
 *	represents a pivot order with identical fill-in as the actual order
 *	(via a depth-first search of the tree), as follows.  If Nv [i] > 0,
 *	then i represents a node in the assembly tree, and the parent of i is
 *	Pe [i], or EMPTY if i is a root.  If Nv [i] = 0, then (i, Pe [i])
 *	represents an edge in a subtree, the root of which is a node in the
 *	assembly tree.  Note that i refers to a row/column in the original
 *	matrix, not the permuted matrix.
 *
 * Info:  A c_float array of size AMD_INFO.  If present, (that is, not NULL),
 *	then statistics about the ordering are returned in the Info array.
 *	See amd.h for a description.

 * ----------------------------------------------------------------------------
 * INPUT/MODIFIED (undefined on output):
 * ----------------------------------------------------------------------------
 *
 * Len:  An integer array of size n.  On input, Len [i] holds the number of
 *	entries in row i of the matrix, excluding the diagonal.  The contents
 *	of Len are undefined on output.
 *
 * Iw:  An integer array of size iwlen.  On input, Iw [0..pfree-1] holds the
 *	description of each row i in the matrix.  The matrix must be symmetric,
 *	and both upper and lower triangular parts must be present.  The
 *	diagonal must not be present.  Row i is held as follows:
 *
 *	    Len [i]:  the length of the row i data structure in the Iw array.
 *	    Iw [Pe [i] ... Pe [i] + Len [i] - 1]:
 *		the list of column indices for nonzeros in row i (simple
 *		supervariables), excluding the diagonal.  All supervariables
 *		start with one row/column each (supervariable i is just row i).
 *		If Len [i] is zero on input, then Pe [i] is ignored on input.
 *
 *	    Note that the rows need not be in any particular order, and there
 *	    may be empty space between the rows.
 *
 *	During execution, the supervariable i experiences fill-in.  This is
 *	represented by placing in i a list of the elements that cause fill-in
 *	in supervariable i:
 *
 *	    Len [i]:  the length of supervariable i in the Iw array.
 *	    Iw [Pe [i] ... Pe [i] + Elen [i] - 1]:
 *		the list of elements that contain i.  This list is kept short
 *		by removing absorbed elements.
 *	    Iw [Pe [i] + Elen [i] ... Pe [i] + Len [i] - 1]:
 *		the list of supervariables in i.  This list is kept short by
 *		removing nonprincipal variables, and any entry j that is also
 *		contained in at least one of the elements (j in Le) in the list
 *		for i (e in row i).
 *
 *	When supervariable i is selected as pivot, we create an element e of
 *	the same name (e=i):
 *
 *	    Len [e]:  the length of element e in the Iw array.
 *	    Iw [Pe [e] ... Pe [e] + Len [e] - 1]:
 *		the list of supervariables in element e.
 *
 *	An element represents the fill-in that occurs when supervariable i is
 *	selected as pivot (which represents the selection of row i and all
 *	non-principal variables whose principal variable is i).  We use the
 *	term Le to denote the set of all supervariables in element e.  Absorbed
 *	supervariables and elements are pruned from these lists when
 *	computationally convenient.
 *
 *  CAUTION:  THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION.
 *  The contents of Iw are undefined on output.

 * ----------------------------------------------------------------------------
 * OUTPUT (need not be set on input):
 * ----------------------------------------------------------------------------
 *
 * Nv:  An integer array of size n.  During execution, ABS (Nv [i]) is equal to
 *	the number of rows that are represented by the principal supervariable
 *	i.  If i is a nonprincipal or dense variable, then Nv [i] = 0.
 *	Initially, Nv [i] = 1 for all i.  Nv [i] < 0 signifies that i is a
 *	principal variable in the pattern Lme of the current pivot element me.
 *	After element me is constructed, Nv [i] is set back to a positive
 *	value.
 *
 *	On output, Nv [i] holds the number of pivots represented by super
 *	row/column i of the original matrix, or Nv [i] = 0 for non-principal
 *	rows/columns.  Note that i refers to a row/column in the original
 *	matrix, not the permuted matrix.
 *
 * Elen:  An integer array of size n.  See the description of Iw above.  At the
 *	start of execution, Elen [i] is set to zero for all rows i.  During
 *	execution, Elen [i] is the number of elements in the list for
 *	supervariable i.  When e becomes an element, Elen [e] = FLIP (esize) is
 *	set, where esize is the size of the element (the number of pivots, plus
 *	the number of nonpivotal entries).  Thus Elen [e] < EMPTY.
 *	Elen (i) = EMPTY set when variable i becomes nonprincipal.
 *
 *	For variables, Elen (i) >= EMPTY holds until just before the
 *	postordering and permutation vectors are computed.  For elements,
 *	Elen [e] < EMPTY holds.
 *
 *	On output, Elen [i] is the degree of the row/column in the Cholesky
 *	factorization of the permuted matrix, corresponding to the original row
 *	i, if i is a super row/column.  It is equal to EMPTY if i is
 *	non-principal.  Note that i refers to a row/column in the original
 *	matrix, not the permuted matrix.
 *
 *	Note that the contents of Elen on output differ from the Fortran
 *	version (Elen holds the inverse permutation in the Fortran version,
 *	which is instead returned in the Next array in this C version,
 *	described below).
 *
 * Last: In a degree list, Last [i] is the supervariable preceding i, or EMPTY
 *	if i is the head of the list.  In a hash bucket, Last [i] is the hash
 *	key for i.
 *
 *	Last [Head [hash]] is also used as the head of a hash bucket if
 *	Head [hash] contains a degree list (see the description of Head,
 *	below).
 *
 *	On output, Last [0..n-1] holds the permutation.  That is, if
 *	i = Last [k], then row i is the kth pivot row (where k ranges from 0 to
 *	n-1).  Row Last [k] of A is the kth row in the permuted matrix, PAP'.
 *
 * Next: Next [i] is the supervariable following i in a link list, or EMPTY if
 *	i is the last in the list.  Used for two kinds of lists:  degree lists
 *	and hash buckets (a supervariable can be in only one kind of list at a
 *	time).
 *
 *	On output Next [0..n-1] holds the inverse permutation. 	That is, if
 *	k = Next [i], then row i is the kth pivot row. Row i of A appears as
 *	the (Next[i])-th row in the permuted matrix, PAP'.
 *
 *	Note that the contents of Next on output differ from the Fortran
 *	version (Next is undefined on output in the Fortran version).

 * ----------------------------------------------------------------------------
 * LOCAL WORKSPACE (not input or output - used only during execution):
 * ----------------------------------------------------------------------------
 *
 * Degree:  An integer array of size n.  If i is a supervariable, then
 *	Degree [i] holds the current approximation of the external degree of
 *	row i (an upper bound).  The external degree is the number of nonzeros
 *	in row i, minus ABS (Nv [i]), the diagonal part.  The bound is equal to
 *	the exact external degree if Elen [i] is less than or equal to two.
 *
 *	We also use the term "external degree" for elements e to refer to
 *	|Le \ Lme|.  If e is an element, then Degree [e] is |Le|, which is the
 *	degree of the off-diagonal part of the element e (not including the
 *	diagonal part).
 *
 * Head:   An integer array of size n.  Head is used for degree lists.
 *	Head [deg] is the first supervariable in a degree list.  All
 *	supervariables i in a degree list Head [deg] have the same approximate
 *	degree, namely, deg = Degree [i].  If the list Head [deg] is empty then
 *	Head [deg] = EMPTY.
 *
 *	During supervariable detection Head [hash] also serves as a pointer to
 *	a hash bucket.  If Head [hash] >= 0, there is a degree list of degree
 *	hash.  The hash bucket head pointer is Last [Head [hash]].  If
 *	Head [hash] = EMPTY, then the degree list and hash bucket are both
 *	empty.  If Head [hash] < EMPTY, then the degree list is empty, and
 *	FLIP (Head [hash]) is the head of the hash bucket.  After supervariable
 *	detection is complete, all hash buckets are empty, and the
 *	(Last [Head [hash]] = EMPTY) condition is restored for the non-empty
 *	degree lists.
 *
 * W:  An integer array of size n.  The flag array W determines the status of
 *	elements and variables, and the external degree of elements.
 *
 *	for elements:
 *	    if W [e] = 0, then the element e is absorbed.
 *	    if W [e] >= wflg, then W [e] - wflg is the size of the set
 *		|Le \ Lme|, in terms of nonzeros (the sum of ABS (Nv [i]) for
 *		each principal variable i that is both in the pattern of
 *		element e and NOT in the pattern of the current pivot element,
 *		me).
 *	    if wflg > W [e] > 0, then e is not absorbed and has not yet been
 *		seen in the scan of the element lists in the computation of
 *		|Le\Lme| in Scan 1 below.
 *
 *	for variables:
 *	    during supervariable detection, if W [j] != wflg then j is
 *	    not in the pattern of variable i.
 *
 *	The W array is initialized by setting W [i] = 1 for all i, and by
 *	setting wflg = 2.  It is reinitialized if wflg becomes too large (to
 *	ensure that wflg+n does not cause integer overflow).

 * ----------------------------------------------------------------------------
 * LOCAL INTEGERS:
 * ----------------------------------------------------------------------------
 */


version (DLONG){
    void AMD_2
    (
        SuiteSparse_long n,		/* A is n-by-n, where n > 0 */
        SuiteSparse_long *Pe,		/* Pe [0..n-1]: index in Iw of row i on input */
        SuiteSparse_long *Iw,		/* workspace of size iwlen. Iw [0..pfree-1] holds the matrix on input */
        SuiteSparse_long *Len,	/* Len [0..n-1]: length for row/column i on input */
        SuiteSparse_long iwlen,		/* length of Iw. iwlen >= pfree + n */
        SuiteSparse_long pfree,		/* Iw [pfree ... iwlen-1] is empty on input */
        /* 7 size-n workspaces, not defined on input: */
		SuiteSparse_long *Nv,		/* the size of each supernode on output */
        SuiteSparse_long *Next,	/* the output inverse permutation */
        SuiteSparse_long *Last,	/* the output permutation */
        SuiteSparse_long *Head,
        SuiteSparse_long *Elen,/* the size columns of L for each supernode */
        SuiteSparse_long *Degree,
        SuiteSparse_long *W,
		/* control parameters and output statistics */
        c_float *Control,	/* array of size AMD_CONTROL */
        c_float *Info	/* array of size AMD_INFO */
    )
	{

    Int deg;
	Int degme;
	Int dext;
	Int lemax;
	Int e;
	Int elenme;
	Int eln;
	Int i;
	Int ilast;
	Int inext;
	Int j;
	Int jlast;
	Int jnext;
	Int k;
	Int knt1;
	Int knt2;
	Int knt3;
	Int lenj;
	Int ln;
	Int me;
	Int mindeg;
	Int nel;
	Int nleft;
	Int nvi;
	Int nvj;
	Int nvpiv;
	Int slenme;
	Int wbig;
	Int we;
	Int wflg;
	Int wnvi;
	Int ok;
	Int ndense;
	Int ncmpa;
	Int dense;
	Int aggressive;

    //unsigned Int hash ;	    /* unsigned, so that hash % n is well defined.*/
	UInt hash ;	    /* unsigned, so that hash % n is well defined.*/

/*
 * deg:		the degree of a variable or element
 * degme:	size, |Lme|, of the current element, me (= Degree [me])
 * dext:	external degree, |Le \ Lme|, of some element e
 * lemax:	largest |Le| seen so far (called dmax in Fortran version)
 * e:		an element
 * elenme:	the length, Elen [me], of element list of pivotal variable
 * eln:		the length, Elen [...], of an element list
 * hash:	the computed value of the hash function
 * i:		a supervariable
 * ilast:	the entry in a link list preceding i
 * inext:	the entry in a link list following i
 * j:		a supervariable
 * jlast:	the entry in a link list preceding j
 * jnext:	the entry in a link list, or path, following j
 * k:		the pivot order of an element or variable
 * knt1:	loop counter used during element construction
 * knt2:	loop counter used during element construction
 * knt3:	loop counter used during compression
 * lenj:	Len [j]
 * ln:		length of a supervariable list
 * me:		current supervariable being eliminated, and the current
 *		    element created by eliminating that supervariable
 * mindeg:	current minimum degree
 * nel:		number of pivots selected so far
 * nleft:	n - nel, the number of nonpivotal rows/columns remaining
 * nvi:		the number of variables in a supervariable i (= Nv [i])
 * nvj:		the number of variables in a supervariable j (= Nv [j])
 * nvpiv:	number of pivots in current element
 * slenme:	number of variables in variable list of pivotal variable
 * wbig:	= (INT_MAX - n) for the int version, (SuiteSparse_long_max - n)
 *                  for the SuiteSparse_long version.  wflg is not allowed to
 *                  be >= wbig.
 * we:		W [e]
 * wflg:	used for flagging the W array.  See description of Iw.
 * wnvi:	wflg - Nv [i]
 * x:		either a supervariable or an element
 *
 * ok:		true if supervariable j can be absorbed into i
 * ndense:	number of "dense" rows/columns
 * dense:	rows/columns with initial degree > dense are considered "dense"
 * aggressive:	true if aggressive absorption is being performed
 * ncmpa:	number of garbage collections

 * ----------------------------------------------------------------------------
 * LOCAL DOUBLES, used for statistical output only (except for alpha):
 * ----------------------------------------------------------------------------
 */

    c_float f, r, ndiv, s, nms_lu, nms_ldl, dmax, alpha, lnz, lnzme ;

/*
 * f:		nvpiv
 * r:		degme + nvpiv
 * ndiv:	number of divisions for LU or LDL' factorizations
 * s:		number of multiply-subtract pairs for LU factorization, for the
 *		    current element me
 * nms_lu	number of multiply-subtract pairs for LU factorization
 * nms_ldl	number of multiply-subtract pairs for LDL' factorization
 * dmax:	the largest number of entries in any column of L, including the
 *		    diagonal
 * alpha:	"dense" degree ratio
 * lnz:		the number of nonzeros in L (excluding the diagonal)
 * lnzme:	the number of nonzeros in L (excl. the diagonal) for the
 *		    current element me

 * ----------------------------------------------------------------------------
 * LOCAL "POINTERS" (indices into the Iw array)
 * ----------------------------------------------------------------------------
*/

    Int p, p1, p2, p3, p4, pdst, pend, pj, pme, pme1, pme2, pn, psrc ;

/*
 * Any parameter (Pe [...] or pfree) or local variable starting with "p" (for
 * Pointer) is an index into Iw, and all indices into Iw use variables starting
 * with "p."  The only exception to this rule is the iwlen input argument.
 *
 * p:           pointer into lots of things
 * p1:          Pe [i] for some variable i (start of element list)
 * p2:          Pe [i] + Elen [i] -  1 for some variable i
 * p3:          index of first supervariable in clean list
 * p4:
 * pdst:        destination pointer, for compression
 * pend:        end of memory to compress
 * pj:          pointer into an element or variable
 * pme:         pointer into the current element (pme1...pme2)
 * pme1:        the current element, me, is stored in Iw [pme1...pme2]
 * pme2:        the end of the current element
 * pn:          pointer into a "clean" variable, also used to compress
 * psrc:        source pointer, for compression
*/

/* ========================================================================= */
/*  INITIALIZATIONS */
/* ========================================================================= */

    /* Note that this restriction on iwlen is slightly more restrictive than
     * what is actually required in AMD_2.  AMD_2 can operate with no elbow
     * room at all, but it will be slow.  For better performance, at least
     * size-n elbow room is enforced. */
    ASSERT (iwlen >= pfree + n) ;
    ASSERT (n > 0) ;

    /* initialize output statistics */
    lnz = 0 ;
    ndiv = 0 ;
    nms_lu = 0 ;
    nms_ldl = 0 ;
    dmax = 1 ;
    me = EMPTY ;

    mindeg = 0 ;
    ncmpa = 0 ;
    nel = 0 ;
    lemax = 0 ;

    /* get control parameters */
    if (Control != cast(c_float *) NULL)
    {
	alpha = Control [AMD_DENSE] ;
	aggressive = (Control [AMD_AGGRESSIVE] != 0) ;
    }
    else
    {
	alpha = AMD_DEFAULT_DENSE ;
	aggressive = AMD_DEFAULT_AGGRESSIVE ;
    }
    /* Note: if alpha is NaN, this is undefined: */
    if (alpha < 0)
    {
	/* only remove completely dense rows/columns */
	dense = n-2 ;
    }
    else
    {
	dense = cast(Int) (alpha * sqrt (cast(c_float) n)) ;
    }
    dense = MAX (16, dense) ;
    dense = MIN (n,  dense) ;
    AMD_DEBUG1 ("\n\nAMD (debug), alpha %g, aggr. "~ID~"\n", alpha, aggressive);

    for (i = 0 ; i < n ; i++)
    {
	Last [i] = EMPTY ;
	Head [i] = EMPTY ;
	Next [i] = EMPTY ;
	/* if separate Hhead array is used for hash buckets: *
	Hhead [i] = EMPTY ;
	*/
	Nv [i] = 1 ;
	W [i] = 1 ;
	Elen [i] = 0 ;
	Degree [i] = Len [i] ;
    }

version(NDEBUG){}
else {
    AMD_DEBUG1 ("\n======Nel "~ID~" initial\n", nel);
    AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, Last,
		Head, Elen, Degree, W, -1) ;
} // !NDEBUG

    /* initialize wflg */
    wbig = Int_MAX - n ;
    wflg = clear_flag (0, wbig, W, n) ;

    /* --------------------------------------------------------------------- */
    /* initialize degree lists and eliminate dense and empty rows */
    /* --------------------------------------------------------------------- */

    ndense = 0 ;

    for (i = 0 ; i < n ; i++)
    {
	deg = Degree [i] ;
	ASSERT (deg >= 0 && deg < n) ;
	if (deg == 0)
	{

	    /* -------------------------------------------------------------
	     * we have a variable that can be eliminated at once because
	     * there is no off-diagonal non-zero in its row.  Note that
	     * Nv [i] = 1 for an empty variable i.  It is treated just
	     * the same as an eliminated element i.
	     * ------------------------------------------------------------- */

	    Elen [i] = FLIP (1) ;
	    nel++ ;
	    Pe [i] = EMPTY ;
	    W [i] = 0 ;

	}
	else if (deg > dense)
	{

	    /* -------------------------------------------------------------
	     * Dense variables are not treated as elements, but as unordered,
	     * non-principal variables that have no parent.  They do not take
	     * part in the postorder, since Nv [i] = 0.  Note that the Fortran
	     * version does not have this option.
	     * ------------------------------------------------------------- */

	    AMD_DEBUG1 ("Dense node "~ID~" degree "~ID~"\n", i, deg);
	    ndense++ ;
	    Nv [i] = 0 ;		/* do not postorder this node */
	    Elen [i] = EMPTY ;
	    nel++ ;
	    Pe [i] = EMPTY ;

	}
	else
	{

	    /* -------------------------------------------------------------
	     * place i in the degree list corresponding to its degree
	     * ------------------------------------------------------------- */

	    inext = Head [deg] ;
	    ASSERT (inext >= EMPTY && inext < n) ;
	    if (inext != EMPTY) Last [inext] = i ;
	    Next [i] = inext ;
	    Head [deg] = i ;

	}
    }

/* ========================================================================= */
/* WHILE (selecting pivots) DO */
/* ========================================================================= */

    while (nel < n)
    {

version(NDEBUG){}
else {
	AMD_DEBUG1 ("\n======Nel "~ID~"\n", nel) ;
	//if (AMD_debug >= 2)
	if (true)		// todo : 
	{
	    AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next,
		    Last, Head, Elen, Degree, W, nel) ;
	}
} // !NDEBUG

/* ========================================================================= */
/* GET PIVOT OF MINIMUM DEGREE */
/* ========================================================================= */

	/* ----------------------------------------------------------------- */
	/* find next supervariable for elimination */
	/* ----------------------------------------------------------------- */

	ASSERT (mindeg >= 0 && mindeg < n) ;
	for (deg = mindeg ; deg < n ; deg++)
	{
	    me = Head [deg] ;
	    if (me != EMPTY) break ;
	}
	mindeg = deg ;
	ASSERT (me >= 0 && me < n) ;
	AMD_DEBUG1 ("=================me: "~ID~"\n", me) ;

	/* ----------------------------------------------------------------- */
	/* remove chosen variable from link list */
	/* ----------------------------------------------------------------- */

	inext = Next [me] ;
	ASSERT (inext >= EMPTY && inext < n) ;
	if (inext != EMPTY) Last [inext] = EMPTY ;
	Head [deg] = inext ;

	/* ----------------------------------------------------------------- */
	/* me represents the elimination of pivots nel to nel+Nv[me]-1. */
	/* place me itself as the first in this set. */
	/* ----------------------------------------------------------------- */

	elenme = Elen [me] ;
	nvpiv = Nv [me] ;
	ASSERT (nvpiv > 0) ;
	nel += nvpiv ;

/* ========================================================================= */
/* CONSTRUCT NEW ELEMENT */
/* ========================================================================= */

	/* -----------------------------------------------------------------
	 * At this point, me is the pivotal supervariable.  It will be
	 * converted into the current element.  Scan list of the pivotal
	 * supervariable, me, setting tree pointers and constructing new list
	 * of supervariables for the new element, me.  p is a pointer to the
	 * current position in the old list.
	 * ----------------------------------------------------------------- */

	/* flag the variable "me" as being in Lme by negating Nv [me] */
	Nv [me] = -nvpiv ;
	degme = 0 ;
	ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;

	if (elenme == 0)
	{

	    /* ------------------------------------------------------------- */
	    /* construct the new element in place */
	    /* ------------------------------------------------------------- */

	    pme1 = Pe [me] ;
	    pme2 = pme1 - 1 ;

	    for (p = pme1 ; p <= pme1 + Len [me] - 1 ; p++)
	    {
		i = Iw [p] ;
		ASSERT (i >= 0 && i < n && Nv [i] >= 0) ;
		nvi = Nv [i] ;
		if (nvi > 0)
		{

		    /* ----------------------------------------------------- */
		    /* i is a principal variable not yet placed in Lme. */
		    /* store i in new list */
		    /* ----------------------------------------------------- */

		    /* flag i as being in Lme by negating Nv [i] */
		    degme += nvi ;
		    Nv [i] = -nvi ;
		    Iw [++pme2] = i ;

		    /* ----------------------------------------------------- */
		    /* remove variable i from degree list. */
		    /* ----------------------------------------------------- */

		    ilast = Last [i] ;
		    inext = Next [i] ;
		    ASSERT (ilast >= EMPTY && ilast < n) ;
		    ASSERT (inext >= EMPTY && inext < n) ;
		    if (inext != EMPTY) Last [inext] = ilast ;
		    if (ilast != EMPTY)
		    {
			Next [ilast] = inext ;
		    }
		    else
		    {
			/* i is at the head of the degree list */
			ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
			Head [Degree [i]] = inext ;
		    }
		}
	    }
	}
	else
	{

	    /* ------------------------------------------------------------- */
	    /* construct the new element in empty space, Iw [pfree ...] */
	    /* ------------------------------------------------------------- */

	    p = Pe [me] ;
	    pme1 = pfree ;
	    slenme = Len [me] - elenme ;

	    for (knt1 = 1 ; knt1 <= elenme + 1 ; knt1++)
	    {

		if (knt1 > elenme)
		{
		    /* search the supervariables in me. */
		    e = me ;
		    pj = p ;
		    ln = slenme ;
		    AMD_DEBUG2 ("Search sv: "~ID~" "~ID~" "~ID~"\n", me,pj,ln);
		}
		else
		{
		    /* search the elements in me. */
		    e = Iw [p++] ;
		    ASSERT (e >= 0 && e < n) ;
		    pj = Pe [e] ;
		    ln = Len [e] ;
		    AMD_DEBUG2 ("Search element e "~ID~" in me "~ID~"\n", e,me);
		    ASSERT (Elen [e] < EMPTY && W [e] > 0 && pj >= 0) ;
		}
		ASSERT (ln >= 0 && (ln == 0 || (pj >= 0 && pj < iwlen))) ;

		/* ---------------------------------------------------------
		 * search for different supervariables and add them to the
		 * new list, compressing when necessary. this loop is
		 * executed once for each element in the list and once for
		 * all the supervariables in the list.
		 * --------------------------------------------------------- */

		for (knt2 = 1 ; knt2 <= ln ; knt2++)
		{
		    i = Iw [pj++] ;
		    ASSERT (i >= 0 && i < n && (i == me || Elen [i] >= EMPTY));
		    nvi = Nv [i] ;
		    AMD_DEBUG2 (": "~ID~" "~ID~" "~ID~" "~ID~"\n", i, Elen [i], Nv [i], wflg);

		    if (nvi > 0)
		    {

			/* ------------------------------------------------- */
			/* compress Iw, if necessary */
			/* ------------------------------------------------- */

			if (pfree >= iwlen)
			{

			    AMD_DEBUG1 ("GARBAGE COLLECTION\n");

			    /* prepare for compressing Iw by adjusting pointers
			     * and lengths so that the lists being searched in
			     * the inner and outer loops contain only the
			     * remaining entries. */

			    Pe [me] = p ;
			    Len [me] -= knt1 ;
			    /* check if nothing left of supervariable me */
			    if (Len [me] == 0) Pe [me] = EMPTY ;
			    Pe [e] = pj ;
			    Len [e] = ln - knt2 ;
			    /* nothing left of element e */
			    if (Len [e] == 0) Pe [e] = EMPTY ;

			    ncmpa++ ;	/* one more garbage collection */

			    /* store first entry of each object in Pe */
			    /* FLIP the first entry in each object */
			    for (j = 0 ; j < n ; j++)
			    {
				pn = Pe [j] ;
				if (pn >= 0)
				{
				    ASSERT (pn >= 0 && pn < iwlen) ;
				    Pe [j] = Iw [pn] ;
				    Iw [pn] = FLIP (j) ;
				}
			    }

			    /* psrc/pdst point to source/destination */
			    psrc = 0 ;
			    pdst = 0 ;
			    pend = pme1 - 1 ;

			    while (psrc <= pend)
			    {
				/* search for next FLIP'd entry */
				j = FLIP (Iw [psrc++]) ;
				if (j >= 0)
				{
				    AMD_DEBUG2 ("Got object j: "~ID~"\n", j);
				    Iw [pdst] = Pe [j] ;
				    Pe [j] = pdst++ ;
				    lenj = Len [j] ;
				    /* copy from source to destination */
				    for (knt3 = 0 ; knt3 <= lenj - 2 ; knt3++)
				    {
					Iw [pdst++] = Iw [psrc++] ;
				    }
				}
			    }

			    /* move the new partially-constructed element */
			    p1 = pdst ;
			    for (psrc = pme1 ; psrc <= pfree-1 ; psrc++)
			    {
				Iw [pdst++] = Iw [psrc] ;
			    }
			    pme1 = p1 ;
			    pfree = pdst ;
			    pj = Pe [e] ;
			    p = Pe [me] ;

			}

			/* ------------------------------------------------- */
			/* i is a principal variable not yet placed in Lme */
			/* store i in new list */
			/* ------------------------------------------------- */

			/* flag i as being in Lme by negating Nv [i] */
			degme += nvi ;
			Nv [i] = -nvi ;
			Iw [pfree++] = i ;
			AMD_DEBUG2 ("     s: "~ID~"     nv "~ID~"\n", i, Nv [i]);

			/* ------------------------------------------------- */
			/* remove variable i from degree link list */
			/* ------------------------------------------------- */

			ilast = Last [i] ;
			inext = Next [i] ;
			ASSERT (ilast >= EMPTY && ilast < n) ;
			ASSERT (inext >= EMPTY && inext < n) ;
			if (inext != EMPTY) Last [inext] = ilast ;
			if (ilast != EMPTY)
			{
			    Next [ilast] = inext ;
			}
			else
			{
			    /* i is at the head of the degree list */
			    ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
			    Head [Degree [i]] = inext ;
			}
		    }
		}

		if (e != me)
		{
		    /* set tree pointer and flag to indicate element e is
		     * absorbed into new element me (the parent of e is me) */
		    AMD_DEBUG1 (" Element "~ID~" => "~ID~"\n", e, me);
		    Pe [e] = FLIP (me) ;
		    W [e] = 0 ;
		}
	    }

	    pme2 = pfree - 1 ;
	}

	/* ----------------------------------------------------------------- */
	/* me has now been converted into an element in Iw [pme1..pme2] */
	/* ----------------------------------------------------------------- */

	/* degme holds the external degree of new element */
	Degree [me] = degme ;
	Pe [me] = pme1 ;
	Len [me] = pme2 - pme1 + 1 ;
	ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;

	Elen [me] = FLIP (nvpiv + degme) ;
	/* FLIP (Elen (me)) is now the degree of pivot (including
	 * diagonal part). */

version(NDEBUG){}
else {
	AMD_DEBUG2 ("New element structure: length= "~ID~"\n", pme2-pme1+1);
	for (pme = pme1 ; pme <= pme2 ; pme++) AMD_DEBUG3 (" "~ID~"", Iw[pme]);
	AMD_DEBUG3 ("\n");
} // !NDEBUG

	/* ----------------------------------------------------------------- */
	/* make sure that wflg is not too large. */
	/* ----------------------------------------------------------------- */

	/* With the current value of wflg, wflg+n must not cause integer
	 * overflow */

	wflg = clear_flag (wflg, wbig, W, n) ;

/* ========================================================================= */
/* COMPUTE (W [e] - wflg) = |Le\Lme| FOR ALL ELEMENTS */
/* ========================================================================= */

	/* -----------------------------------------------------------------
	 * Scan 1:  compute the external degrees of previous elements with
	 * respect to the current element.  That is:
	 *       (W [e] - wflg) = |Le \ Lme|
	 * for each element e that appears in any supervariable in Lme.  The
	 * notation Le refers to the pattern (list of supervariables) of a
	 * previous element e, where e is not yet absorbed, stored in
	 * Iw [Pe [e] + 1 ... Pe [e] + Len [e]].  The notation Lme
	 * refers to the pattern of the current element (stored in
	 * Iw [pme1..pme2]).   If aggressive absorption is enabled, and
	 * (W [e] - wflg) becomes zero, then the element e will be absorbed
	 * in Scan 2.
	 * ----------------------------------------------------------------- */

	AMD_DEBUG2 ("me: ");
	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n) ;
	    eln = Elen [i] ;
	    AMD_DEBUG3 (""~ID~" Elen "~ID~": \n", i, eln);
	    if (eln > 0)
	    {
		/* note that Nv [i] has been negated to denote i in Lme: */
		nvi = -Nv [i] ;
		ASSERT (nvi > 0 && Pe [i] >= 0 && Pe [i] < iwlen) ;
		wnvi = wflg - nvi ;
		for (p = Pe [i] ; p <= Pe [i] + eln - 1 ; p++)
		{
		    e = Iw [p] ;
		    ASSERT (e >= 0 && e < n) ;
		    we = W [e] ;
		    AMD_DEBUG4 ("    e "~ID~" we "~ID~" ", e, we);
		    if (we >= wflg)
		    {
			/* unabsorbed element e has been seen in this loop */
			AMD_DEBUG4 ("    unabsorbed, first time seen");
			we -= nvi ;
		    }
		    else if (we != 0)
		    {
			/* e is an unabsorbed element */
			/* this is the first we have seen e in all of Scan 1 */
			AMD_DEBUG4 ("    unabsorbed");
			we = Degree [e] + wnvi ;
		    }
		    AMD_DEBUG4 ("\n");
		    W [e] = we ;
		}
	    }
	}
	AMD_DEBUG2 ("\n");

/* ========================================================================= */
/* DEGREE UPDATE AND ELEMENT ABSORPTION */
/* ========================================================================= */

	/* -----------------------------------------------------------------
	 * Scan 2:  for each i in Lme, sum up the degree of Lme (which is
	 * degme), plus the sum of the external degrees of each Le for the
	 * elements e appearing within i, plus the supervariables in i.
	 * Place i in hash list.
	 * ----------------------------------------------------------------- */

	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n && Nv [i] < 0 && Elen [i] >= 0) ;
	    AMD_DEBUG2 ("Updating: i "~ID~" "~ID~" "~ID~"\n", i, Elen[i], Len [i]);
	    p1 = Pe [i] ;
	    p2 = p1 + Elen [i] - 1 ;
	    pn = p1 ;
	    hash = 0 ;
	    deg = 0 ;
	    ASSERT (p1 >= 0 && p1 < iwlen && p2 >= -1 && p2 < iwlen) ;

	    /* ------------------------------------------------------------- */
	    /* scan the element list associated with supervariable i */
	    /* ------------------------------------------------------------- */

	    /* UMFPACK/MA38-style approximate degree: */
	    if (aggressive)
	    {
		for (p = p1 ; p <= p2 ; p++)
		{
		    e = Iw [p] ;
		    ASSERT (e >= 0 && e < n) ;
		    we = W [e] ;
		    if (we != 0)
		    {
			/* e is an unabsorbed element */
			/* dext = | Le \ Lme | */
			dext = we - wflg ;
			if (dext > 0)
			{
			    deg += dext ;
			    Iw [pn++] = e ;
			    hash += e ;
			    AMD_DEBUG4 (" e: "~ID~" hash = "~ID~"\n",e,hash);
			}
			else
			{
			    /* external degree of e is zero, absorb e into me*/
			    AMD_DEBUG1 (" Element "~ID~" =>"~ID~" (aggressive)\n", e, me);
			    ASSERT (dext == 0) ;
			    Pe [e] = FLIP (me) ;
			    W [e] = 0 ;
			}
		    }
		}
	    }
	    else
	    {
		for (p = p1 ; p <= p2 ; p++)
		{
		    e = Iw [p] ;
		    ASSERT (e >= 0 && e < n) ;
		    we = W [e] ;
		    if (we != 0)
		    {
			/* e is an unabsorbed element */
			dext = we - wflg ;
			ASSERT (dext >= 0) ;
			deg += dext ;
			Iw [pn++] = e ;
			hash += e ;
			AMD_DEBUG4 ("	e: "~ID~" hash = "~ID~"\n",e,hash);
		    }
		}
	    }

	    /* count the number of elements in i (including me): */
	    Elen [i] = pn - p1 + 1 ;

	    /* ------------------------------------------------------------- */
	    /* scan the supervariables in the list associated with i */
	    /* ------------------------------------------------------------- */

	    /* The bulk of the AMD run time is typically spent in this loop,
	     * particularly if the matrix has many dense rows that are not
	     * removed prior to ordering. */
	    p3 = pn ;
	    p4 = p1 + Len [i] ;
	    for (p = p2 + 1 ; p < p4 ; p++)
	    {
		j = Iw [p] ;
		ASSERT (j >= 0 && j < n) ;
		nvj = Nv [j] ;
		if (nvj > 0)
		{
		    /* j is unabsorbed, and not in Lme. */
		    /* add to degree and add to new list */
		    deg += nvj ;
		    Iw [pn++] = j ;
		    hash += j ;
		    AMD_DEBUG4 ("  s: "~ID~" hash "~ID~" Nv[j]= "~ID~"\n", j, hash, nvj) ;
		}
	    }

	    /* ------------------------------------------------------------- */
	    /* update the degree and check for mass elimination */
	    /* ------------------------------------------------------------- */

	    /* with aggressive absorption, deg==0 is identical to the
	     * Elen [i] == 1 && p3 == pn test, below. */
	    ASSERT (IMPLIES (aggressive, (deg==0) == (Elen[i]==1 && p3==pn))) ;

	    if (Elen [i] == 1 && p3 == pn)
	    {

		/* --------------------------------------------------------- */
		/* mass elimination */
		/* --------------------------------------------------------- */

		/* There is nothing left of this node except for an edge to
		 * the current pivot element.  Elen [i] is 1, and there are
		 * no variables adjacent to node i.  Absorb i into the
		 * current pivot element, me.  Note that if there are two or
		 * more mass eliminations, fillin due to mass elimination is
		 * possible within the nvpiv-by-nvpiv pivot block.  It is this
		 * step that causes AMD's analysis to be an upper bound.
		 *
		 * The reason is that the selected pivot has a lower
		 * approximate degree than the true degree of the two mass
		 * eliminated nodes.  There is no edge between the two mass
		 * eliminated nodes.  They are merged with the current pivot
		 * anyway.
		 *
		 * No fillin occurs in the Schur complement, in any case,
		 * and this effect does not decrease the quality of the
		 * ordering itself, just the quality of the nonzero and
		 * flop count analysis.  It also means that the post-ordering
		 * is not an exact elimination tree post-ordering. */

		AMD_DEBUG1 ("  MASS i "~ID~" => parent e "~ID~"\n", i, me);
		Pe [i] = FLIP (me) ;
		nvi = -Nv [i] ;
		degme -= nvi ;
		nvpiv += nvi ;
		nel += nvi ;
		Nv [i] = 0 ;
		Elen [i] = EMPTY ;

	    }
	    else
	    {

		/* --------------------------------------------------------- */
		/* update the upper-bound degree of i */
		/* --------------------------------------------------------- */

		/* the following degree does not yet include the size
		 * of the current element, which is added later: */

		Degree [i] = MIN (Degree [i], deg) ;

		/* --------------------------------------------------------- */
		/* add me to the list for i */
		/* --------------------------------------------------------- */

		/* move first supervariable to end of list */
		Iw [pn] = Iw [p3] ;
		/* move first element to end of element part of list */
		Iw [p3] = Iw [p1] ;
		/* add new element, me, to front of list. */
		Iw [p1] = me ;
		/* store the new length of the list in Len [i] */
		Len [i] = pn - p1 + 1 ;

		/* --------------------------------------------------------- */
		/* place in hash bucket.  Save hash key of i in Last [i]. */
		/* --------------------------------------------------------- */

		/* NOTE: this can fail if hash is negative, because the ANSI C
		 * standard does not define a % b when a and/or b are negative.
		 * That's why hash is defined as an unsigned Int, to avoid this
		 * problem. */
		hash = hash % n ;
		ASSERT ((cast(Int) hash) >= 0 && (cast(Int) hash) < n) ;

		/* if the Hhead array is not used: */
		j = Head [hash] ;
		if (j <= EMPTY)
		{
		    /* degree list is empty, hash head is FLIP (j) */
		    Next [i] = FLIP (j) ;
		    Head [hash] = FLIP (i) ;
		}
		else
		{
		    /* degree list is not empty, use Last [Head [hash]] as
		     * hash head. */
		    Next [i] = Last [j] ;
		    Last [j] = i ;
		}

		/* if a separate Hhead array is used: *
		Next [i] = Hhead [hash] ;
		Hhead [hash] = i ;
		*/

		Last [i] = hash ;
	    }
	}

	Degree [me] = degme ;

	/* ----------------------------------------------------------------- */
	/* Clear the counter array, W [...], by incrementing wflg. */
	/* ----------------------------------------------------------------- */

	/* make sure that wflg+n does not cause integer overflow */
	lemax =  MAX (lemax, degme) ;
	wflg += lemax ;
	wflg = clear_flag (wflg, wbig, W, n) ;
	/*  at this point, W [0..n-1] < wflg holds */

/* ========================================================================= */
/* SUPERVARIABLE DETECTION */
/* ========================================================================= */

	AMD_DEBUG1 ("Detecting supervariables:\n");
	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n) ;
	    AMD_DEBUG2 ("Consider i "~ID~" nv "~ID~"\n", i, Nv [i]);
	    if (Nv [i] < 0)
	    {
		/* i is a principal variable in Lme */

		/* ---------------------------------------------------------
		 * examine all hash buckets with 2 or more variables.  We do
		 * this by examing all unique hash keys for supervariables in
		 * the pattern Lme of the current element, me
		 * --------------------------------------------------------- */

		/* let i = head of hash bucket, and empty the hash bucket */
		ASSERT (Last [i] >= 0 && Last [i] < n) ;
		hash = Last [i] ;

		/* if Hhead array is not used: */
		j = Head [hash] ;
		if (j == EMPTY)
		{
		    /* hash bucket and degree list are both empty */
		    i = EMPTY ;
		}
		else if (j < EMPTY)
		{
		    /* degree list is empty */
		    i = FLIP (j) ;
		    Head [hash] = EMPTY ;
		}
		else
		{
		    /* degree list is not empty, restore Last [j] of head j */
		    i = Last [j] ;
		    Last [j] = EMPTY ;
		}

		/* if separate Hhead array is used: *
		i = Hhead [hash] ;
		Hhead [hash] = EMPTY ;
		*/

		ASSERT (i >= EMPTY && i < n) ;
		AMD_DEBUG2 ("----i "~ID~" hash "~ID~"\n", i, hash);

		while (i != EMPTY && Next [i] != EMPTY)
		{

		    /* -----------------------------------------------------
		     * this bucket has one or more variables following i.
		     * scan all of them to see if i can absorb any entries
		     * that follow i in hash bucket.  Scatter i into w.
		     * ----------------------------------------------------- */

		    ln = Len [i] ;
		    eln = Elen [i] ;
		    ASSERT (ln >= 0 && eln >= 0) ;
		    ASSERT (Pe [i] >= 0 && Pe [i] < iwlen) ;
		    /* do not flag the first element in the list (me) */
		    for (p = Pe [i] + 1 ; p <= Pe [i] + ln - 1 ; p++)
		    {
			ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
			W [Iw [p]] = wflg ;
		    }

		    /* ----------------------------------------------------- */
		    /* scan every other entry j following i in bucket */
		    /* ----------------------------------------------------- */

		    jlast = i ;
		    j = Next [i] ;
		    ASSERT (j >= EMPTY && j < n) ;

		    while (j != EMPTY)
		    {
			/* ------------------------------------------------- */
			/* check if j and i have identical nonzero pattern */
			/* ------------------------------------------------- */

			AMD_DEBUG3 ("compare i "~ID~" and j "~ID~"\n", i,j);

			/* check if i and j have the same Len and Elen */
			ASSERT (Len [j] >= 0 && Elen [j] >= 0) ;
			ASSERT (Pe [j] >= 0 && Pe [j] < iwlen) ;
			ok = (Len [j] == ln) && (Elen [j] == eln) ;
			/* skip the first element in the list (me) */
			for (p = Pe [j] + 1 ; ok && p <= Pe [j] + ln - 1 ; p++)
			{
			    ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
			    if (W [Iw [p]] != wflg) ok = 0 ;
			}
			if (ok)
			{
			    /* --------------------------------------------- */
			    /* found it!  j can be absorbed into i */
			    /* --------------------------------------------- */

			    AMD_DEBUG1 ("found it! j "~ID~" => i "~ID~"\n", j,i);
			    Pe [j] = FLIP (i) ;
			    /* both Nv [i] and Nv [j] are negated since they */
			    /* are in Lme, and the absolute values of each */
			    /* are the number of variables in i and j: */
			    Nv [i] += Nv [j] ;
			    Nv [j] = 0 ;
			    Elen [j] = EMPTY ;
			    /* delete j from hash bucket */
			    ASSERT (j != Next [j]) ;
			    j = Next [j] ;
			    Next [jlast] = j ;

			}
			else
			{
			    /* j cannot be absorbed into i */
			    jlast = j ;
			    ASSERT (j != Next [j]) ;
			    j = Next [j] ;
			}
			ASSERT (j >= EMPTY && j < n) ;
		    }

		    /* -----------------------------------------------------
		     * no more variables can be absorbed into i
		     * go to next i in bucket and clear flag array
		     * ----------------------------------------------------- */

		    wflg++ ;
		    i = Next [i] ;
		    ASSERT (i >= EMPTY && i < n) ;

		}
	    }
	}
	AMD_DEBUG2 ("detect done\n");

/* ========================================================================= */
/* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVARIABLES FROM ELEMENT */
/* ========================================================================= */

	p = pme1 ;
	nleft = n - nel ;
	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n) ;
	    nvi = -Nv [i] ;
	    AMD_DEBUG3 ("Restore i "~ID~" "~ID~"\n", i, nvi);
	    if (nvi > 0)
	    {
		/* i is a principal variable in Lme */
		/* restore Nv [i] to signify that i is principal */
		Nv [i] = nvi ;

		/* --------------------------------------------------------- */
		/* compute the external degree (add size of current element) */
		/* --------------------------------------------------------- */

		deg = Degree [i] + degme - nvi ;
		deg = MIN (deg, nleft - nvi) ;
		ASSERT (IMPLIES (aggressive, deg > 0) && deg >= 0 && deg < n) ;

		/* --------------------------------------------------------- */
		/* place the supervariable at the head of the degree list */
		/* --------------------------------------------------------- */

		inext = Head [deg] ;
		ASSERT (inext >= EMPTY && inext < n) ;
		if (inext != EMPTY) Last [inext] = i ;
		Next [i] = inext ;
		Last [i] = EMPTY ;
		Head [deg] = i ;

		/* --------------------------------------------------------- */
		/* save the new degree, and find the minimum degree */
		/* --------------------------------------------------------- */

		mindeg = MIN (mindeg, deg) ;
		Degree [i] = deg ;

		/* --------------------------------------------------------- */
		/* place the supervariable in the element pattern */
		/* --------------------------------------------------------- */

		Iw [p++] = i ;

	    }
	}
	AMD_DEBUG2("restore done\n");

/* ========================================================================= */
/* FINALIZE THE NEW ELEMENT */
/* ========================================================================= */

	AMD_DEBUG2 ("ME = "~ID~" DONE\n", me);
	Nv [me] = nvpiv ;
	/* save the length of the list for the new element me */
	Len [me] = p - pme1 ;
	if (Len [me] == 0)
	{
	    /* there is nothing left of the current pivot element */
	    /* it is a root of the assembly tree */
	    Pe [me] = EMPTY ;
	    W [me] = 0 ;
	}
	if (elenme != 0)
	{
	    /* element was not constructed in place: deallocate part of */
	    /* it since newly nonprincipal variables may have been removed */
	    pfree = p ;
	}

	/* The new element has nvpiv pivots and the size of the contribution
	 * block for a multifrontal method is degme-by-degme, not including
	 * the "dense" rows/columns.  If the "dense" rows/columns are included,
	 * the frontal matrix is no larger than
	 * (degme+ndense)-by-(degme+ndense).
	 */

	if (Info != cast(c_float *) NULL)
	{
	    f = nvpiv ;
	    r = degme + ndense ;
	    dmax = MAX (dmax, f + r) ;

	    /* number of nonzeros in L (excluding the diagonal) */
	    lnzme = f*r + (f-1)*f/2 ;
	    lnz += lnzme ;

	    /* number of divide operations for LDL' and for LU */
	    ndiv += lnzme ;

	    /* number of multiply-subtract pairs for LU */
	    s = f*r*r + r*(f-1)*f + (f-1)*f*(2*f-1)/6 ;
	    nms_lu += s ;

	    /* number of multiply-subtract pairs for LDL' */
	    nms_ldl += (s + lnzme)/2 ;
	}

version(NDEBUG){}
else {
	AMD_DEBUG2 ("finalize done nel "~ID~" n "~ID~"\n   ::::\n", nel, n);
	for (pme = Pe [me] ; pme <= Pe [me] + Len [me] - 1 ; pme++)
	{
	      AMD_DEBUG3 (" "~ID~"", Iw [pme]);
	}
	AMD_DEBUG3 ("\n");
} // !NDEBUG

    }

/* ========================================================================= */
/* DONE SELECTING PIVOTS */
/* ========================================================================= */

    if (Info != cast(c_float *) NULL)
    {

	/* count the work to factorize the ndense-by-ndense submatrix */
	f = ndense ;
	dmax = MAX (dmax, cast(c_float) ndense) ;

	/* number of nonzeros in L (excluding the diagonal) */
	lnzme = (f-1)*f/2 ;
	lnz += lnzme ;

	/* number of divide operations for LDL' and for LU */
	ndiv += lnzme ;

	/* number of multiply-subtract pairs for LU */
	s = (f-1)*f*(2*f-1)/6 ;
	nms_lu += s ;

	/* number of multiply-subtract pairs for LDL' */
	nms_ldl += (s + lnzme)/2 ;

	/* number of nz's in L (excl. diagonal) */
	Info [AMD_LNZ] = lnz ;

	/* number of divide ops for LU and LDL' */
	Info [AMD_NDIV] = ndiv ;

	/* number of multiply-subtract pairs for LDL' */
	Info [AMD_NMULTSUBS_LDL] = nms_ldl ;

	/* number of multiply-subtract pairs for LU */
	Info [AMD_NMULTSUBS_LU] = nms_lu ;

	/* number of "dense" rows/columns */
	Info [AMD_NDENSE] = ndense ;

	/* largest front is dmax-by-dmax */
	Info [AMD_DMAX] = dmax ;

	/* number of garbage collections in AMD */
	Info [AMD_NCMPA] = ncmpa ;

	/* successful ordering */
	Info [AMD_STATUS] = AMD_OK ;
    }

/* ========================================================================= */
/* POST-ORDERING */
/* ========================================================================= */

/* -------------------------------------------------------------------------
 * Variables at this point:
 *
 * Pe: holds the elimination tree.  The parent of j is FLIP (Pe [j]),
 *	or EMPTY if j is a root.  The tree holds both elements and
 *	non-principal (unordered) variables absorbed into them.
 *	Dense variables are non-principal and unordered.
 *
 * Elen: holds the size of each element, including the diagonal part.
 *	FLIP (Elen [e]) > 0 if e is an element.  For unordered
 *	variables i, Elen [i] is EMPTY.
 *
 * Nv: Nv [e] > 0 is the number of pivots represented by the element e.
 *	For unordered variables i, Nv [i] is zero.
 *
 * Contents no longer needed:
 *	W, Iw, Len, Degree, Head, Next, Last.
 *
 * The matrix itself has been destroyed.
 *
 * n: the size of the matrix.
 * No other scalars needed (pfree, iwlen, etc.)
 * ------------------------------------------------------------------------- */

    /* restore Pe */
    for (i = 0 ; i < n ; i++)
    {
	Pe [i] = FLIP (Pe [i]) ;
    }

    /* restore Elen, for output information, and for postordering */
    for (i = 0 ; i < n ; i++)
    {
	Elen [i] = FLIP (Elen [i]) ;
    }

/* Now the parent of j is Pe [j], or EMPTY if j is a root.  Elen [e] > 0
 * is the size of element e.  Elen [i] is EMPTY for unordered variable i. */

version(NDEBUG){}
else {
    AMD_DEBUG2 ("\nTree:\n");
    for (i = 0 ; i < n ; i++)
    {
	AMD_DEBUG2 (" "~ID~" parent: "~ID~"   ", i, Pe [i]);
	ASSERT (Pe [i] >= EMPTY && Pe [i] < n) ;
	if (Nv [i] > 0)
	{
	    /* this is an element */
	    e = i ;
	    AMD_DEBUG2 (" element, size is "~ID~"\n", Elen [i]);
	    ASSERT (Elen [e] > 0) ;
	}
	AMD_DEBUG2 ("\n");
    }
    AMD_DEBUG2 ("\nelements:\n");
    for (e = 0 ; e < n ; e++)
    {
	if (Nv [e] > 0)
	{
	    AMD_DEBUG3 ("Element e= "~ID~" size "~ID~" nv "~ID~" \n", e, Elen [e], Nv [e]);
	}
    }
    AMD_DEBUG2 ("\nvariables:\n");
    for (i = 0 ; i < n ; i++)
    {
	Int cnt ;
	if (Nv [i] == 0)
	{
	    AMD_DEBUG3 ("i unordered: "~ID~"\n", i);
	    j = Pe [i] ;
	    cnt = 0 ;
	    AMD_DEBUG3 ("  j: "~ID~"\n", j);
	    if (j == EMPTY)
	    {
		AMD_DEBUG3 ("	i is a dense variable\n");
	    }
	    else
	    {
		ASSERT (j >= 0 && j < n) ;
		while (Nv [j] == 0)
		{
		    AMD_DEBUG3 ("	j : "~ID~"\n", j);
		    j = Pe [j] ;
		    AMD_DEBUG3 ("	j:: "~ID~"\n", j);
		    cnt++ ;
		    if (cnt > n) break ;
		}
		e = j ;
		AMD_DEBUG3 ("	got to e: "~ID~"\n", e);
	    }
	}
    }
} // !NDEBUG

/* ========================================================================= */
/* compress the paths of the variables */
/* ========================================================================= */

    for (i = 0 ; i < n ; i++)
    {
	if (Nv [i] == 0)
	{

	    /* -------------------------------------------------------------
	     * i is an un-ordered row.  Traverse the tree from i until
	     * reaching an element, e.  The element, e, was the principal
	     * supervariable of i and all nodes in the path from i to when e
	     * was selected as pivot.
	     * ------------------------------------------------------------- */

	    AMD_DEBUG1 ("Path compression, i unordered: "~ID~"\n", i);
	    j = Pe [i] ;
	    ASSERT (j >= EMPTY && j < n) ;
	    AMD_DEBUG3 ("	j: "~ID~"\n", j);
	    if (j == EMPTY)
	    {
		/* Skip a dense variable.  It has no parent. */
		AMD_DEBUG3 (("      i is a dense variable\n")) ;
		continue ;
	    }

	    /* while (j is a variable) */
	    while (Nv [j] == 0)
	    {
		AMD_DEBUG3 ("		j : "~ID~"\n", j);
		j = Pe [j] ;
		AMD_DEBUG3 ("		j:: "~ID~"\n", j);
		ASSERT (j >= 0 && j < n) ;
	    }
	    /* got to an element e */
	    e = j ;
	    AMD_DEBUG3 ("got to e: "~ID~"\n", e);

	    /* -------------------------------------------------------------
	     * traverse the path again from i to e, and compress the path
	     * (all nodes point to e).  Path compression allows this code to
	     * compute in O(n) time.
	     * ------------------------------------------------------------- */

	    j = i ;
	    /* while (j is a variable) */
	    while (Nv [j] == 0)
	    {
		jnext = Pe [j] ;
		AMD_DEBUG3 ("j "~ID~" jnext "~ID~"\n", j, jnext);
		Pe [j] = e ;
		j = jnext ;
		ASSERT (j >= 0 && j < n) ;
	    }
	}
    }

/* ========================================================================= */
/* postorder the assembly tree */
/* ========================================================================= */

    AMD_postorder (n, Pe, Nv, Elen,
	W,			/* output order */
	Head, Next, Last) ;	/* workspace */

/* ========================================================================= */
/* compute output permutation and inverse permutation */
/* ========================================================================= */

    /* W [e] = k means that element e is the kth element in the new
     * order.  e is in the range 0 to n-1, and k is in the range 0 to
     * the number of elements.  Use Head for inverse order. */

    for (k = 0 ; k < n ; k++)
    {
	Head [k] = EMPTY ;
	Next [k] = EMPTY ;
    }
    for (e = 0 ; e < n ; e++)
    {
	k = W [e] ;
	ASSERT ((k == EMPTY) == (Nv [e] == 0)) ;
	if (k != EMPTY)
	{
	    ASSERT (k >= 0 && k < n) ;
	    Head [k] = e ;
	}
    }

    /* construct output inverse permutation in Next,
     * and permutation in Last */
    nel = 0 ;
    for (k = 0 ; k < n ; k++)
    {
	e = Head [k] ;
	if (e == EMPTY) break ;
	ASSERT (e >= 0 && e < n && Nv [e] > 0) ;
	Next [e] = nel ;
	nel += Nv [e] ;
    }
    ASSERT (nel == n - ndense) ;

    /* order non-principal variables (dense, & those merged into supervar's) */
    for (i = 0 ; i < n ; i++)
    {
	if (Nv [i] == 0)
	{
	    e = Pe [i] ;
	    ASSERT (e >= EMPTY && e < n) ;
	    if (e != EMPTY)
	    {
		/* This is an unordered variable that was merged
		 * into element e via supernode detection or mass
		 * elimination of i when e became the pivot element.
		 * Place i in order just before e. */
		ASSERT (Next [i] == EMPTY && Nv [e] > 0) ;
		Next [i] = Next [e] ;
		Next [e]++ ;
	    }
	    else
	    {
		/* This is a dense unordered variable, with no parent.
		 * Place it last in the output order. */
		Next [i] = nel++ ;
	    }
	}
    }
    ASSERT (nel == n) ;

    AMD_DEBUG2 (("\n\nPerm:\n")) ;
    for (i = 0 ; i < n ; i++)
    {
	k = Next [i] ;
	ASSERT (k >= 0 && k < n) ;
	Last [k] = i ;
	AMD_DEBUG2 ("   perm ["~ID~"] = "~ID~"\n", k, i);
    }
}

} // version DLONG
else {
    void AMD_2
    (
        int n,		/* A is n-by-n, where n > 0 */
        int *Pe,		/* Pe [0..n-1]: index in Iw of row i on input */
        int *Iw,		/* workspace of size iwlen. Iw [0..pfree-1] holds the matrix on input */
        int *Len,	/* Len [0..n-1]: length for row/column i on input */
        int iwlen,		/* length of Iw. iwlen >= pfree + n */
        int pfree,		/* Iw [pfree ... iwlen-1] is empty on input */
        /* 7 size-n workspaces, not defined on input: */
		int *Nv,		/* the size of each supernode on output */
        int *Next,	/* the output inverse permutation */
        int *Last,	/* the output permutation */
        int *Head,
        int *Elen,/* the size columns of L for each supernode */
        int *Degree,
        int *W,
		/* control parameters and output statistics */
        c_float *Control,	/* array of size AMD_CONTROL */
        c_float *Info	/* array of size AMD_INFO */
    )

	{

    Int deg;
	Int degme;
	Int dext;
	Int lemax;
	Int e;
	Int elenme;
	Int eln;
	Int i;
	Int ilast;
	Int inext;
	Int j;
	Int jlast;
	Int jnext;
	Int k;
	Int knt1;
	Int knt2;
	Int knt3;
	Int lenj;
	Int ln;
	Int me;
	Int mindeg;
	Int nel;
	Int nleft;
	Int nvi;
	Int nvj;
	Int nvpiv;
	Int slenme;
	Int wbig;
	Int we;
	Int wflg;
	Int wnvi;
	Int ok;
	Int ndense;
	Int ncmpa;
	Int dense;
	Int aggressive;

    //unsigned Int hash ;	    /* unsigned, so that hash % n is well defined.*/
	UInt hash ;	    /* unsigned, so that hash % n is well defined.*/

/*
 * deg:		the degree of a variable or element
 * degme:	size, |Lme|, of the current element, me (= Degree [me])
 * dext:	external degree, |Le \ Lme|, of some element e
 * lemax:	largest |Le| seen so far (called dmax in Fortran version)
 * e:		an element
 * elenme:	the length, Elen [me], of element list of pivotal variable
 * eln:		the length, Elen [...], of an element list
 * hash:	the computed value of the hash function
 * i:		a supervariable
 * ilast:	the entry in a link list preceding i
 * inext:	the entry in a link list following i
 * j:		a supervariable
 * jlast:	the entry in a link list preceding j
 * jnext:	the entry in a link list, or path, following j
 * k:		the pivot order of an element or variable
 * knt1:	loop counter used during element construction
 * knt2:	loop counter used during element construction
 * knt3:	loop counter used during compression
 * lenj:	Len [j]
 * ln:		length of a supervariable list
 * me:		current supervariable being eliminated, and the current
 *		    element created by eliminating that supervariable
 * mindeg:	current minimum degree
 * nel:		number of pivots selected so far
 * nleft:	n - nel, the number of nonpivotal rows/columns remaining
 * nvi:		the number of variables in a supervariable i (= Nv [i])
 * nvj:		the number of variables in a supervariable j (= Nv [j])
 * nvpiv:	number of pivots in current element
 * slenme:	number of variables in variable list of pivotal variable
 * wbig:	= (INT_MAX - n) for the int version, (SuiteSparse_long_max - n)
 *                  for the SuiteSparse_long version.  wflg is not allowed to
 *                  be >= wbig.
 * we:		W [e]
 * wflg:	used for flagging the W array.  See description of Iw.
 * wnvi:	wflg - Nv [i]
 * x:		either a supervariable or an element
 *
 * ok:		true if supervariable j can be absorbed into i
 * ndense:	number of "dense" rows/columns
 * dense:	rows/columns with initial degree > dense are considered "dense"
 * aggressive:	true if aggressive absorption is being performed
 * ncmpa:	number of garbage collections

 * ----------------------------------------------------------------------------
 * LOCAL DOUBLES, used for statistical output only (except for alpha):
 * ----------------------------------------------------------------------------
 */

    c_float f, r, ndiv, s, nms_lu, nms_ldl, dmax, alpha, lnz, lnzme ;

/*
 * f:		nvpiv
 * r:		degme + nvpiv
 * ndiv:	number of divisions for LU or LDL' factorizations
 * s:		number of multiply-subtract pairs for LU factorization, for the
 *		    current element me
 * nms_lu	number of multiply-subtract pairs for LU factorization
 * nms_ldl	number of multiply-subtract pairs for LDL' factorization
 * dmax:	the largest number of entries in any column of L, including the
 *		    diagonal
 * alpha:	"dense" degree ratio
 * lnz:		the number of nonzeros in L (excluding the diagonal)
 * lnzme:	the number of nonzeros in L (excl. the diagonal) for the
 *		    current element me

 * ----------------------------------------------------------------------------
 * LOCAL "POINTERS" (indices into the Iw array)
 * ----------------------------------------------------------------------------
*/

    Int p, p1, p2, p3, p4, pdst, pend, pj, pme, pme1, pme2, pn, psrc ;

/*
 * Any parameter (Pe [...] or pfree) or local variable starting with "p" (for
 * Pointer) is an index into Iw, and all indices into Iw use variables starting
 * with "p."  The only exception to this rule is the iwlen input argument.
 *
 * p:           pointer into lots of things
 * p1:          Pe [i] for some variable i (start of element list)
 * p2:          Pe [i] + Elen [i] -  1 for some variable i
 * p3:          index of first supervariable in clean list
 * p4:
 * pdst:        destination pointer, for compression
 * pend:        end of memory to compress
 * pj:          pointer into an element or variable
 * pme:         pointer into the current element (pme1...pme2)
 * pme1:        the current element, me, is stored in Iw [pme1...pme2]
 * pme2:        the end of the current element
 * pn:          pointer into a "clean" variable, also used to compress
 * psrc:        source pointer, for compression
*/

/* ========================================================================= */
/*  INITIALIZATIONS */
/* ========================================================================= */

    /* Note that this restriction on iwlen is slightly more restrictive than
     * what is actually required in AMD_2.  AMD_2 can operate with no elbow
     * room at all, but it will be slow.  For better performance, at least
     * size-n elbow room is enforced. */
    ASSERT (iwlen >= pfree + n) ;
    ASSERT (n > 0) ;

    /* initialize output statistics */
    lnz = 0 ;
    ndiv = 0 ;
    nms_lu = 0 ;
    nms_ldl = 0 ;
    dmax = 1 ;
    me = EMPTY ;

    mindeg = 0 ;
    ncmpa = 0 ;
    nel = 0 ;
    lemax = 0 ;

    /* get control parameters */
    if (Control != cast(c_float *) NULL)
    {
	alpha = Control [AMD_DENSE] ;
	aggressive = (Control [AMD_AGGRESSIVE] != 0) ;
    }
    else
    {
	alpha = AMD_DEFAULT_DENSE ;
	aggressive = AMD_DEFAULT_AGGRESSIVE ;
    }
    /* Note: if alpha is NaN, this is undefined: */
    if (alpha < 0)
    {
	/* only remove completely dense rows/columns */
	dense = n-2 ;
    }
    else
    {
	dense = cast(Int) (alpha * sqrt (cast(c_float) n)) ;
    }
    dense = MAX (16, dense) ;
    dense = MIN (n,  dense) ;
    AMD_DEBUG1 ("\n\nAMD (debug), alpha %g, aggr. "~ID~"\n", alpha, aggressive);

    for (i = 0 ; i < n ; i++)
    {
	Last [i] = EMPTY ;
	Head [i] = EMPTY ;
	Next [i] = EMPTY ;
	/* if separate Hhead array is used for hash buckets: *
	Hhead [i] = EMPTY ;
	*/
	Nv [i] = 1 ;
	W [i] = 1 ;
	Elen [i] = 0 ;
	Degree [i] = Len [i] ;
    }

version(NDEBUG){}
else {
    AMD_DEBUG1 ("\n======Nel "~ID~" initial\n", nel);
    AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, Last,
		Head, Elen, Degree, W, -1) ;
} // !NDEBUG

    /* initialize wflg */
    wbig = Int_MAX - n ;
    wflg = clear_flag (0, wbig, W, n) ;

    /* --------------------------------------------------------------------- */
    /* initialize degree lists and eliminate dense and empty rows */
    /* --------------------------------------------------------------------- */

    ndense = 0 ;

    for (i = 0 ; i < n ; i++)
    {
	deg = Degree [i] ;
	ASSERT (deg >= 0 && deg < n) ;
	if (deg == 0)
	{

	    /* -------------------------------------------------------------
	     * we have a variable that can be eliminated at once because
	     * there is no off-diagonal non-zero in its row.  Note that
	     * Nv [i] = 1 for an empty variable i.  It is treated just
	     * the same as an eliminated element i.
	     * ------------------------------------------------------------- */

	    Elen [i] = FLIP (1) ;
	    nel++ ;
	    Pe [i] = EMPTY ;
	    W [i] = 0 ;

	}
	else if (deg > dense)
	{

	    /* -------------------------------------------------------------
	     * Dense variables are not treated as elements, but as unordered,
	     * non-principal variables that have no parent.  They do not take
	     * part in the postorder, since Nv [i] = 0.  Note that the Fortran
	     * version does not have this option.
	     * ------------------------------------------------------------- */

	    AMD_DEBUG1 ("Dense node "~ID~" degree "~ID~"\n", i, deg);
	    ndense++ ;
	    Nv [i] = 0 ;		/* do not postorder this node */
	    Elen [i] = EMPTY ;
	    nel++ ;
	    Pe [i] = EMPTY ;

	}
	else
	{

	    /* -------------------------------------------------------------
	     * place i in the degree list corresponding to its degree
	     * ------------------------------------------------------------- */

	    inext = Head [deg] ;
	    ASSERT (inext >= EMPTY && inext < n) ;
	    if (inext != EMPTY) Last [inext] = i ;
	    Next [i] = inext ;
	    Head [deg] = i ;

	}
    }

/* ========================================================================= */
/* WHILE (selecting pivots) DO */
/* ========================================================================= */

    while (nel < n)
    {

version(NDEBUG){}
else {
	AMD_DEBUG1 ("\n======Nel "~ID~"\n", nel) ;
	//if (AMD_debug >= 2)
	if (true)		// todo : 
	{
	    AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next,
		    Last, Head, Elen, Degree, W, nel) ;
	}
} // !NDEBUG

/* ========================================================================= */
/* GET PIVOT OF MINIMUM DEGREE */
/* ========================================================================= */

	/* ----------------------------------------------------------------- */
	/* find next supervariable for elimination */
	/* ----------------------------------------------------------------- */

	ASSERT (mindeg >= 0 && mindeg < n) ;
	for (deg = mindeg ; deg < n ; deg++)
	{
	    me = Head [deg] ;
	    if (me != EMPTY) break ;
	}
	mindeg = deg ;
	ASSERT (me >= 0 && me < n) ;
	AMD_DEBUG1 ("=================me: "~ID~"\n", me) ;

	/* ----------------------------------------------------------------- */
	/* remove chosen variable from link list */
	/* ----------------------------------------------------------------- */

	inext = Next [me] ;
	ASSERT (inext >= EMPTY && inext < n) ;
	if (inext != EMPTY) Last [inext] = EMPTY ;
	Head [deg] = inext ;

	/* ----------------------------------------------------------------- */
	/* me represents the elimination of pivots nel to nel+Nv[me]-1. */
	/* place me itself as the first in this set. */
	/* ----------------------------------------------------------------- */

	elenme = Elen [me] ;
	nvpiv = Nv [me] ;
	ASSERT (nvpiv > 0) ;
	nel += nvpiv ;

/* ========================================================================= */
/* CONSTRUCT NEW ELEMENT */
/* ========================================================================= */

	/* -----------------------------------------------------------------
	 * At this point, me is the pivotal supervariable.  It will be
	 * converted into the current element.  Scan list of the pivotal
	 * supervariable, me, setting tree pointers and constructing new list
	 * of supervariables for the new element, me.  p is a pointer to the
	 * current position in the old list.
	 * ----------------------------------------------------------------- */

	/* flag the variable "me" as being in Lme by negating Nv [me] */
	Nv [me] = -nvpiv ;
	degme = 0 ;
	ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;

	if (elenme == 0)
	{

	    /* ------------------------------------------------------------- */
	    /* construct the new element in place */
	    /* ------------------------------------------------------------- */

	    pme1 = Pe [me] ;
	    pme2 = pme1 - 1 ;

	    for (p = pme1 ; p <= pme1 + Len [me] - 1 ; p++)
	    {
		i = Iw [p] ;
		ASSERT (i >= 0 && i < n && Nv [i] >= 0) ;
		nvi = Nv [i] ;
		if (nvi > 0)
		{

		    /* ----------------------------------------------------- */
		    /* i is a principal variable not yet placed in Lme. */
		    /* store i in new list */
		    /* ----------------------------------------------------- */

		    /* flag i as being in Lme by negating Nv [i] */
		    degme += nvi ;
		    Nv [i] = -nvi ;
		    Iw [++pme2] = i ;

		    /* ----------------------------------------------------- */
		    /* remove variable i from degree list. */
		    /* ----------------------------------------------------- */

		    ilast = Last [i] ;
		    inext = Next [i] ;
		    ASSERT (ilast >= EMPTY && ilast < n) ;
		    ASSERT (inext >= EMPTY && inext < n) ;
		    if (inext != EMPTY) Last [inext] = ilast ;
		    if (ilast != EMPTY)
		    {
			Next [ilast] = inext ;
		    }
		    else
		    {
			/* i is at the head of the degree list */
			ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
			Head [Degree [i]] = inext ;
		    }
		}
	    }
	}
	else
	{

	    /* ------------------------------------------------------------- */
	    /* construct the new element in empty space, Iw [pfree ...] */
	    /* ------------------------------------------------------------- */

	    p = Pe [me] ;
	    pme1 = pfree ;
	    slenme = Len [me] - elenme ;

	    for (knt1 = 1 ; knt1 <= elenme + 1 ; knt1++)
	    {

		if (knt1 > elenme)
		{
		    /* search the supervariables in me. */
		    e = me ;
		    pj = p ;
		    ln = slenme ;
		    AMD_DEBUG2 ("Search sv: "~ID~" "~ID~" "~ID~"\n", me,pj,ln);
		}
		else
		{
		    /* search the elements in me. */
		    e = Iw [p++] ;
		    ASSERT (e >= 0 && e < n) ;
		    pj = Pe [e] ;
		    ln = Len [e] ;
		    AMD_DEBUG2 ("Search element e "~ID~" in me "~ID~"\n", e,me);
		    ASSERT (Elen [e] < EMPTY && W [e] > 0 && pj >= 0) ;
		}
		ASSERT (ln >= 0 && (ln == 0 || (pj >= 0 && pj < iwlen))) ;

		/* ---------------------------------------------------------
		 * search for different supervariables and add them to the
		 * new list, compressing when necessary. this loop is
		 * executed once for each element in the list and once for
		 * all the supervariables in the list.
		 * --------------------------------------------------------- */

		for (knt2 = 1 ; knt2 <= ln ; knt2++)
		{
		    i = Iw [pj++] ;
		    ASSERT (i >= 0 && i < n && (i == me || Elen [i] >= EMPTY));
		    nvi = Nv [i] ;
		    AMD_DEBUG2 (": "~ID~" "~ID~" "~ID~" "~ID~"\n", i, Elen [i], Nv [i], wflg);

		    if (nvi > 0)
		    {

			/* ------------------------------------------------- */
			/* compress Iw, if necessary */
			/* ------------------------------------------------- */

			if (pfree >= iwlen)
			{

			    AMD_DEBUG1 ("GARBAGE COLLECTION\n");

			    /* prepare for compressing Iw by adjusting pointers
			     * and lengths so that the lists being searched in
			     * the inner and outer loops contain only the
			     * remaining entries. */

			    Pe [me] = p ;
			    Len [me] -= knt1 ;
			    /* check if nothing left of supervariable me */
			    if (Len [me] == 0) Pe [me] = EMPTY ;
			    Pe [e] = pj ;
			    Len [e] = ln - knt2 ;
			    /* nothing left of element e */
			    if (Len [e] == 0) Pe [e] = EMPTY ;

			    ncmpa++ ;	/* one more garbage collection */

			    /* store first entry of each object in Pe */
			    /* FLIP the first entry in each object */
			    for (j = 0 ; j < n ; j++)
			    {
				pn = Pe [j] ;
				if (pn >= 0)
				{
				    ASSERT (pn >= 0 && pn < iwlen) ;
				    Pe [j] = Iw [pn] ;
				    Iw [pn] = FLIP (j) ;
				}
			    }

			    /* psrc/pdst point to source/destination */
			    psrc = 0 ;
			    pdst = 0 ;
			    pend = pme1 - 1 ;

			    while (psrc <= pend)
			    {
				/* search for next FLIP'd entry */
				j = FLIP (Iw [psrc++]) ;
				if (j >= 0)
				{
				    AMD_DEBUG2 ("Got object j: "~ID~"\n", j);
				    Iw [pdst] = Pe [j] ;
				    Pe [j] = pdst++ ;
				    lenj = Len [j] ;
				    /* copy from source to destination */
				    for (knt3 = 0 ; knt3 <= lenj - 2 ; knt3++)
				    {
					Iw [pdst++] = Iw [psrc++] ;
				    }
				}
			    }

			    /* move the new partially-constructed element */
			    p1 = pdst ;
			    for (psrc = pme1 ; psrc <= pfree-1 ; psrc++)
			    {
				Iw [pdst++] = Iw [psrc] ;
			    }
			    pme1 = p1 ;
			    pfree = pdst ;
			    pj = Pe [e] ;
			    p = Pe [me] ;

			}

			/* ------------------------------------------------- */
			/* i is a principal variable not yet placed in Lme */
			/* store i in new list */
			/* ------------------------------------------------- */

			/* flag i as being in Lme by negating Nv [i] */
			degme += nvi ;
			Nv [i] = -nvi ;
			Iw [pfree++] = i ;
			AMD_DEBUG2 ("     s: "~ID~"     nv "~ID~"\n", i, Nv [i]);

			/* ------------------------------------------------- */
			/* remove variable i from degree link list */
			/* ------------------------------------------------- */

			ilast = Last [i] ;
			inext = Next [i] ;
			ASSERT (ilast >= EMPTY && ilast < n) ;
			ASSERT (inext >= EMPTY && inext < n) ;
			if (inext != EMPTY) Last [inext] = ilast ;
			if (ilast != EMPTY)
			{
			    Next [ilast] = inext ;
			}
			else
			{
			    /* i is at the head of the degree list */
			    ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
			    Head [Degree [i]] = inext ;
			}
		    }
		}

		if (e != me)
		{
		    /* set tree pointer and flag to indicate element e is
		     * absorbed into new element me (the parent of e is me) */
		    AMD_DEBUG1 (" Element "~ID~" => "~ID~"\n", e, me);
		    Pe [e] = FLIP (me) ;
		    W [e] = 0 ;
		}
	    }

	    pme2 = pfree - 1 ;
	}

	/* ----------------------------------------------------------------- */
	/* me has now been converted into an element in Iw [pme1..pme2] */
	/* ----------------------------------------------------------------- */

	/* degme holds the external degree of new element */
	Degree [me] = degme ;
	Pe [me] = pme1 ;
	Len [me] = pme2 - pme1 + 1 ;
	ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;

	Elen [me] = FLIP (nvpiv + degme) ;
	/* FLIP (Elen (me)) is now the degree of pivot (including
	 * diagonal part). */

version(NDEBUG){}
else {
	AMD_DEBUG2 ("New element structure: length= "~ID~"\n", pme2-pme1+1);
	for (pme = pme1 ; pme <= pme2 ; pme++) AMD_DEBUG3 (" "~ID~"", Iw[pme]);
	AMD_DEBUG3 ("\n");
} // !NDEBUG

	/* ----------------------------------------------------------------- */
	/* make sure that wflg is not too large. */
	/* ----------------------------------------------------------------- */

	/* With the current value of wflg, wflg+n must not cause integer
	 * overflow */

	wflg = clear_flag (wflg, wbig, W, n) ;

/* ========================================================================= */
/* COMPUTE (W [e] - wflg) = |Le\Lme| FOR ALL ELEMENTS */
/* ========================================================================= */

	/* -----------------------------------------------------------------
	 * Scan 1:  compute the external degrees of previous elements with
	 * respect to the current element.  That is:
	 *       (W [e] - wflg) = |Le \ Lme|
	 * for each element e that appears in any supervariable in Lme.  The
	 * notation Le refers to the pattern (list of supervariables) of a
	 * previous element e, where e is not yet absorbed, stored in
	 * Iw [Pe [e] + 1 ... Pe [e] + Len [e]].  The notation Lme
	 * refers to the pattern of the current element (stored in
	 * Iw [pme1..pme2]).   If aggressive absorption is enabled, and
	 * (W [e] - wflg) becomes zero, then the element e will be absorbed
	 * in Scan 2.
	 * ----------------------------------------------------------------- */

	AMD_DEBUG2 ("me: ");
	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n) ;
	    eln = Elen [i] ;
	    AMD_DEBUG3 (""~ID~" Elen "~ID~": \n", i, eln);
	    if (eln > 0)
	    {
		/* note that Nv [i] has been negated to denote i in Lme: */
		nvi = -Nv [i] ;
		ASSERT (nvi > 0 && Pe [i] >= 0 && Pe [i] < iwlen) ;
		wnvi = wflg - nvi ;
		for (p = Pe [i] ; p <= Pe [i] + eln - 1 ; p++)
		{
		    e = Iw [p] ;
		    ASSERT (e >= 0 && e < n) ;
		    we = W [e] ;
		    AMD_DEBUG4 ("    e "~ID~" we "~ID~" ", e, we);
		    if (we >= wflg)
		    {
			/* unabsorbed element e has been seen in this loop */
			AMD_DEBUG4 ("    unabsorbed, first time seen");
			we -= nvi ;
		    }
		    else if (we != 0)
		    {
			/* e is an unabsorbed element */
			/* this is the first we have seen e in all of Scan 1 */
			AMD_DEBUG4 ("    unabsorbed");
			we = Degree [e] + wnvi ;
		    }
		    AMD_DEBUG4 ("\n");
		    W [e] = we ;
		}
	    }
	}
	AMD_DEBUG2 ("\n");

/* ========================================================================= */
/* DEGREE UPDATE AND ELEMENT ABSORPTION */
/* ========================================================================= */

	/* -----------------------------------------------------------------
	 * Scan 2:  for each i in Lme, sum up the degree of Lme (which is
	 * degme), plus the sum of the external degrees of each Le for the
	 * elements e appearing within i, plus the supervariables in i.
	 * Place i in hash list.
	 * ----------------------------------------------------------------- */

	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n && Nv [i] < 0 && Elen [i] >= 0) ;
	    AMD_DEBUG2 ("Updating: i "~ID~" "~ID~" "~ID~"\n", i, Elen[i], Len [i]);
	    p1 = Pe [i] ;
	    p2 = p1 + Elen [i] - 1 ;
	    pn = p1 ;
	    hash = 0 ;
	    deg = 0 ;
	    ASSERT (p1 >= 0 && p1 < iwlen && p2 >= -1 && p2 < iwlen) ;

	    /* ------------------------------------------------------------- */
	    /* scan the element list associated with supervariable i */
	    /* ------------------------------------------------------------- */

	    /* UMFPACK/MA38-style approximate degree: */
	    if (aggressive)
	    {
		for (p = p1 ; p <= p2 ; p++)
		{
		    e = Iw [p] ;
		    ASSERT (e >= 0 && e < n) ;
		    we = W [e] ;
		    if (we != 0)
		    {
			/* e is an unabsorbed element */
			/* dext = | Le \ Lme | */
			dext = we - wflg ;
			if (dext > 0)
			{
			    deg += dext ;
			    Iw [pn++] = e ;
			    hash += e ;
			    AMD_DEBUG4 (" e: "~ID~" hash = "~ID~"\n",e,hash);
			}
			else
			{
			    /* external degree of e is zero, absorb e into me*/
			    AMD_DEBUG1 (" Element "~ID~" =>"~ID~" (aggressive)\n", e, me);
			    ASSERT (dext == 0) ;
			    Pe [e] = FLIP (me) ;
			    W [e] = 0 ;
			}
		    }
		}
	    }
	    else
	    {
		for (p = p1 ; p <= p2 ; p++)
		{
		    e = Iw [p] ;
		    ASSERT (e >= 0 && e < n) ;
		    we = W [e] ;
		    if (we != 0)
		    {
			/* e is an unabsorbed element */
			dext = we - wflg ;
			ASSERT (dext >= 0) ;
			deg += dext ;
			Iw [pn++] = e ;
			hash += e ;
			AMD_DEBUG4 ("	e: "~ID~" hash = "~ID~"\n",e,hash);
		    }
		}
	    }

	    /* count the number of elements in i (including me): */
	    Elen [i] = pn - p1 + 1 ;

	    /* ------------------------------------------------------------- */
	    /* scan the supervariables in the list associated with i */
	    /* ------------------------------------------------------------- */

	    /* The bulk of the AMD run time is typically spent in this loop,
	     * particularly if the matrix has many dense rows that are not
	     * removed prior to ordering. */
	    p3 = pn ;
	    p4 = p1 + Len [i] ;
	    for (p = p2 + 1 ; p < p4 ; p++)
	    {
		j = Iw [p] ;
		ASSERT (j >= 0 && j < n) ;
		nvj = Nv [j] ;
		if (nvj > 0)
		{
		    /* j is unabsorbed, and not in Lme. */
		    /* add to degree and add to new list */
		    deg += nvj ;
		    Iw [pn++] = j ;
		    hash += j ;
		    AMD_DEBUG4 ("  s: "~ID~" hash "~ID~" Nv[j]= "~ID~"\n", j, hash, nvj) ;
		}
	    }

	    /* ------------------------------------------------------------- */
	    /* update the degree and check for mass elimination */
	    /* ------------------------------------------------------------- */

	    /* with aggressive absorption, deg==0 is identical to the
	     * Elen [i] == 1 && p3 == pn test, below. */
	    ASSERT (IMPLIES (aggressive, (deg==0) == (Elen[i]==1 && p3==pn))) ;

	    if (Elen [i] == 1 && p3 == pn)
	    {

		/* --------------------------------------------------------- */
		/* mass elimination */
		/* --------------------------------------------------------- */

		/* There is nothing left of this node except for an edge to
		 * the current pivot element.  Elen [i] is 1, and there are
		 * no variables adjacent to node i.  Absorb i into the
		 * current pivot element, me.  Note that if there are two or
		 * more mass eliminations, fillin due to mass elimination is
		 * possible within the nvpiv-by-nvpiv pivot block.  It is this
		 * step that causes AMD's analysis to be an upper bound.
		 *
		 * The reason is that the selected pivot has a lower
		 * approximate degree than the true degree of the two mass
		 * eliminated nodes.  There is no edge between the two mass
		 * eliminated nodes.  They are merged with the current pivot
		 * anyway.
		 *
		 * No fillin occurs in the Schur complement, in any case,
		 * and this effect does not decrease the quality of the
		 * ordering itself, just the quality of the nonzero and
		 * flop count analysis.  It also means that the post-ordering
		 * is not an exact elimination tree post-ordering. */

		AMD_DEBUG1 ("  MASS i "~ID~" => parent e "~ID~"\n", i, me);
		Pe [i] = FLIP (me) ;
		nvi = -Nv [i] ;
		degme -= nvi ;
		nvpiv += nvi ;
		nel += nvi ;
		Nv [i] = 0 ;
		Elen [i] = EMPTY ;

	    }
	    else
	    {

		/* --------------------------------------------------------- */
		/* update the upper-bound degree of i */
		/* --------------------------------------------------------- */

		/* the following degree does not yet include the size
		 * of the current element, which is added later: */

		Degree [i] = MIN (Degree [i], deg) ;

		/* --------------------------------------------------------- */
		/* add me to the list for i */
		/* --------------------------------------------------------- */

		/* move first supervariable to end of list */
		Iw [pn] = Iw [p3] ;
		/* move first element to end of element part of list */
		Iw [p3] = Iw [p1] ;
		/* add new element, me, to front of list. */
		Iw [p1] = me ;
		/* store the new length of the list in Len [i] */
		Len [i] = pn - p1 + 1 ;

		/* --------------------------------------------------------- */
		/* place in hash bucket.  Save hash key of i in Last [i]. */
		/* --------------------------------------------------------- */

		/* NOTE: this can fail if hash is negative, because the ANSI C
		 * standard does not define a % b when a and/or b are negative.
		 * That's why hash is defined as an unsigned Int, to avoid this
		 * problem. */
		hash = hash % n ;
		ASSERT ((cast(Int) hash) >= 0 && (cast(Int) hash) < n) ;

		/* if the Hhead array is not used: */
		j = Head [hash] ;
		if (j <= EMPTY)
		{
		    /* degree list is empty, hash head is FLIP (j) */
		    Next [i] = FLIP (j) ;
		    Head [hash] = FLIP (i) ;
		}
		else
		{
		    /* degree list is not empty, use Last [Head [hash]] as
		     * hash head. */
		    Next [i] = Last [j] ;
		    Last [j] = i ;
		}

		/* if a separate Hhead array is used: *
		Next [i] = Hhead [hash] ;
		Hhead [hash] = i ;
		*/

		Last [i] = hash ;
	    }
	}

	Degree [me] = degme ;

	/* ----------------------------------------------------------------- */
	/* Clear the counter array, W [...], by incrementing wflg. */
	/* ----------------------------------------------------------------- */

	/* make sure that wflg+n does not cause integer overflow */
	lemax =  MAX (lemax, degme) ;
	wflg += lemax ;
	wflg = clear_flag (wflg, wbig, W, n) ;
	/*  at this point, W [0..n-1] < wflg holds */

/* ========================================================================= */
/* SUPERVARIABLE DETECTION */
/* ========================================================================= */

	AMD_DEBUG1 ("Detecting supervariables:\n");
	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n) ;
	    AMD_DEBUG2 ("Consider i "~ID~" nv "~ID~"\n", i, Nv [i]);
	    if (Nv [i] < 0)
	    {
		/* i is a principal variable in Lme */

		/* ---------------------------------------------------------
		 * examine all hash buckets with 2 or more variables.  We do
		 * this by examing all unique hash keys for supervariables in
		 * the pattern Lme of the current element, me
		 * --------------------------------------------------------- */

		/* let i = head of hash bucket, and empty the hash bucket */
		ASSERT (Last [i] >= 0 && Last [i] < n) ;
		hash = Last [i] ;

		/* if Hhead array is not used: */
		j = Head [hash] ;
		if (j == EMPTY)
		{
		    /* hash bucket and degree list are both empty */
		    i = EMPTY ;
		}
		else if (j < EMPTY)
		{
		    /* degree list is empty */
		    i = FLIP (j) ;
		    Head [hash] = EMPTY ;
		}
		else
		{
		    /* degree list is not empty, restore Last [j] of head j */
		    i = Last [j] ;
		    Last [j] = EMPTY ;
		}

		/* if separate Hhead array is used: *
		i = Hhead [hash] ;
		Hhead [hash] = EMPTY ;
		*/

		ASSERT (i >= EMPTY && i < n) ;
		AMD_DEBUG2 ("----i "~ID~" hash "~ID~"\n", i, hash);

		while (i != EMPTY && Next [i] != EMPTY)
		{

		    /* -----------------------------------------------------
		     * this bucket has one or more variables following i.
		     * scan all of them to see if i can absorb any entries
		     * that follow i in hash bucket.  Scatter i into w.
		     * ----------------------------------------------------- */

		    ln = Len [i] ;
		    eln = Elen [i] ;
		    ASSERT (ln >= 0 && eln >= 0) ;
		    ASSERT (Pe [i] >= 0 && Pe [i] < iwlen) ;
		    /* do not flag the first element in the list (me) */
		    for (p = Pe [i] + 1 ; p <= Pe [i] + ln - 1 ; p++)
		    {
			ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
			W [Iw [p]] = wflg ;
		    }

		    /* ----------------------------------------------------- */
		    /* scan every other entry j following i in bucket */
		    /* ----------------------------------------------------- */

		    jlast = i ;
		    j = Next [i] ;
		    ASSERT (j >= EMPTY && j < n) ;

		    while (j != EMPTY)
		    {
			/* ------------------------------------------------- */
			/* check if j and i have identical nonzero pattern */
			/* ------------------------------------------------- */

			AMD_DEBUG3 ("compare i "~ID~" and j "~ID~"\n", i,j);

			/* check if i and j have the same Len and Elen */
			ASSERT (Len [j] >= 0 && Elen [j] >= 0) ;
			ASSERT (Pe [j] >= 0 && Pe [j] < iwlen) ;
			ok = (Len [j] == ln) && (Elen [j] == eln) ;
			/* skip the first element in the list (me) */
			for (p = Pe [j] + 1 ; ok && p <= Pe [j] + ln - 1 ; p++)
			{
			    ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
			    if (W [Iw [p]] != wflg) ok = 0 ;
			}
			if (ok)
			{
			    /* --------------------------------------------- */
			    /* found it!  j can be absorbed into i */
			    /* --------------------------------------------- */

			    AMD_DEBUG1 ("found it! j "~ID~" => i "~ID~"\n", j,i);
			    Pe [j] = FLIP (i) ;
			    /* both Nv [i] and Nv [j] are negated since they */
			    /* are in Lme, and the absolute values of each */
			    /* are the number of variables in i and j: */
			    Nv [i] += Nv [j] ;
			    Nv [j] = 0 ;
			    Elen [j] = EMPTY ;
			    /* delete j from hash bucket */
			    ASSERT (j != Next [j]) ;
			    j = Next [j] ;
			    Next [jlast] = j ;

			}
			else
			{
			    /* j cannot be absorbed into i */
			    jlast = j ;
			    ASSERT (j != Next [j]) ;
			    j = Next [j] ;
			}
			ASSERT (j >= EMPTY && j < n) ;
		    }

		    /* -----------------------------------------------------
		     * no more variables can be absorbed into i
		     * go to next i in bucket and clear flag array
		     * ----------------------------------------------------- */

		    wflg++ ;
		    i = Next [i] ;
		    ASSERT (i >= EMPTY && i < n) ;

		}
	    }
	}
	AMD_DEBUG2 ("detect done\n");

/* ========================================================================= */
/* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVARIABLES FROM ELEMENT */
/* ========================================================================= */

	p = pme1 ;
	nleft = n - nel ;
	for (pme = pme1 ; pme <= pme2 ; pme++)
	{
	    i = Iw [pme] ;
	    ASSERT (i >= 0 && i < n) ;
	    nvi = -Nv [i] ;
	    AMD_DEBUG3 ("Restore i "~ID~" "~ID~"\n", i, nvi);
	    if (nvi > 0)
	    {
		/* i is a principal variable in Lme */
		/* restore Nv [i] to signify that i is principal */
		Nv [i] = nvi ;

		/* --------------------------------------------------------- */
		/* compute the external degree (add size of current element) */
		/* --------------------------------------------------------- */

		deg = Degree [i] + degme - nvi ;
		deg = MIN (deg, nleft - nvi) ;
		ASSERT (IMPLIES (aggressive, deg > 0) && deg >= 0 && deg < n) ;

		/* --------------------------------------------------------- */
		/* place the supervariable at the head of the degree list */
		/* --------------------------------------------------------- */

		inext = Head [deg] ;
		ASSERT (inext >= EMPTY && inext < n) ;
		if (inext != EMPTY) Last [inext] = i ;
		Next [i] = inext ;
		Last [i] = EMPTY ;
		Head [deg] = i ;

		/* --------------------------------------------------------- */
		/* save the new degree, and find the minimum degree */
		/* --------------------------------------------------------- */

		mindeg = MIN (mindeg, deg) ;
		Degree [i] = deg ;

		/* --------------------------------------------------------- */
		/* place the supervariable in the element pattern */
		/* --------------------------------------------------------- */

		Iw [p++] = i ;

	    }
	}
	AMD_DEBUG2("restore done\n");

/* ========================================================================= */
/* FINALIZE THE NEW ELEMENT */
/* ========================================================================= */

	AMD_DEBUG2 ("ME = "~ID~" DONE\n", me);
	Nv [me] = nvpiv ;
	/* save the length of the list for the new element me */
	Len [me] = p - pme1 ;
	if (Len [me] == 0)
	{
	    /* there is nothing left of the current pivot element */
	    /* it is a root of the assembly tree */
	    Pe [me] = EMPTY ;
	    W [me] = 0 ;
	}
	if (elenme != 0)
	{
	    /* element was not constructed in place: deallocate part of */
	    /* it since newly nonprincipal variables may have been removed */
	    pfree = p ;
	}

	/* The new element has nvpiv pivots and the size of the contribution
	 * block for a multifrontal method is degme-by-degme, not including
	 * the "dense" rows/columns.  If the "dense" rows/columns are included,
	 * the frontal matrix is no larger than
	 * (degme+ndense)-by-(degme+ndense).
	 */

	if (Info != cast(c_float *) NULL)
	{
	    f = nvpiv ;
	    r = degme + ndense ;
	    dmax = MAX (dmax, f + r) ;

	    /* number of nonzeros in L (excluding the diagonal) */
	    lnzme = f*r + (f-1)*f/2 ;
	    lnz += lnzme ;

	    /* number of divide operations for LDL' and for LU */
	    ndiv += lnzme ;

	    /* number of multiply-subtract pairs for LU */
	    s = f*r*r + r*(f-1)*f + (f-1)*f*(2*f-1)/6 ;
	    nms_lu += s ;

	    /* number of multiply-subtract pairs for LDL' */
	    nms_ldl += (s + lnzme)/2 ;
	}

version(NDEBUG){}
else {
	AMD_DEBUG2 ("finalize done nel "~ID~" n "~ID~"\n   ::::\n", nel, n);
	for (pme = Pe [me] ; pme <= Pe [me] + Len [me] - 1 ; pme++)
	{
	      AMD_DEBUG3 (" "~ID~"", Iw [pme]);
	}
	AMD_DEBUG3 ("\n");
} // !NDEBUG

    }

/* ========================================================================= */
/* DONE SELECTING PIVOTS */
/* ========================================================================= */

    if (Info != cast(c_float *) NULL)
    {

	/* count the work to factorize the ndense-by-ndense submatrix */
	f = ndense ;
	dmax = MAX (dmax, cast(c_float) ndense) ;

	/* number of nonzeros in L (excluding the diagonal) */
	lnzme = (f-1)*f/2 ;
	lnz += lnzme ;

	/* number of divide operations for LDL' and for LU */
	ndiv += lnzme ;

	/* number of multiply-subtract pairs for LU */
	s = (f-1)*f*(2*f-1)/6 ;
	nms_lu += s ;

	/* number of multiply-subtract pairs for LDL' */
	nms_ldl += (s + lnzme)/2 ;

	/* number of nz's in L (excl. diagonal) */
	Info [AMD_LNZ] = lnz ;

	/* number of divide ops for LU and LDL' */
	Info [AMD_NDIV] = ndiv ;

	/* number of multiply-subtract pairs for LDL' */
	Info [AMD_NMULTSUBS_LDL] = nms_ldl ;

	/* number of multiply-subtract pairs for LU */
	Info [AMD_NMULTSUBS_LU] = nms_lu ;

	/* number of "dense" rows/columns */
	Info [AMD_NDENSE] = ndense ;

	/* largest front is dmax-by-dmax */
	Info [AMD_DMAX] = dmax ;

	/* number of garbage collections in AMD */
	Info [AMD_NCMPA] = ncmpa ;

	/* successful ordering */
	Info [AMD_STATUS] = AMD_OK ;
    }

/* ========================================================================= */
/* POST-ORDERING */
/* ========================================================================= */

/* -------------------------------------------------------------------------
 * Variables at this point:
 *
 * Pe: holds the elimination tree.  The parent of j is FLIP (Pe [j]),
 *	or EMPTY if j is a root.  The tree holds both elements and
 *	non-principal (unordered) variables absorbed into them.
 *	Dense variables are non-principal and unordered.
 *
 * Elen: holds the size of each element, including the diagonal part.
 *	FLIP (Elen [e]) > 0 if e is an element.  For unordered
 *	variables i, Elen [i] is EMPTY.
 *
 * Nv: Nv [e] > 0 is the number of pivots represented by the element e.
 *	For unordered variables i, Nv [i] is zero.
 *
 * Contents no longer needed:
 *	W, Iw, Len, Degree, Head, Next, Last.
 *
 * The matrix itself has been destroyed.
 *
 * n: the size of the matrix.
 * No other scalars needed (pfree, iwlen, etc.)
 * ------------------------------------------------------------------------- */

    /* restore Pe */
    for (i = 0 ; i < n ; i++)
    {
	Pe [i] = FLIP (Pe [i]) ;
    }

    /* restore Elen, for output information, and for postordering */
    for (i = 0 ; i < n ; i++)
    {
	Elen [i] = FLIP (Elen [i]) ;
    }

/* Now the parent of j is Pe [j], or EMPTY if j is a root.  Elen [e] > 0
 * is the size of element e.  Elen [i] is EMPTY for unordered variable i. */

version(NDEBUG){}
else {
    AMD_DEBUG2 ("\nTree:\n");
    for (i = 0 ; i < n ; i++)
    {
	AMD_DEBUG2 (" "~ID~" parent: "~ID~"   ", i, Pe [i]);
	ASSERT (Pe [i] >= EMPTY && Pe [i] < n) ;
	if (Nv [i] > 0)
	{
	    /* this is an element */
	    e = i ;
	    AMD_DEBUG2 (" element, size is "~ID~"\n", Elen [i]);
	    ASSERT (Elen [e] > 0) ;
	}
	AMD_DEBUG2 ("\n");
    }
    AMD_DEBUG2 ("\nelements:\n");
    for (e = 0 ; e < n ; e++)
    {
	if (Nv [e] > 0)
	{
	    AMD_DEBUG3 ("Element e= "~ID~" size "~ID~" nv "~ID~" \n", e, Elen [e], Nv [e]);
	}
    }
    AMD_DEBUG2 ("\nvariables:\n");
    for (i = 0 ; i < n ; i++)
    {
	Int cnt ;
	if (Nv [i] == 0)
	{
	    AMD_DEBUG3 ("i unordered: "~ID~"\n", i);
	    j = Pe [i] ;
	    cnt = 0 ;
	    AMD_DEBUG3 ("  j: "~ID~"\n", j);
	    if (j == EMPTY)
	    {
		AMD_DEBUG3 ("	i is a dense variable\n");
	    }
	    else
	    {
		ASSERT (j >= 0 && j < n) ;
		while (Nv [j] == 0)
		{
		    AMD_DEBUG3 ("	j : "~ID~"\n", j);
		    j = Pe [j] ;
		    AMD_DEBUG3 ("	j:: "~ID~"\n", j);
		    cnt++ ;
		    if (cnt > n) break ;
		}
		e = j ;
		AMD_DEBUG3 ("	got to e: "~ID~"\n", e);
	    }
	}
    }
} // !NDEBUG

/* ========================================================================= */
/* compress the paths of the variables */
/* ========================================================================= */

    for (i = 0 ; i < n ; i++)
    {
	if (Nv [i] == 0)
	{

	    /* -------------------------------------------------------------
	     * i is an un-ordered row.  Traverse the tree from i until
	     * reaching an element, e.  The element, e, was the principal
	     * supervariable of i and all nodes in the path from i to when e
	     * was selected as pivot.
	     * ------------------------------------------------------------- */

	    AMD_DEBUG1 ("Path compression, i unordered: "~ID~"\n", i);
	    j = Pe [i] ;
	    ASSERT (j >= EMPTY && j < n) ;
	    AMD_DEBUG3 ("	j: "~ID~"\n", j);
	    if (j == EMPTY)
	    {
		/* Skip a dense variable.  It has no parent. */
		AMD_DEBUG3 (("      i is a dense variable\n")) ;
		continue ;
	    }

	    /* while (j is a variable) */
	    while (Nv [j] == 0)
	    {
		AMD_DEBUG3 ("		j : "~ID~"\n", j);
		j = Pe [j] ;
		AMD_DEBUG3 ("		j:: "~ID~"\n", j);
		ASSERT (j >= 0 && j < n) ;
	    }
	    /* got to an element e */
	    e = j ;
	    AMD_DEBUG3 ("got to e: "~ID~"\n", e);

	    /* -------------------------------------------------------------
	     * traverse the path again from i to e, and compress the path
	     * (all nodes point to e).  Path compression allows this code to
	     * compute in O(n) time.
	     * ------------------------------------------------------------- */

	    j = i ;
	    /* while (j is a variable) */
	    while (Nv [j] == 0)
	    {
		jnext = Pe [j] ;
		AMD_DEBUG3 ("j "~ID~" jnext "~ID~"\n", j, jnext);
		Pe [j] = e ;
		j = jnext ;
		ASSERT (j >= 0 && j < n) ;
	    }
	}
    }

/* ========================================================================= */
/* postorder the assembly tree */
/* ========================================================================= */

    AMD_postorder (n, Pe, Nv, Elen,
	W,			/* output order */
	Head, Next, Last) ;	/* workspace */

/* ========================================================================= */
/* compute output permutation and inverse permutation */
/* ========================================================================= */

    /* W [e] = k means that element e is the kth element in the new
     * order.  e is in the range 0 to n-1, and k is in the range 0 to
     * the number of elements.  Use Head for inverse order. */

    for (k = 0 ; k < n ; k++)
    {
	Head [k] = EMPTY ;
	Next [k] = EMPTY ;
    }
    for (e = 0 ; e < n ; e++)
    {
	k = W [e] ;
	ASSERT ((k == EMPTY) == (Nv [e] == 0)) ;
	if (k != EMPTY)
	{
	    ASSERT (k >= 0 && k < n) ;
	    Head [k] = e ;
	}
    }

    /* construct output inverse permutation in Next,
     * and permutation in Last */
    nel = 0 ;
    for (k = 0 ; k < n ; k++)
    {
	e = Head [k] ;
	if (e == EMPTY) break ;
	ASSERT (e >= 0 && e < n && Nv [e] > 0) ;
	Next [e] = nel ;
	nel += Nv [e] ;
    }
    ASSERT (nel == n - ndense) ;

    /* order non-principal variables (dense, & those merged into supervar's) */
    for (i = 0 ; i < n ; i++)
    {
	if (Nv [i] == 0)
	{
	    e = Pe [i] ;
	    ASSERT (e >= EMPTY && e < n) ;
	    if (e != EMPTY)
	    {
		/* This is an unordered variable that was merged
		 * into element e via supernode detection or mass
		 * elimination of i when e became the pivot element.
		 * Place i in order just before e. */
		ASSERT (Next [i] == EMPTY && Nv [e] > 0) ;
		Next [i] = Next [e] ;
		Next [e]++ ;
	    }
	    else
	    {
		/* This is a dense unordered variable, with no parent.
		 * Place it last in the output order. */
		Next [i] = nel++ ;
	    }
	}
    }
    ASSERT (nel == n) ;

    AMD_DEBUG2 (("\n\nPerm:\n")) ;
    for (i = 0 ; i < n ; i++)
    {
	k = Next [i] ;
	ASSERT (k >= 0 && k < n) ;
	Last [k] = i ;
	AMD_DEBUG2 ("   perm ["~ID~"] = "~ID~"\n", k, i);
    }
}

}	// version DLONG