The R Project SVN R-packages

/* ========================================================================== */
/* === UMF_kernel_init ====================================================== */
/* ========================================================================== */

/* -------------------------------------------------------------------------- */
/* UMFPACK Version 4.5, 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) ;
}