--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umf_kernel_init.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,1056 @@
+/* ========================================================================== */
+/* === UMF_kernel_init ====================================================== */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.  CISE Dept,   */
+/* Univ. of Florida.  All Rights Reserved.  See ../Doc/License for License.   */
+/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
+/* -------------------------------------------------------------------------- */
+
+/*
+    Initialize the kernel: scale the matrix, load the initial elements, and
+    build the tuple lists.
+
+    Returns TRUE if successful, FALSE if out of memory or if the pattern has
+    changed since UMFPACK_*symbolic.  UMFPACK_numeric allocates at least enough
+    space for UMF_kernel_init to succeed; otherwise it does not call
+    UMF_kernel_init.  So an out-of-memory condition means that the pattern must
+    have gotten larger.
+*/
+
+#include "umf_internal.h"
+#include "umf_tuple_lengths.h"
+#include "umf_build_tuples.h"
+#include "umf_mem_init_memoryspace.h"
+#include "umf_mem_alloc_element.h"
+#include "umf_mem_alloc_head_block.h"
+#include "umf_mem_alloc_tail_block.h"
+#include "umf_mem_free_tail_block.h"
+#include "umf_scale.h"
+
+/* ========================================================================== */
+/* === packsp =============================================================== */
+/* ========================================================================== */
+
+/* remove zero or small entries from a column of L or a row of U */
+
+PRIVATE Int packsp	/* returns new value of pnew */
+(
+    Int pnew,		/* index into Memory of next free space */
+    Int *p_p,		/* ptr to index of old pattern in Memory on input,
+			   new pattern on output */
+    Int *p_len,		/* ptr to length of old pattern on input,
+			   new pattern on output */
+    Int drop,		/* TRUE if small nonzero entries are to be dropped */
+    double droptol,	/* the drop tolerance */
+    Unit *Memory	/* contains the sparse vector on input and output */
+)
+{
+    Entry x ;
+    double s ;
+    Entry *Bx, *Bx2 ;
+    Int p, i, len, len_new, *Bi, *Bi2 ;
+
+    /* get the pointers to the sparse vector, and its length */
+    p = *p_p ;
+    len = *p_len ;
+    Bi = (Int   *) (Memory + p) ; p += UNITS (Int,   len) ;
+    Bx = (Entry *) (Memory + p) ; p += UNITS (Entry, len) ;
+    DEBUGm4 (("  p "ID" len "ID" pnew "ID"\n", p, len, pnew)) ;
+
+    /* the vector resides in Bi [0..len-1] and Bx [0..len-1] */
+
+    /* first, compact the vector in place */
+    len_new = 0 ;
+    for (p = 0 ; p < len ; p++)
+    {
+	i = Bi [p] ;
+	x = Bx [p] ;
+	DEBUGm4 (("    old vector: i "ID" value: ", i)) ;
+	EDEBUGk (-4, x) ;
+	DEBUGm4 (("\n")) ;
+	ASSERT (i >= 0) ;
+	/* skip if zero or below drop tolerance */
+	if (IS_ZERO (x)) continue ;
+	if (drop)
+	{
+	    APPROX_ABS (s, x) ;
+	    if (s <= droptol) continue ;
+	}
+	/* store the value back into the vector */
+	if (len_new != p)
+	{
+	    Bi [len_new] = i ;
+	    Bx [len_new] = x ;
+	}
+	len_new++ ;
+    }
+    ASSERT (len_new <= len) ;
+
+    /* the vector is now in Bi [0..len_new-1] and Bx [0..len_new-1] */
+
+#ifndef NDEBUG
+    for (p = 0 ; p < len_new ; p++)
+    {
+	DEBUGm4 (("    new vector: i "ID" value: ", Bi [p])) ;
+	EDEBUGk (-4, Bx [p]) ;
+	DEBUGm4 (("\n")) ;
+	ASSERT (Bi [p] >= 0) ;
+    }
+#endif
+
+    /* allocate new space for the compacted vector */
+    *p_p = pnew ;
+    *p_len = len_new ;
+    Bi2 = (Int   *) (Memory + pnew) ; pnew += UNITS (Int,   len_new) ;
+    Bx2 = (Entry *) (Memory + pnew) ; pnew += UNITS (Entry, len_new) ;
+    DEBUGm4 (("  pnew "ID" len_new "ID"\n", pnew, len_new)) ;
+
+    /* shift the vector upwards, into its new space */
+    for (p = 0 ; p < len_new ; p++)
+    {
+	Bi2 [p] = Bi [p] ;
+    }
+    for (p = 0 ; p < len_new ; p++)
+    {
+	Bx2 [p] = Bx [p] ;
+    }
+
+#ifndef NDEBUG
+    for (p = 0 ; p < len_new ; p++)
+    {
+	DEBUGm4 (("    packed vec: i "ID" value: ", Bi2 [p])) ;
+	EDEBUGk (-4, Bx2 [p]) ;
+	DEBUGm4 (("\n")) ;
+	ASSERT (Bi2 [p] >= 0) ;
+    }
+#endif
+
+    /* return the pointer to the space just after the new vector */
+    return (pnew) ;
+}
+
+
+/* ========================================================================== */
+/* === UMF_kernel_init ====================================================== */
+/* ========================================================================== */
+
+GLOBAL Int UMF_kernel_init
+(
+    const Int Ap [ ],		/* user's input matrix (not modified) */
+    const Int Ai [ ],
+    const double Ax [ ],
+#ifdef COMPLEX
+    const double Az [ ],
+#endif
+    NumericType *Numeric,
+    WorkType *Work,
+    SymbolicType *Symbolic
+)
+{
+    /* ---------------------------------------------------------------------- */
+    /* local variables */
+    /* ---------------------------------------------------------------------- */
+
+    Entry x, pivot_value ;
+    double unused = 0, rsmin, rsmax, rs, droptol ;
+    Entry *D, *C, *Lval, **Rpx ;
+    double *Rs ;
+    Int row, k, oldcol, size, e, p1, p2, p, nz, *Rows, *Cols, *E, i, *Upos,
+	*Lpos, n_row, n_col, *Wp, *Cperm_init, *Frpos, *Fcpos, *Row_degree, nn,
+	*Row_tlen, *Col_degree, *Col_tlen, oldrow, newrow, ilast, *Wrp,
+	*Rperm_init, col, n_inner, prefer_diagonal, *Diagonal_map, nempty,
+	*Diagonal_imap, fixQ, rdeg, cdeg, nempty_col, *Esize, esize, pnew,
+	*Lip, *Uip, *Lilen, *Uilen, llen, pa, *Cdeg, *Rdeg, n1, clen, do_scale,
+	lnz, unz, lip, uip, k1, *Rperm, *Cperm, pivcol, *Li, lilen, drop,
+	**Rpi, nempty_row, dense_row_threshold, empty_elements, rpi, rpx ;
+    Element *ep ;
+    Unit *Memory ;
+#ifdef COMPLEX
+    Int split = SPLIT (Az) ;
+#endif
+#ifndef NRECIPROCAL
+    Int do_recip = FALSE ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* get parameters */
+    /* ---------------------------------------------------------------------- */
+
+    DEBUG0 (("KERNEL INIT\n")) ;
+
+    n_row = Symbolic->n_row ;
+    n_col = Symbolic->n_col ;
+    nn = MAX (n_row, n_col) ;
+    n_inner = MIN (n_row, n_col) ;
+    nempty_col = Symbolic->nempty_col ;
+    nempty_row = Symbolic->nempty_row ;
+    nempty = MIN (nempty_row, nempty_col) ;
+    ASSERT (n_row > 0 && n_col > 0) ;
+    Cperm_init = Symbolic->Cperm_init ;
+    Rperm_init = Symbolic->Rperm_init ;
+    Cdeg = Symbolic->Cdeg ;
+    Rdeg = Symbolic->Rdeg ;
+    n1 = Symbolic->n1 ;
+    dense_row_threshold = Symbolic->dense_row_threshold ;
+    DEBUG0 (("Singletons: "ID"\n", n1)) ;
+    Work->nforced = 0 ;
+    Work->ndiscard = 0 ;
+    Work->noff_diagonal = 0 ;
+
+    nz = Ap [n_col] ;
+    if (nz < 0 || Ap [0] != 0 || nz != Symbolic->nz)
+    {
+	DEBUGm4 (("nz or Ap [0] bad\n")) ;
+	return (FALSE) ;	/* pattern changed */
+    }
+
+    prefer_diagonal = Symbolic->prefer_diagonal ;
+    Diagonal_map = Work->Diagonal_map ;
+    Diagonal_imap = Work->Diagonal_imap ;
+
+    /* ---------------------------------------------------------------------- */
+    /* initialize the Numeric->Memory space for LU, elements, and tuples */
+    /* ---------------------------------------------------------------------- */
+
+    UMF_mem_init_memoryspace (Numeric) ;
+    DEBUG1 (("Kernel init head usage, before allocs: "ID"\n", Numeric->ihead)) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* initialize the Work and Numeric objects */
+    /* ---------------------------------------------------------------------- */
+
+    /* current front is empty */
+    Work->fnpiv = 0 ;
+    Work->fncols = 0 ;
+    Work->fnrows = 0 ;
+    Work->fncols_max = 0 ;
+    Work->fnrows_max = 0 ;
+    Work->fnzeros = 0 ;
+    Work->fcurr_size = 0 ;
+    Work->fnr_curr = 0 ;
+    Work->fnc_curr = 0 ;
+
+    Work->nz = nz ;
+    Work->prior_element = EMPTY ;
+    Work->ulen = 0 ;
+    Work->llen = 0 ;
+    Work->npiv = n1 ;
+    Work->frontid = 0 ;
+    Work->nextcand = n1 ;
+
+    Memory = Numeric->Memory ;
+    Rperm = Numeric->Rperm ;
+    Cperm = Numeric->Cperm ;
+    Row_degree = Numeric->Rperm ;
+    Col_degree = Numeric->Cperm ;
+    /* Row_tuples = Numeric->Uip ; not needed */
+    Row_tlen   = Numeric->Uilen ;
+    /* Col_tuples = Numeric->Lip ; not needed */
+    Col_tlen   = Numeric->Lilen ;
+
+    Lip = Numeric->Lip ;
+    Uip = Numeric->Uip ;
+    Lilen = Numeric->Lilen ;
+    Uilen = Numeric->Uilen ;
+
+    Frpos = Work->Frpos ;
+    Fcpos = Work->Fcpos ;
+    Wp = Work->Wp ;
+    Wrp = Work->Wrp ;
+
+    D = Numeric->D ;
+    Upos = Numeric->Upos ;
+    Lpos = Numeric->Lpos ;
+    for (k = 0 ; k < n_inner ; k++)
+    {
+	CLEAR (D [k]) ;
+    }
+
+    Rs = Numeric->Rs ;
+
+    for (row = 0 ; row <= n_row ; row++)
+    {
+	Lpos [row] = EMPTY ;
+	/* Row_tuples [row] = 0 ; set in UMF_build_tuples */
+	/* Row_degree [row] = 0 ; initialized below */
+	Row_tlen [row] = 0 ;
+	/* Frpos [row] = EMPTY ;  do this later */
+    }
+
+    for (col = 0 ; col <= n_col ; col++)
+    {
+	Upos [col] = EMPTY ;
+	/* Col_tuples [col] = 0 ; set in UMF_build_tuples */
+	/* Col_degree [col] = 0 ; initialized below */
+	Col_tlen [col] = 0 ;
+	Fcpos [col] = EMPTY ;
+	Wrp [col] = 0 ;
+    }
+    Work->Wrpflag = 1 ;
+
+    /* When cleared, Wp [0..nn] is < 0 */
+    for (i = 0 ; i <= nn ; i++)
+    {
+	Wp [i] = EMPTY ;
+    }
+    /* In col search, Wp [row] is set to a position, which is >= 0. */
+
+    /* When cleared, Wrp [0..n_col] is < Wrpflag */
+    /* In row search, Wrp [col] is set to Wrpflag. */
+
+    /* no need to initialize Wm, Wio, Woi, and Woo */
+
+    /* clear the external degree counters */
+    Work->cdeg0 = 1 ;
+    Work->rdeg0 = 1 ;
+
+    fixQ = Symbolic->fixQ ;
+
+    E = Work->E ;
+
+    Numeric->n_row = n_row ;
+    Numeric->n_col = n_col ;
+    Numeric->npiv = 0 ;
+    Numeric->nnzpiv = 0 ;
+    Numeric->min_udiag = 0.0 ;
+    Numeric->max_udiag = 0.0 ;
+    Numeric->rcond = 0.0 ;
+    Numeric->isize = 0 ;
+    Numeric->nLentries = 0 ;
+    Numeric->nUentries = 0 ;
+    Numeric->lnz = 0 ;
+    Numeric->unz = 0 ;
+    Numeric->all_lnz = 0 ;
+    Numeric->all_unz = 0 ;
+    Numeric->maxfrsize = 0 ;
+    Numeric->maxnrows = 0 ;
+    Numeric->maxncols = 0 ;
+    Numeric->flops = 0. ;
+    Numeric->n1 = n1 ;
+    droptol = Numeric->droptol ;
+    drop = (droptol > 0) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* compute the scale factors, if requested, and check the input matrix */
+    /* ---------------------------------------------------------------------- */
+
+    /* UMFPACK_SCALE_SUM: Rs [i] = sum of the absolute values in row i.
+     * UMFPACK_SCALE_MAX: Rs [i] = max of the absolute values in row i.
+     *
+     * If A is complex, an approximate abs is used (|xreal| + |ximag|).
+     *
+     * If min (Rs [0..n_row]) >= RECIPROCAL_TOLERANCE, then the scale
+     * factors are inverted, and the rows of A are multiplied by the scale
+     * factors.  Otherwise, the rows are divided by the scale factors.  If
+     * NRECIPROCAL is defined, then the rows are always divided by the scale
+     * factors.
+     *
+     * For MATLAB (either built-in routine or mexFunction), or for gcc,
+     * the rows are always divided by the scale factors.
+     */
+
+    do_scale = (Numeric->scale != UMFPACK_SCALE_NONE) ;
+
+    if (do_scale)
+    {
+	int do_max = Numeric->scale == UMFPACK_SCALE_MAX ;
+	for (row = 0 ; row < n_row ; row++)
+	{
+	    Rs [row] = 0.0 ;
+	}
+	for (col = 0 ; col < n_col ; col++)
+	{
+	    ilast = EMPTY ;
+	    p1 = Ap [col] ;
+	    p2 = Ap [col+1] ;
+	    if (p1 > p2)
+	    {
+		/* invalid matrix */
+		DEBUGm4 (("invalid matrix (Ap)\n")) ;
+		return (FALSE) ;
+	    }
+	    for (p = p1 ; p < p2 ; p++)
+	    {
+		Entry aij ;
+		double value ;
+		row = Ai [p] ;
+		if (row <= ilast || row >= n_row)
+		{
+		    /* invalid matrix, columns must be sorted, no duplicates */
+		    DEBUGm4 (("invalid matrix (Ai)\n")) ;
+		    return (FALSE) ;
+		}
+		ASSIGN (aij, Ax, Az, p, split) ;
+		APPROX_ABS (value, aij) ;
+		rs = Rs [row] ;
+		if (!SCALAR_IS_NAN (rs))
+		{
+		    if (SCALAR_IS_NAN (value))
+		    {
+			/* if any entry in the row is NaN, then the scale factor
+			 * is NaN too (for now) and then set to 1.0 below */
+			Rs [row] = value ;
+		    }
+		    else if (do_max)
+		    {
+			Rs [row] = MAX (rs, value) ;
+		    }
+		    else
+		    {
+			Rs [row] += value ;
+		    }
+		}
+		DEBUG4 (("i "ID" j "ID" value %g,  Rs[i]: %g\n",
+		    row, col, value, Rs[row])) ;
+		ilast = row ;
+	    }
+	}
+	DEBUG2 (("Rs[0] = %30.20e\n", Rs [0])) ;
+	for (row = 0 ; row < n_row ; row++)
+	{
+	    rs = Rs [row] ;
+	    if (SCALAR_IS_ZERO (rs) || SCALAR_IS_NAN (rs))
+	    {
+		/* don't scale a completely zero row, or one with NaN's */
+		Rs [row] = 1.0 ;
+	    }
+	}
+	rsmin = Rs [0] ;
+	rsmax = Rs [0] ;
+	for (row = 0 ; row < n_row ; row++)
+	{
+	    DEBUG2 (("sum %30.20e ", Rs [row])) ;
+	    rsmin = MIN (rsmin, Rs [row]) ;
+	    rsmax = MAX (rsmax, Rs [row]) ;
+	    DEBUG2 (("Rs["ID"] = %30.20e\n", row, Rs [row])) ;
+	}
+#ifndef NRECIPROCAL
+	/* multiply by the reciprocal if Rs is not too small */
+	do_recip = (rsmin >= RECIPROCAL_TOLERANCE) ;
+	if (do_recip)
+	{
+	    /* invert the scale factors */
+	    for (row = 0 ; row < n_row ; row++)
+	    {
+		Rs [row] = 1.0 / Rs [row] ;
+	    }
+	}
+#endif
+    }
+    else
+    {
+	/* no scaling, rsmin and rsmax not computed */
+	rsmin = -1 ;
+	rsmax = -1 ;
+#ifndef NRECIPROCAL
+	do_recip = FALSE ;
+#endif
+	/* check the input matrix */
+	if (!AMD_valid (n_row, n_col, Ap, Ai))
+	{
+	    /* matrix is invalid */
+	    return (FALSE) ;
+	}
+    }
+
+    Numeric->rsmin = rsmin ;
+    Numeric->rsmax = rsmax ;
+#ifndef NRECIPROCAL
+    Numeric->do_recip = do_recip ;
+#else
+    Numeric->do_recip = FALSE ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* construct the inverse row Rperm_init permutation (use Frpos as temp) */
+    /* ---------------------------------------------------------------------- */
+
+    DEBUG3 (("\n\n===================LOAD_MATRIX:\n")) ;
+
+    for (newrow = 0 ; newrow < n_row ; newrow++)
+    {
+	oldrow = Rperm_init [newrow] ;
+	ASSERT (oldrow >= 0 && oldrow < n_row) ;
+	Frpos [oldrow] = newrow ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* construct the diagonal imap if doing symmetric pivoting */
+    /* ---------------------------------------------------------------------- */
+
+    if (prefer_diagonal)
+    {
+	ASSERT (n_row == n_col) ;
+	ASSERT (nempty_col == Symbolic->nempty_row) ;
+	ASSERT (nempty_col == nempty) ;
+	for (i = 0 ; i < nn ; i++)
+	{
+	    Diagonal_map [i] = EMPTY ;
+	    Diagonal_imap [i] = EMPTY ;
+	}
+	for (k = n1 ; k < nn - nempty ; k++)
+	{
+	    newrow = Symbolic->Diagonal_map [k] ;
+	    Diagonal_map [k] = newrow ;
+	    Diagonal_imap [newrow] = k ;
+	}
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate O (n_row) workspace at the tail end of Memory */
+    /* ---------------------------------------------------------------------- */
+
+    rpi = UMF_mem_alloc_tail_block (Numeric, UNITS (Int *, n_row+1)) ;
+    rpx = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry *, n_row+1)) ;
+    if (!rpi || !rpx)
+    {
+	/* :: pattern change (out of memory for Rpx, Rpx) :: */
+	/* out of memory, which can only mean that the pattern has changed */
+	return (FALSE) ;	/* pattern changed */
+    }
+    Rpi = (Int   **) (Memory + rpx) ;
+    Rpx = (Entry **) (Memory + rpi) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate the LU factors for the columns of the singletons */
+    /* ---------------------------------------------------------------------- */
+
+    DEBUG1 (("Allocating singletons:\n")) ;
+    for (k = 0 ; k < n1 ; k++)
+    {
+	lnz = Cdeg [k] - 1 ;
+	unz = Rdeg [k] - 1 ;
+
+	DEBUG1 (("Singleton k "ID" pivrow "ID" pivcol "ID" cdeg "ID" rdeg "
+	    ID"\n", k, Rperm_init [k], Cperm_init [k], Cdeg [k], Rdeg [k])) ;
+	ASSERT (unz >= 0 && lnz >= 0 && (lnz == 0 || unz == 0)) ;
+	DEBUG1 (("   lnz "ID" unz "ID"\n", lnz, unz)) ;
+
+	size = UNITS (Int, lnz) + UNITS (Entry, lnz)
+	     + UNITS (Int, unz) + UNITS (Entry, unz) ;
+	p = UMF_mem_alloc_head_block (Numeric, size) ;
+	DEBUG1 (("Kernel init head usage: "ID"\n", Numeric->ihead)) ;
+	if (!p)
+	{
+	    /* :: pattern change (out of memory for singletons) :: */
+	    DEBUG0 (("Pattern has gotten larger - kernel init failed\n")) ;
+	    return (FALSE) ;	/* pattern changed */
+	}
+
+	Numeric->all_lnz += lnz ;
+	Numeric->all_unz += unz ;
+
+	/* allocate the column of L */
+	lip = p ;
+	p += UNITS (Int, lnz) ;
+	p += UNITS (Entry, lnz) ;
+
+	/* allocate the row of U */
+	uip = p ;
+	Rpi [k] = (Int *) (Memory + p) ;
+	p += UNITS (Int, unz) ;
+	Rpx [k] = (Entry *) (Memory + p) ;
+	/* p += UNITS (Entry, unz) ; (not needed) */
+
+	/* a single column of L (no Lchains) */
+	Lip [k] = lip ;
+	Lilen [k] = lnz ;
+
+	/* a single row of L (no Uchains) */
+	Uip [k] = uip ;
+	Uilen [k] = unz ;
+
+	Wp [k] = unz ;
+
+	/* save row and column inverse permutation */
+	k1 = ONES_COMPLEMENT (k) ;
+	Rperm [k] = k1 ;			/* aliased with Row_degree */
+	Cperm [k] = k1 ;			/* aliased with Col_degree */
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* current frontal matrix is empty */
+    /* ---------------------------------------------------------------------- */
+
+    e = 0 ;
+    E [e] = 0 ;
+    Work->Flublock = (Entry *) NULL ;
+    Work->Flblock  = (Entry *) NULL ;
+    Work->Fublock  = (Entry *) NULL ;
+    Work->Fcblock  = (Entry *) NULL ;
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate the column elements */
+    /* ---------------------------------------------------------------------- */
+
+    Esize = Symbolic->Esize ;
+    empty_elements = FALSE  ;
+    for (k = n1 ; k < n_col - nempty_col ; k++)
+    {
+	e = k - n1 + 1 ;
+	ASSERT (e < Work->elen) ;
+	esize = Esize ? Esize [k-n1] : Cdeg [k] ;
+	if (esize > 0)
+	{
+	    /* allocate an element for this column */
+	    E [e] = UMF_mem_alloc_element (Numeric, esize, 1, &Rows, &Cols, &C,
+		&size, &ep) ;
+	    if (E [e] <= 0)
+	    {
+		/* :: pattern change (out of memory for column elements) :: */
+		return (FALSE) ;	/* pattern has changed */
+	    }
+	    Cols [0] = k ;
+	    DEBUG0 (("Got column element e "ID" esize "ID"\n", e, esize)) ;
+	}
+	else
+	{
+	    /* all rows in this column are dense, or empty */
+	    E [e] = 0 ;
+	    empty_elements = TRUE  ;
+	    DEBUG0 (("column element e is empty "ID"\n", e)) ;
+	}
+    }
+    DEBUG0 (("e "ID" n_col "ID" nempty_col "ID" n1 "ID"\n", e, n_col,
+		nempty_col, n1)) ;
+    ASSERT (e == n_col - nempty_col - n1) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate the row elements for dense rows of A (if any) */
+    /* ---------------------------------------------------------------------- */
+
+    if (Esize)
+    {
+	for (k = n1 ; k < n_row - nempty_row ; k++)
+	{
+	    rdeg = Rdeg [k] ;
+	    if (rdeg > dense_row_threshold)
+	    {
+		/* allocate an element for this dense row */
+		e++ ;
+		ASSERT (e < Work->elen) ;
+		E [e] = UMF_mem_alloc_element (Numeric, 1, rdeg, &Rows, &Cols,
+		    &C, &size, &ep) ;
+		if (E [e] <= 0)
+		{
+		    /* :: pattern change (out of memory for row elements) :: */
+		    return (FALSE) ;	/* pattern has changed */
+		}
+		Rows [0] = k ;
+		Rpi [k] = Cols ;
+		Rpx [k] = C ;
+		Wp [k] = rdeg ;
+		DEBUG0 (("Got row element e "ID" rdeg "ID"\n", e, rdeg)) ;
+	    }
+	}
+    }
+
+    /* elements are currently in the range 0 to e */
+    Work->nel = e ;
+
+    /* ---------------------------------------------------------------------- */
+    /* create the first n1 columns of L and U */
+    /* ---------------------------------------------------------------------- */
+
+    for (k = 0 ; k < n1 ; k++)
+    {
+	pivcol = Cperm_init [k] ;
+	p2 = Ap [pivcol+1] ;
+
+	/* get the kth column of L */
+	p = Lip [k] ;
+	Li = (Int *) (Memory + p) ;
+	lilen = Lilen [k] ;
+	p += UNITS (Int, lilen) ;
+	Lval = (Entry *) (Memory + p) ;
+
+	llen = 0 ;
+
+	for (pa = Ap [pivcol] ; pa < p2 ; pa++)
+	{
+	    oldrow = Ai [pa] ;
+	    newrow = Frpos [oldrow] ;
+	    ASSIGN (x, Ax, Az, pa, split) ;
+
+	    /* scale the value using the scale factors, Rs */
+	    if (do_scale)
+	    {
+#ifndef NRECIPROCAL
+		if (do_recip)
+		{
+		    SCALE (x, Rs [oldrow]) ;
+		}
+		else
+#endif
+		{
+		    SCALE_DIV (x, Rs [oldrow]) ;
+		}
+	    }
+
+	    if (newrow == k)
+	    {
+		/* this is the pivot entry itself */
+		ASSERT (oldrow == Rperm_init [k]) ;
+		D [k] = x ;
+	    }
+	    else if (newrow < k)
+	    {
+		/* this entry goes in a row of U */
+		DEBUG1 (("Singleton row of U: k "ID" newrow "ID"\n",
+		    k, newrow)) ;
+		if (--(Wp [newrow]) < 0)
+		{
+		    /* :: pattern change (singleton row too long) :: */
+		    DEBUGm4 (("bad U singleton row (too long)\n")) ;
+		    return (FALSE) ;	/* pattern changed */
+		}
+		*(Rpi [newrow]++) = k ;
+		*(Rpx [newrow]++) = x ;
+	    }
+	    else
+	    {
+		/* this entry goes in a column of L */
+		DEBUG1 (("Singleton col of L: k "ID" newrow "ID"\n",
+		    k, newrow)) ;
+		if (llen >= lilen)
+		{
+		    DEBUGm4 (("bad L singleton col (too long)\n")) ;
+		    return (FALSE) ;	/* pattern changed */
+		}
+		Li   [llen] = newrow ;
+		Lval [llen] = x ;
+		llen++ ;
+	    }
+	}
+
+	if (llen != lilen)
+	{
+	    /* :: pattern change (singleton column too long) :: */
+	    DEBUGm4 (("bad L singleton col (too short)\n")) ;
+	    return (FALSE) ;	/* pattern changed */
+	}
+
+	/* scale the column of L */
+	if (llen > 0)
+	{
+	    pivot_value = D [k] ;
+	    UMF_scale (llen, pivot_value, Lval) ;
+	}
+
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate the elements and copy the columns of A */
+    /* ---------------------------------------------------------------------- */
+
+    /* also apply the row and column pre-ordering.  */
+    for (k = n1 ; k < n_col ; k++)
+    {
+	/* The newcol is k, which is what the name of the column is in the
+	 * UMFPACK kernel.  The user's name for the column is oldcol. */
+	oldcol = Cperm_init [k] ;
+
+	ASSERT (oldcol >= 0 && oldcol < n_col) ;
+
+	p2 = Ap [oldcol+1] ;
+
+	cdeg = Cdeg [k] ;
+	ASSERT (cdeg >= 0) ;
+	ASSERT (IMPLIES (
+	    (Symbolic->ordering != UMFPACK_ORDERING_GIVEN) && n1 > 0,
+	    cdeg > 1 || cdeg == 0)) ;
+
+	/* if fixQ: set Col_degree to 0 for the NON_PIVOTAL_COL macro */
+	Col_degree [k] = fixQ ? 0 : cdeg ;
+
+	/* get the element for this column (if any) */
+	e = k - n1 + 1 ;
+	if (k < n_col - nempty_col)
+	{
+	    esize = Esize ? Esize [k-n1] : cdeg ;
+	    if (E [e])
+	    {
+		Int ncols, nrows ;
+		Unit *pp ;
+		pp = Memory + E [e] ;
+		GET_ELEMENT (ep, pp, Cols, Rows, ncols, nrows, C) ;
+		ASSERT (ncols == 1) ;
+		ASSERT (nrows == esize) ;
+		ASSERT (Cols [0] == k) ;
+	    }
+	}
+	else
+	{
+	    ASSERT (cdeg == 0) ;
+	    esize = 0 ;
+	}
+
+	clen = 0 ;
+
+	for (pa = Ap [oldcol] ; pa < p2 ; pa++)
+	{
+	    oldrow = Ai [pa] ;
+	    newrow = Frpos [oldrow] ;
+	    ASSIGN (x, Ax, Az, pa, split) ;
+
+	    /* scale the value using the scale factors, Rs */
+	    if (do_scale)
+	    {
+#ifndef NRECIPROCAL
+		if (do_recip)
+		{
+		    /* multiply by the reciprocal */
+		    SCALE (x, Rs [oldrow]) ;
+		}
+		else
+#endif
+		{
+		    /* divide instead */
+		    SCALE_DIV (x, Rs [oldrow]) ;
+		}
+	    }
+
+	    rdeg = Rdeg [newrow] ;
+	    if (newrow < n1 || rdeg > dense_row_threshold)
+	    {
+		/* this entry goes in a row of U or into a dense row */
+		DEBUG1 (("Singleton/dense row of U: k "ID" newrow "ID"\n",
+		    k, newrow)) ;
+		if (--(Wp [newrow]) < 0)
+		{
+		    DEBUGm4 (("bad row of U or A (too long)\n")) ;
+		    return (FALSE) ;	/* pattern changed */
+		}
+		*(Rpi [newrow]++) = k ;
+		*(Rpx [newrow]++) = x ;
+	    }
+	    else
+	    {
+		/* this entry goes in an initial element */
+		DEBUG1 (("In element k "ID" e "ID" newrow "ID"\n",
+		    k, e, newrow)) ;
+		if (clen >= esize)
+		{
+		    DEBUGm4 (("bad A column (too long)\n")) ;
+		    return (FALSE) ;	/* pattern changed */
+		}
+		ASSERT (E [e]) ;
+		ASSERT (k < n_col - nempty_col) ;
+		Rows [clen] = newrow ;
+		C    [clen] = x ;
+		clen++ ;
+#ifndef NDEBUG
+		if (Diagonal_map && (newrow == Diagonal_map [k]))
+		{
+		    DEBUG0 (("Diagonal: old: row "ID" col "ID" : "
+			"new: row "ID" col "ID" : ",
+			oldrow, oldcol, newrow, k)) ;
+		    EDEBUGk (0, x) ;
+		}
+#endif
+	    }
+	}
+
+	if (clen != esize)
+	{
+	    /* :: pattern change (singleton column too short) :: */
+	    DEBUGm4 (("bad A column (too short)\n")) ;
+	    return (FALSE) ;	/* pattern changed */
+	}
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* free the Rpi and Rpx workspace at the tail end of memory */
+    /* ---------------------------------------------------------------------- */
+
+    UMF_mem_free_tail_block (Numeric, rpi) ;
+    UMF_mem_free_tail_block (Numeric, rpx) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* prune zeros and small entries from the singleton rows and columns */
+    /* ---------------------------------------------------------------------- */
+
+    if (n1 > 0)
+    {
+	pnew = Lip [0] ;
+	ASSERT (pnew == 1) ;
+	for (k = 0 ; k < n1 ; k++)
+	{
+	    DEBUGm4 (("\nPrune singleton L col "ID"\n", k)) ;
+	    pnew = packsp (pnew, &Lip [k], &Lilen [k], drop, droptol, Memory) ;
+	    Numeric->lnz += Lilen [k] ;
+	    DEBUGm4 (("\nPrune singleton U row "ID"\n", k)) ;
+	    pnew = packsp (pnew, &Uip [k], &Uilen [k], drop, droptol, Memory) ;
+	    Numeric->unz += Uilen [k] ;
+	}
+	/* free the unused space at the head of memory */
+	Numeric->ihead = pnew ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* initialize row degrees */
+    /* ---------------------------------------------------------------------- */
+
+    for (k = 0 ; k < n1 ; k++)
+    {
+	if (Wp [k] != 0)
+	{
+	    /* :: pattern change (singleton row too short) :: */
+	    DEBUGm4 (("bad U singleton row (too short)\n")) ;
+	    return (FALSE) ;	/* pattern changed */
+	}
+    }
+
+    for (k = n1 ; k < n_row ; k++)
+    {
+	DEBUG1 (("Initial row degree k "ID" oldrow "ID" Rdeg "ID"\n",
+	    k, Rperm_init [k], Rdeg [k])) ;
+	rdeg = Rdeg [k] ;
+	Row_degree [k] = rdeg ;
+	if (rdeg > dense_row_threshold && Wp [k] != 0)
+	{
+	    /* :: pattern change (dense row too short) :: */
+	    DEBUGm4 (("bad dense row (too short)\n")) ;
+	    return (FALSE) ;	/* pattern changed */
+	}
+    }
+
+#ifndef NDEBUG
+    if (prefer_diagonal)
+    {
+	Entry aij ;
+	Int *InvCperm, newcol ;
+	UMF_dump_diagonal_map (Diagonal_map, Diagonal_imap, n1, nn, nempty) ;
+	InvCperm = (Int *) malloc (n_col * sizeof (Int)) ;
+	ASSERT (InvCperm != (Int *) NULL) ;
+	for (newcol = 0 ; newcol < n_col ; newcol++)
+	{
+	    oldcol = Cperm_init [newcol] ;
+	    InvCperm [oldcol] = newcol ;
+	}
+	DEBUGm3 (("Diagonal of P2*A:\n")) ;
+	for (oldcol = 0 ; oldcol < n_col ; oldcol++)
+	{
+	    newcol = InvCperm [oldcol] ;
+	    for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
+	    {
+		oldrow = Ai [p] ;
+		newrow = Frpos [oldrow] ;
+		ASSIGN (aij, Ax, Az, p, split) ;
+		if (newrow == Diagonal_map [newcol])
+		{
+		    DEBUG0 (("old row "ID" col "ID" new row "ID" col "ID,
+			oldrow, oldcol, newrow, newcol)) ;
+		    EDEBUGk (0, aij) ;
+		    DEBUG0 ((" scaled ")) ;
+		    if (do_scale)
+		    {
+#ifndef NRECIPROCAL
+			if (do_recip)
+			{
+			    SCALE (aij, Rs [oldrow]) ;
+			}
+			else
+#endif
+			{
+			    SCALE_DIV (aij, Rs [oldrow]) ;
+			}
+		    }
+		    EDEBUGk (0, aij) ;
+		    DEBUG0 (("\n")) ;
+		}
+	    }
+	}
+	free (InvCperm) ;
+    }
+#endif
+
+    Col_degree [n_col] = 0 ;
+
+    /* ---------------------------------------------------------------------- */
+    /* pack the element name space */
+    /* ---------------------------------------------------------------------- */
+
+    if (empty_elements)
+    {
+	Int e2 = 0 ;
+	DEBUG0 (("\n\n============= Packing element space\n")) ;
+	for (e = 1 ; e <= Work->nel ; e++)
+	{
+	    if (E [e])
+	    {
+		e2++ ;
+		E [e2] = E [e] ;
+	    }
+	}
+	Work->nel = e2 ;
+    }
+
+#ifndef NDEBUG
+    DEBUG0 (("Number of initial elements: "ID"\n", Work->nel)) ;
+    for (e = 0 ; e <= Work->nel ; e++) UMF_dump_element (Numeric, Work,e,TRUE) ;
+#endif
+
+    for (e = Work->nel + 1 ; e < Work->elen ; e++)
+    {
+	E [e] = 0 ;
+    }
+
+    /* Frpos no longer needed */
+    for (row = 0 ; row <= n_row ; row++)
+    {
+	Frpos [row] = EMPTY ;
+    }
+
+    /* clear Wp */
+    for (i = 0 ; i <= nn ; i++)
+    {
+	Wp [i] = EMPTY ;
+    }
+
+    DEBUG1 (("Kernel init head usage: "ID"\n", Numeric->ihead)) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* build the tuple lists */
+    /* ---------------------------------------------------------------------- */
+
+    /* if the memory usage changes, then the pattern has changed */
+
+    (void) UMF_tuple_lengths (Numeric, Work, &unused) ;
+    if (!UMF_build_tuples (Numeric, Work))
+    {
+	/* :: pattern change (out of memory in umf_build_tuples) :: */
+	/* We ran out of memory, which can only mean that */
+	/* the pattern (Ap and or Ai) has changed (gotten larger). */
+	DEBUG0 (("Pattern has gotten larger - build tuples failed\n")) ;
+	return (FALSE) ;	/* pattern changed */
+    }
+
+    Numeric->init_usage = Numeric->max_usage ;
+
+    /* ---------------------------------------------------------------------- */
+    /* construct the row merge sets */
+    /* ---------------------------------------------------------------------- */
+
+    for (i = 0 ; i <= Symbolic->nfr ; i++)
+    {
+	Work->Front_new1strow [i] = Symbolic->Front_1strow [i] ;
+    }
+
+#ifndef NDEBUG
+    UMF_dump_rowmerge (Numeric, Symbolic, Work) ;
+    DEBUG6 (("Column form of original matrix:\n")) ;
+    UMF_dump_col_matrix (Ax,
+#ifdef COMPLEX
+	Az,
+#endif
+	Ai, Ap, n_row, n_col, nz) ;
+    UMF_dump_memory (Numeric) ;
+    UMF_dump_matrix (Numeric, Work, FALSE) ;
+    DEBUG0 (("kernel init done...\n")) ;
+#endif
+
+    return (TRUE) ;
+}